diff --git a/matlab/DiocotronStability.m b/matlab/DiocotronStability.m index fcd054c..b936cbe 100644 --- a/matlab/DiocotronStability.m +++ b/matlab/DiocotronStability.m @@ -1,41 +1,45 @@ rbm=0.005; rbp=0.0128; -%rbm=0.09; -%rbp=0.113; -%b=0.16; -%a=0.001; +rbm=0.09; +rbp=0.011; +rbm=0.09; +rbp=0.11; +% b=0.16; +% a=0.001; a=M.rgrid(1); b=M.rgrid(end); Q=mean(M.Er(1,floor(length(M.zgrid)/2),end-30:end),3)/2*a; %Q=mean(M.Er(1,100,end-30:end),3)/2*a; -omega_d=M.omepe^2/(2*M.omece); +Q=-(M.potout-M.potinn)*M.phinorm/(2*log(a/b)); +omepe=sqrt(9.92e14*M.qe^2/M.eps_0/M.me); +omega_d=omepe^2/(2*M.omece); if( a~=0) - omega_q=2*Q/(M.B0*(rbm^2)); + omega_q=-2*Q/(M.B0*(rbm^2)); else omega_q=0; end Levaluated=10; stability=zeros(Levaluated,1); for l=1:Levaluated stability(l)=(-l*(1-omega_q/omega_d)*(1-(rbm/rbp)^2)*(1-(a/b)^(2*l))+ ... 2*(1+(a/b)^(2*l))-(1+(rbm/rbp)^(2*l))*((a/rbm)^(2*l)+(rbp/b)^(2*l)))^2 -... 4*(rbm/rbp)^(2*l)*(1-(rbp/b)^(2*l))^2*(1-(a/rbm)^(2*l))^2; if(stability(l) <0) fprintf(2,"Diocotron mode l=%d is potentially unstable\n",l); end end fig=figure; plot(1:Levaluated,stability) hold on plot([1 Levaluated],[0 0],'k') for l=1:Levaluated if stability(l)<0 plot(l,stability(l),'rx') end end xlabel('azimuthal mode number') ylabel('stability condition') title(sprintf('r_b^-=%1.2g [m] r_b^+=%1.2g [m] a=%1.2g [m] b=%1.2g [m]',rbm,rbp,a,b)) M.savegraph(fig,sprintf('%s/%s_diocstab',M.folder,M.name),[15,10]) \ No newline at end of file diff --git a/matlab/Testflux.m b/matlab/Testflux.m index 945861b..ae479ef 100644 --- a/matlab/Testflux.m +++ b/matlab/Testflux.m @@ -1,101 +1,139 @@ -%% -rAthetpos=30; +%% Computes the electric potential profile along the magnetic field line passing at z=0 and r=rgrid(rAthetpos) -fieldstep=length(M.t2d);%-100:5:length(M.t2d); +rAthetpos=15+27; +rAthetpos=29; + +fieldstep=length(M.t2d)-500:10:length(M.t2d); Vpar=M.fluidUR(:,:,fieldstep).*(M.Br./M.B)' + (M.Bz./M.B)'.*M.fluidUZ(:,:,fieldstep); N=M.N(:,:,fieldstep); Gammapar=mean(N.*Vpar,3); - levels=M.rAthet(rAthetpos,floor(M.nz/2)+1)*[1 1]; rAthetposend=find(levels(1)>=M.rAthet(:,1),1,'last'); rleft=sqrt(levels(1)/(M.rAthet(rAthetposend,1)/M.rgrid(rAthetposend)^2)); -thefig=figure; + +%% triple subplot showing the magnetic field line of interest on the: +% particle density, parallel electron flux and parallel velocity surface plots +thefig=figure('Name', sprintf('%s flux line ',M.name)); sp(1)=subplot(1,3,1); contourf(M.zgrid,M.rgrid,mean(N,3),'edgecolor','none'); hold on; contour(M.zgrid,M.rgrid,M.rAthet,levels,'r-','linewidth',1.5,'Displayname','Magnetic field lines'); c=colorbar; c.Label.String='n_e [m^{-3}]'; + sp(2)=subplot(1,3,2); contourf(M.zgrid,M.rgrid,Gammapar,'edgecolor','none'); hold on contour(M.zgrid,M.rgrid,M.rAthet,levels,'r-','linewidth',1.5,'Displayname','Magnetic field lines'); c=colorbar; c.Label.String='\Gamma_{par} [m^{-2}s^{-1}]'; + sp(3)=subplot(1,3,3); contourf(M.zgrid,M.rgrid,mean(Vpar,3),'edgecolor','none'); hold on contour(M.zgrid,M.rgrid,M.rAthet,levels,'r-','linewidth',1.5,'Displayname','Magnetic field line'); c=colorbar; c.Label.String='V_{par} [ms^{-1}]'; linkaxes(sp); xlabel(sp,'z [m]') ylabel(sp,'r [m]') ylim(sp,[M.rgrid(1) M.rgrid(sum(M.nnr(1:2)))]) M.savegraph(thefig,sprintf('%s/%s_finalwellfluxline',M.folder,M.name),[12,10]) -%% -dN=1e13*M.weight; -dN=0; -rp=9.85e-3; -rm=7e-3; -V=2*pi*(rp^2-rm^2); -n2=mean(mean(mean(N(rAthetposend+(0:1),[1 end],:),3))); -vpar2=(dN/V/n2)^2; -actvpar2=mean(mean(abs(mean(Vpar(rAthetposend+(0:1),[1 end],:),3))))^2; +%% Plot of the magnetic field line on the density alone +thefig=figure('Name', sprintf('%s dens line ',M.name)); +sp=gca; +contourf(M.zgrid,M.rgrid,mean(N,3),'edgecolor','none'); +hold on; +contour(M.zgrid,M.rgrid,M.rAthet,levels,'r-','linewidth',1.5,'Displayname','Magnetic field lines'); +c=colorbar; +c.Label.String='n_e [m^{-3}]'; +xlabel(sp,'z [m]') +ylabel(sp,'r [m]') +ylim(sp,[M.rgrid(1) M.rgrid(sum(M.nnr(1:2)))]) +M.savegraph(thefig,sprintf('%s/%s_finalwellline',M.folder,M.name),[9,8]) + + + +%% Computation of a prediction for the maximum potential hill height in steady state, +% depending on the external electric field, the particle maxwellian source temperature and the magnetic field mirror ratio and amplitude b=M.rgrid(end); a=M.rgrid(1); vd1=((M.potout-M.potinn)*M.phinorm/M.rgrid(rAthetpos))/M.Bz(floor(M.nz/2)+1,rAthetpos)/log(b/a); %vd1=-M.Er(rAthetpos+5,M.nz/2+1,end)/M.Bz(M.nz/2+1,rAthetpos+5) vinit=sqrt(M.kb*10000/M.me) %vinit=sqrt(80*M.qe/M.me) vd2=((M.potout-M.potinn)*M.phinorm/rleft)/M.Bz(1,rAthetposend)/log(b/a); -deltaphi=-0.5*M.me/M.qe*(vpar2-(vd1+vinit)^2*(1-M.Rcurv)) -ratioparper=vpar2/((vd1+vinit)^2*(1-M.Rcurv)) - -%% +deltaphi=0.5*M.me/M.qe*((vd1+vinit)^2*(1-M.Rcurv)) +%ratioparper=vpar2/((vd1+vinit)^2*(1-M.Rcurv)) -thefig=figure; +%% Plot of the parallel density profile and the well along the magnetic field line of interest +thefig=figure('Name', sprintf('%s final well',M.name)); %fieldstep=fieldstep(end); plotaxes(1)=subplot(1,2,1,'Parent',thefig); plotaxes(2)=subplot(1,2,2,'Parent',thefig); dens=mean(M.N(:,:,fieldstep),3); model=M.potentialwellmodel(fieldstep); z=model.z; r=model.r; pot=model.pot; rathet=model.rathet; [Zmesh,Rmesh]=meshgrid(M.zgrid,M.rAthet(:,floor(M.nz/2)+1)); densplot=zeros(size(Zmesh,1),size(Zmesh,2),length(fieldstep)); for i=1:length(fieldstep) densplot(:,:,i)=griddata(Zmesh,M.rAthet,dens,Zmesh,Rmesh,'natural'); end plot(plotaxes(1),M.zgrid,mean(densplot(rAthetpos,1:end,:),3)); xlim(plotaxes(1),[M.zgrid(1) M.zgrid(end)]) xlabel(plotaxes(1),'z [m]') ylabel(plotaxes(1),'n [m^{-3}]') title(plotaxes(1),'Density') plotpot=zeros(size(Zmesh,1),size(Zmesh,2),length(fieldstep)); for i=1:length(fieldstep) plotpot(:,:,i)=griddata(z,rathet,pot(:,i),Zmesh,Rmesh,'natural'); end plotpot=mean(plotpot,3); plot(plotaxes(2),M.zgrid(1:end),plotpot(rAthetpos,1:end)); hold on plot(plotaxes(2),M.zgrid(1:end),deltaphi*ones(size(M.zgrid)),'k--','displayname','Analytical') xlim(plotaxes(2),[M.zgrid(1) M.zgrid(end)]) xlabel(plotaxes(2),'z [m]') ylabel(plotaxes(2),'depth [eV]') title(plotaxes(2),'Well') sgtitle(thefig,sprintf('r(0)=%1.2f [mm] t= %1.2f [ns]',M.rgrid(rAthetpos)*1e3,M.t2d(fieldstep(end))*1e9)) M.savegraph(thefig,sprintf('%s/%s_finalwell',M.folder,M.name),[15,10]) +%% Plot of the final electric potential well along the magnetic field line +thefig=figure('Name', sprintf('%s final well only',M.name)); +theplots=gobjects(2,1); +theplots(1)=plot(M.zgrid(1:end),-plotpot(rAthetpos,1:end),'displayname','Simulation'); +hold on +theplots(2)=plot(M.zgrid(1:end),-deltaphi*ones(size(M.zgrid)),'k--','displayname','\Delta\phi_{21} Prediction'); +xlim([M.zgrid(1) M.zgrid(end)]) +xlabel('z [m]') +ylabel('\phi(z)-\phi(0) [V]') +legend + +ylimits=ylim; +%annotation('textarrow',[M.zgrid(5) M.zgrid(1)]/(M.zgrid(end)-M.zgrid(1))+0.5,[-plotpot(rAthetpos,1)*1.05 -plotpot(rAthetpos,1)]/diff(ylimits),'Linewidth',2) +p1=[M.zgrid(30),-plotpot(rAthetpos,1)*1.02]; +p2=[M.zgrid(1),-plotpot(rAthetpos,1)]; +dp=p2-p1; +quiver(p1(1), p1(2), dp(1), dp(2),1,'k','showarrowhead','off') +text(p1(1),p1(2),'\Delta\phi_{21}') +title(sprintf('r(0)=%1.2f [mm] \\Delta\\phi_{ba}= %1.2g[kV] R=%1.1f',M.rgrid(rAthetpos)*1e3,(M.potout-M.potinn)*M.phinorm,M.Rcurv)) + +legend(theplots,'location','north') +M.savegraph(thefig,sprintf('%s/%s_finalwellonly',M.folder,M.name),[9,8]) + + + diff --git a/matlab/WellVideosaving.m b/matlab/WellVideosaving.m index 6006adc..ad3b9cf 100644 --- a/matlab/WellVideosaving.m +++ b/matlab/WellVideosaving.m @@ -1,108 +1,108 @@ fieldstep=1; -maxN=10e15; -Minwell=-1700; +maxN=5e12; +Minwell=-10; tend=min(8000,length(M.t2d)); tend=length(M.t2d); plot3d=false; rAthetpos=25; thefig=figure('Position',[0 0 1600 900]); if plot3d filename=[M.name,'_well3D.avi']; else filename=[M.name,'_well2D.avi']; end videowriterobj=VideoWriter([M.folder,'/',filename]); videowriterobj.FrameRate=10; open(videowriterobj); plotaxes(1)=subplot(1,2,1,'Parent',thefig); plotaxes(2)=subplot(1,2,2,'Parent',thefig); if plot3d sf=surface(plotaxes(1),M.zgrid,M.rgrid,M.N(:,:,fieldstep),'edgecolor','none'); xlim(plotaxes(1),[M.zgrid(1) M.zgrid(end)]) ylim(plotaxes(1),[M.rgrid(1) M.rgrid(sum(M.nnr(1:2))+5)]) xlabel(plotaxes(1),'z [m]') ylabel(plotaxes(1),'r [m]') title(plotaxes(1),'Density') c = colorbar(plotaxes(1)); c.Label.String= 'n[m^{-3}]'; caxis(plotaxes(1),[-Inf maxN]); well=M.PotentialWell(fieldstep,false); sf=surface(plotaxes(2),M.zgrid(2:end),M.rgrid(2:end),well(2:end,2:end),'edgecolor','none'); xlim(plotaxes(2),[M.zgrid(1) M.zgrid(end)]) ylim(plotaxes(2),[M.rgrid(1) M.rgrid(sum(M.nnr(1:2))+5)]) xlabel(plotaxes(2),'z [m]') ylabel(plotaxes(2),'r [m]') title(plotaxes(2),'Well') c = colorbar(plotaxes(2)); c.Label.String= 'depth [V]'; view(plotaxes(2),3) colormap(plotaxes(2),'jet'); caxis(plotaxes(2),[Minwell 0]); zlim(plotaxes(2),[Minwell 0]); for i=fieldstep:5:tend sgtitle(thefig,sprintf('t= %1.2f [ns]',M.t2d(i)*1e9)) dens=M.N(:,:,i); well=M.PotentialWell(i,false); well=well(2:end,2:end); plotaxes(1).Children.ZData=dens; plotaxes(1).Children.CData=dens; plotaxes(2).Children.ZData=well; plotaxes(2).Children.CData=well; writeVideo(videowriterobj,getframe(thefig)); end else dens=M.N(:,:,fieldstep); model=M.potentialwellmodel(fieldstep); z=model.z; r=model.r; pot=model.pot; rathet=model.rathet; [Zmesh,Rmesh]=meshgrid(M.zgrid,M.rAthet(:,1)); dens=griddata(Zmesh,M.rAthet,dens,Zmesh,Rmesh); plot(plotaxes(1),M.zgrid,dens(rAthetpos,1:end)); xlim(plotaxes(1),[M.zgrid(1) M.zgrid(end)]) ylim(plotaxes(1),[0 maxN]); xlabel(plotaxes(1),'z [m]') ylabel(plotaxes(1),'n [m^{-3}]') title(plotaxes(1),'Density') grid(plotaxes, 'on') pot=griddata(z,rathet,pot,Zmesh,Rmesh); plot(plotaxes(2),M.zgrid(1:end),pot(rAthetpos,1:end)); xlim(plotaxes(2),[M.zgrid(1) M.zgrid(end)]) ylim(plotaxes(2),[Minwell 0]); xlabel(plotaxes(2),'z [m]') ylabel(plotaxes(2),'depth [eV]') title(plotaxes(2),'Well') grid on for i=fieldstep:20:tend sgtitle(thefig,sprintf('rA_\\theta=%1.3g [Tm^2] t= %1.2f [ns]',Rmesh(rAthetpos,1),M.t2d(i)*1e9)) dens=M.N(:,:,i); model=M.potentialwellmodel(i); z=model.z; r=model.r; pot=model.pot; rathet=model.rathet; [Zmesh,Rmesh]=meshgrid(M.zgrid,M.rAthet(:,1)); dens=griddata(Zmesh,M.rAthet,dens,Zmesh,Rmesh); pot=griddata(z,rathet,pot,Zmesh,Rmesh); plotaxes(1).Children.YData=dens(rAthetpos,1:end); plotaxes(2).Children.YData=pot(rAthetpos,1:end); writeVideo(videowriterobj,getframe(thefig)); end end close(videowriterobj); \ No newline at end of file diff --git a/matlab/dispespicFields.m b/matlab/dispespicFields.m index 22ccada..61c9dbf 100644 --- a/matlab/dispespicFields.m +++ b/matlab/dispespicFields.m @@ -1,195 +1,195 @@ function dispespicFields(M,logdensity,showgrid,fixed) %dispespicFields Allows to display the time evolution of the density, electric potential and electric fields % M is of class espic2dhdf5 and contains the simulation results fieldstep=1; if nargin <2 logdensity=false; showgrid=false; fixed=false; end if nargin <3 showgrid=false; fixed=false; end if nargin <4 fixed=false; end fixed=fi(fixed); f=uifigure('Name',sprintf('Grid data %s',M.name)); mf=uipanel(f,'Position',[5 50 f.Position(3)-10 f.Position(4)-55]); mf.AutoResizeChildren='off'; m=uipanel(f,'Position',[5 5 f.Position(3)-10 40]); sgtitle(mf,sprintf('step=%d t=%0.5e s',fieldstep*M.it1,M.t2d(fieldstep))) sld = uislider(m,'Position',[10 30 0.6*m.Position(3) 3]); sld.Value=fieldstep; sld.Limits=[1 size(M.t2d,1)]; edt = uieditfield(m,'numeric','Limits',[1 size(M.t2d,1)],'Value',1); edt.Position=[sld.Position(1)+sld.Position(3)+25 5 40 20]; edt.RoundFractionalValues='on'; MaxN=0; Printbt=uibutton(m,'Position',[edt.Position(1)+edt.Position(3)+10 5 40 20],'Text', 'Save'); Play=uibutton(m,'Position',[Printbt.Position(1)+Printbt.Position(3)+10 5 40 20],'Text', 'Play'); Pause=uibutton(m,'Position',[Play.Position(1)+Play.Position(3)+10 5 40 20],'Text', 'Pause'); %Playbt=uibutton(m,'Position',[Printbt.Position(1)+Printbt.Position(3)+10 5 40 20],'Text', 'Play/Pause'); stop=false; sld.ValueChangingFcn={@updatefigdata,edt,mf}; edt.ValueChangedFcn={@updatefigdata,sld,mf}; Printbt.ButtonPushedFcn={@plotGridButtonPushed}; Play.ButtonPushedFcn={@plotPlayButtonPushed}; Pause.ButtonPushedFcn={@PauseButtonPushed}; PlotEspic2dgriddata(mf,M,fieldstep); function plotPlayButtonPushed(btn,ax) stop=false; i=sld.Value; while ~stop edt.Value=i; sld.Value=i; updatesubplotsdata(i,mf); pause(0.01) i=sld.Value; i=i+10; if(i>sld.Limits(2)) stop=true; end end end function PauseButtonPushed(btn,ax) stop = true; end function PlotEspic2dgriddata(fig,M,fieldstep) %PlotEspic2dgriddata Plot the 2d data of espic2d at time step fieldstep sgtitle(fig,sprintf('step=%d t=%0.5e s',(fieldstep-1)*M.it1,M.t2d(fieldstep))) ax1=subplot(2,2,1,'Parent',fig); sf=surface(ax1,M.zgrid,M.rgrid,M.N(:,:,fieldstep),'edgecolor','none'); if(showgrid) set(sf,'edgecolor','black') end xlim(ax1,[M.zgrid(1) M.zgrid(end)]) ylim(ax1,[M.rgrid(1) M.rgrid(end)]) xlabel(ax1,'z [m]') ylabel(ax1,'r [m]') title(ax1,'Position') c = colorbar(ax1); c.Label.String= 'n[m^{-3}]'; %c.Limits=[0 max(M.N(:))]; if(isboolean(fixed) && fixed) climits=caxis(ax1); MaxN=climits(2); elseif(~isboolean(fixed)) MaxN=fixed.data; caxis(ax1,[-Inf MaxN]); end view(ax1,2) if logdensity set(ax1,'zscale','log') set(ax1,'colorscale','log') end ax2=subplot(2,2,2,'Parent',fig); - contourf(ax2,M.zgrid,M.rgrid,M.pot(:,:,fieldstep)); + surface(ax2,M.zgrid,M.rgrid,M.pot(:,:,fieldstep),'edgecolor','none'); %plot(ax2,M.zgrid,data.pot(:,5)) xlabel(ax2,'z [m]') ylabel(ax2,'r [m]') colormap(ax2,'jet') c = colorbar(ax2); c.Label.String= '\Phi [V]'; title(ax2,'Es potential') grid(ax2, 'on') ax3=subplot(2,2,3,'Parent',fig); contourf(ax3,M.zgrid,M.rgrid,M.Epar(fieldstep)) xlabel(ax3,'z [m]') ylabel(ax3,'r [m]') c = colorbar(ax3); c.Label.String= 'E_{par} [V/m]'; title(ax3,'Parallel Electric field') grid(ax3, 'on') ax4=subplot(2,2,4,'Parent',fig); contourf(ax4,M.zgrid,M.rgrid,M.Eperp(fieldstep)) xlabel(ax4,'z [m]') ylabel(ax4,'r [m]') c = colorbar(ax4); c.Label.String= 'E_{per} [V/m]'; title(ax4,'Perpendicular Electric field') grid(ax4, 'on') linkaxes([ax1,ax2,ax3,ax4]); end function plotGridButtonPushed(btn,ax) %UNTITLED2 Summary of this function goes here % Detailed explanation goes here f=figure(); PlotEspic2dgriddata(f,M,sld.Value); f.PaperOrientation='landscape'; [~, name, ~] = fileparts(M.file); print(f,sprintf('%sGrid%d',name,sld.Value),'-dpdf','-fillpage') end function updatefigdata(control, event, Othercontrol, fig) if strcmp(event.EventName,'ValueChanged') fieldstep=floor(control.Value); control.Value=fieldstep; else fieldstep=floor(event.Value); end Othercontrol.Value=fieldstep; updatesubplotsdata(fieldstep, fig); end function updatesubplotsdata(fieldstep, fig) sgtitle(fig,sprintf('step=%d t=%0.5e s',(fieldstep-1)*M.it1,M.t2d(floor(fieldstep)))) %% update Position histogram ax1=fig.Children(end); dens=M.N(:,:,fieldstep); ax1.Children.ZData=dens; ax1.Children.CData=dens; if(isboolean(fixed) && fixed) climits=caxis(ax1); MaxN=climits(2); nmax=max(dens(:)); if(nmax>MaxN) MaxN=nmax; end caxis(ax1,[0 MaxN]); elseif(~isboolean(fixed)) MaxN=fixed.data; caxis(ax1,[-Inf MaxN]); end % view(ax1,2) %% update ES Potential fig.Children(end-2).Children(1).ZData=M.pot(:,:,fieldstep); %% update Radial electric field fig.Children(end-4).Children(1).ZData=M.Epar(fieldstep); %% update Axial electric field fig.Children(end-6).Children(1).ZData=M.Eperp(fieldstep); drawnow limitrate end end diff --git a/matlab/distribution_function.m b/matlab/distribution_function.m index 5921055..cc91e1e 100644 --- a/matlab/distribution_function.m +++ b/matlab/distribution_function.m @@ -1,357 +1,339 @@ %% Show the particles velocity histogram at the position set using ginput % The histogram is compared between timesteppart and end -timesteppart=length(M.tpart); +timeinit=1; +timesteppart=length(M.tpart)-1; +%timesteppart=1; +timestepNinit=find(M.tpart(timeinit)==M.t2d); +timestepNend=find(M.tpart(timesteppart)==M.t2d); + +gcs=false; % use guiding center reference system for Vperp coordinates + +zleftlim=1; + +%Rindex=15+28; +Rindex=28;%43;%29; +Zindex=193; fig=figure; ax1=gca; +% get the relevant timesteps between 2d and particles variables timestepN=find(M.tpart(end)==M.t2d,1); +% Plot the density at the end of simulation to select the region of interest Ndistrib=M.N(:,:,timestepN); h=contourf(ax1,M.zgrid,M.rgrid,Ndistrib); %set(h, 'EdgeColor', 'none'); hold on [r,z]=find(Ndistrib~=0); xlim(ax1,[M.zgrid(min(z)) M.zgrid(max(z))]) ylim(ax1,[M.rgrid(min(r)) M.rgrid(max(r))]) xlabel(ax1,'Z [m]') ylabel(ax1,'R [m]') c = colorbar(ax1); c.Label.String= 'n [m^{-3}]'; view(ax1,2) -%set(ax1,'colorscale','log') -% [x,y]=ginput(1); -% Zindex=find(x>M.zgrid,1,'last'); -% Rindex=find(y>M.rgrid,1,'last'); -Rindex=29 -Zindex=193; +%% Define the points of interest for ploting the velocity distribution +[x,y]=ginput(1); +[~,Zindex]=min(abs(x-M.zgrid)); +[~,Rindex]=min(abs(y-M.rgrid)); Z=(M.zgrid(Zindex)); R=(M.rgrid(Rindex)); + +% Show this point on the density map plot(Z,R,'rx','Markersize',12); % Rindex=16; Zindex=64; % Rindex=28; Zindex=floor(size(M.zgrid,1)/2)+1; -Rp=M.R(:,1,false); -Zp=M.Z(:,1,false); -Indices=Rp>=M.rgrid(Rindex) & Rp=M.zgrid(Zindex) & Zp=M.rgrid(Rindex) & Rend=M.zgrid(Zindex) & Zend=M.rgrid(Rindex)-deltarm & Rp=M.zgrid(Zindex)-deltaz & Zp=M.rgrid(Rindex)-deltarm & Rend=M.zgrid(Zindex)-deltaz & Zend1) -h1=histogram(Vr,'Binwidth',binwidth,'Normalization','count','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9)); +h1=histogram(Vr,'Binwidth',binwidth,'Normalization','count','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(timeinit)*1e9)); %h1=histfit(ax1,Vr); set(h1,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9)); end hold on h1=histogram(Vrend,'Binwidth',binwidth,'Normalization','count','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9)); %h1=histfit(ax1,Vrend); -set(h1,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9)); +set(h1,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(timesteppart)*1e9)); ylabel('counts') -xlabel('\beta_r') - +xlabel('v_r [m/s]') grid on +ylimits=ylim(); +Veb=(-M.Er(Rindex,Zindex,end)*M.Bz(Zindex,Rindex)+M.Ez(Rindex,Zindex,end)*M.Br(Zindex,Rindex))/M.B(Zindex,Rindex)^2; +plot(-Veb*[1, 1],ylimits,'--k','DisplayName','-V_{ExB}','linewidth',1.5) +plot(Veb*[1, 1],ylimits,'--k','DisplayName','V_{ExB}','linewidth',1.5) ax2=subplot(1,3,2); if(length(Vthet)>1) -h1=histogram(Vthet(:,1),'Binwidth',binwidth,'Normalization','count','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9)); +h1=histogram(Vthet(:,1),'Binwidth',binwidth,'Normalization','count','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(timeinit)*1e9)); set(h1,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9)); end %h1=histfit(ax2,Vthet(:,1)); hold on h1=histogram(Vthetend,'Binwidth',binwidth,'Normalization','count','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9)); %h1=histfit(ax2,Vthetend); -set(h1,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9)); +set(h1,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(timesteppart)*1e9)); + ylabel('counts') -xlabel('\beta_\theta') +xlabel('v_\theta [m/s]') legend(ax2,'location','northoutside','orientation','vertical') grid on +ylimits=ylim(); +Veb=(-M.Er(Rindex,Zindex,end)*M.Bz(Zindex,Rindex)+M.Ez(Rindex,Zindex,end)*M.Br(Zindex,Rindex))/M.B(Zindex,Rindex)^2; +plot(Veb*[1, 1],ylimits,'--k','DisplayName','V_{ExB}','linewidth',1.5) +plot(-Veb*[1, 1],ylimits,'--k','DisplayName','V_{ExB}','linewidth',1.5) ax3=subplot(1,3,3); -%h1=histogram(Vz(:,1),'Binwidth',binwidth,'Normalization','probability','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9)); +%h1=histogram(Vz(:,1),'Binwidth',binwidth,'Normalization','probability','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(timeinit)*1e9)); if(length(Vz)>1) -h1=histfit(ax3,Vz(:,1)); -set(h1,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9)); +h2=histfit(ax3,Vz(:,1)); +set(h2,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9)); +h2(1).FaceAlpha=0.6; +h2(1).FaceColor=[0.00,0.45,0.74]; +h2(2).Color='blue'; end hold on %h1=histogram(Vzend,'Binwidth',binwidth,'Normalization','probability','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9)); -h1=histfit(ax3,Vzend); -set(h1,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9)); +h2=histfit(ax3,Vzend); +set(h2,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(timesteppart)*1e9)); +h2(1).FaceAlpha=0.6; +h2(1).FaceColor=[0.85,0.33,0.10]; +h2(2).Color='red'; ylabel('counts') -xlabel('\beta_z') +xlabel('v_z [m/s]') +ax3.Children=ax3.Children([1,3,2,4]); grid on f=gcf; -sgtitle(sprintf('R=%1.2e[m] Z=%1.2e[m] dt=%1.2e[ns]',M.rgrid(Rindex+1),M.zgrid(Zindex+1),M.dt*1e9)) +sgtitle(sprintf('R=%1.2e[m] Z=%1.2e[m] dt=%1.2e[ns]',M.rgrid(Rindex),M.zgrid(Zindex),M.dt*1e9)) f.PaperOrientation='landscape'; f.PaperUnits='centimeters'; papsize=[16 10]; f.PaperSize=papsize; -[~, name, ~] = fileparts(M.file); -print(f,sprintf('%sParts_V_RZ',name),'-dpdf','-fillpage') -savefig(f,sprintf('%sParts_V_RZ',name)) - - -%% Shows the particles phase space at position (rindex, zindex) -% and times tinit and timesteppart -timeinit=1; -timestepNinit=find(M.tpart(timeinit)==M.t2d); -timesteppart=length(M.tpart); -timestepNend=find(M.tpart(timesteppart)==M.t2d); - -Rp=M.R(:,timeinit,false); -Zp=M.Z(:,timeinit,false); -Indices=Rp>=M.rgrid(Rindex) & Rp=M.zgrid(Zindex) & Zp=M.rgrid(Rindex) & Rend=M.zgrid(Zindex) & Zend1) +h1=histogram(Vperp,'Binwidth',binwidth,'Normalization','count','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(timeinit)*1e9)); +%h1=histfit(ax1,Vr); +set(h1,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9)); +end hold on -p(2)=scatter(Vparend,Vperpend,'.','displayname',legendend); -xlabel('v_{par}/c') -ylabel('v_\perp/c') -legend('location','southoutside') +h1=histogram(Vperpend,'Binwidth',binwidth,'Normalization','count','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9)); +%h1=histfit(ax1,Vrend); +set(h1,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(timesteppart)*1e9)); +ylabel('counts') +grid on +ylimits=ylim(); +plot(mean(Vthetend)*[1, 1],ylimits,'--k','DisplayName','','linewidth',1.5) +if gcs + xlabel('v_{perp}^* [m/s]') + plot(mean(Vperpend)*[1, 1],ylimits,'--g','DisplayName','','linewidth',1.5) +else + xlabel('v_{perp} [m/s]') + plot(mean(Vperpend)*[1, 1],ylimits,'--g','DisplayName','','linewidth',1.5) +end +plot(Veb*[1, 1],ylimits,'--b','DisplayName','V_{ExB}','linewidth',1.5) +legend(ax1,'location','northoutside','orientation','vertical') -zleftlim=1; -vpar=linspace(1,M.vlight,1000); -dN=1e13*M.weight; -dN=0; -rp=9.85e-3; -rm=7e-3; -V=2*pi*(rp^2-rm^2); -N=M.N(:,:,timestepNend); -levels=M.rAthet(Rindex,Zindex)*[1 1]; -rAthetposend=find(levels(1)>=M.rAthet(:,zleftlim),1,'last'); -rleft=sqrt(levels(1)/(M.rAthet(rAthetposend,1)/M.rgrid(rAthetposend)^2)); -rposleft=find(M.rgrid>=rleft,1,'first'); -n2=mean(mean(N(rposleft,[1 end],:),3)); -vpar2=(dN/V/n2)^2; -if n2==0 - vpar2=0; +ax2=subplot(1,3,2); +if(length(Vpar)>1) +h1=histogram(Vpar(:,1),'Binwidth',binwidth,'Normalization','count','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(timeinit)*1e9)); +set(h1,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9)); end -b=M.rgrid(end); -a=M.rgrid(1); -vd1=((M.potout-M.potinn)*M.phinorm/M.rgrid(Rindex))/M.Bz(Zindex,Rindex)/log(b/a); -vd1=-M.Er(Rindex,Zindex,1)/M.Bz(Zindex,Rindex); -%vd1=-M.Er(Rindex,Zindex,timestepNend)/M.Bz(Zindex,Rindex); -vinit=sqrt(M.kb*10000/M.me); -%vinit=sqrt(M.qe*80/M.me) -%vinit=sqrt(120*M.qe/M.me) -%vd2=((M.potout-M.potinn)*M.phinorm/M.rgrid(25))/M.Bz(1,25)/log(b/a); -[Zmesh,Rmesh]=meshgrid(M.zgrid,M.rAthet(:,floor(M.nz/2)+1)); -Rcurv=griddata(Zmesh,M.rAthet,M.B',M.zgrid(zleftlim),M.rAthet(Rindex,Zindex),'natural')/M.B(Zindex,Rindex); -deltaphicomp=-0.5*M.me/M.qe*(vpar2-(vd1+vinit)^2*(1-M.Rcurv)) -ratioparper=vpar2/((vd1+vinit)^2*(1-M.Rcurv)) - -Ndistrib=M.N(:,:,timestepNend); -model=M.potentialwellmodel(timestepNend); -z=model.z; -r=model.r; -pot=model.pot; -rathet=model.rathet; -Zeval=[M.zgrid(zleftlim) M.zgrid(Zindex)]; -Psieval=M.rAthet(Rindex,Zindex)*[1 1]; -phis=griddata(Zmesh,M.rAthet,M.pot(:,:,end),Zeval,Psieval,'natural'); -deltaphi=-diff(phis) - -%deltaphi=-215; -R=griddata(Zmesh,M.rAthet,M.B',M.zgrid(zleftlim),M.rAthet(Rindex,Zindex),'natural')/M.B(Zindex,Rindex); -vper=sqrt(2*M.qe/M.me*deltaphi/(R-1)+vpar.^2/(R-1))/M.vlight; -vpar=vpar/M.vlight; -vpar=vpar(real(vper)~=0); -vper=vper(real(vper)~=0); -xlimits=xlim; -ylimits=ylim; -p(3)=plot(vpar,vper,'b-','displayname','Loss parabolla simul'); hold on -plot(-vpar,vper,'b-') - -vpar=linspace(1,M.vlight,1000); -vper=sqrt(-2*M.qe/M.me*deltaphicomp/(R-1)+vpar.^2/(R-1))/M.vlight; -vpar=vpar/M.vlight; -vpar=vpar(real(vper)~=0); -vper=vper(real(vper)~=0); -p(4)=plot(vpar,vper,'k--','displayname','Loss parabolla prediction'); -hold on -plot(-vpar,vper,'k--') -xlim(xlimits) -ylim(ylimits) -legend(p) -axis equal - -% subplot(2,2,2) -% scatter(Zp,Vpar,'.','displayname','Init') -% hold on -% scatter(Zend,Vparend,'.','displayname','End') -% xlabel('z [m]') -% ylabel('\beta_{par}') -% -% subplot(2,2,3) -% scatter(Vperp,Rp,'.','displayname','Init') -% hold on -% scatter(Vperpend,Rend,'.','displayname','End') -% xlabel('\beta_\perp') -% ylabel('R [m]') - -ax1=subplot(1,2,2); -h=contourf(ax1,M.zgrid,M.rgrid,Ndistrib); -hold on -[r,z]=find(Ndistrib~=0); -xlim(ax1,[M.zgrid(min(z)) M.zgrid(max(z))]) -ylim(ax1,[M.rgrid(min(r)) M.rgrid(max(r))]) -xlabel(ax1,'Z [m]') -ylabel(ax1,'R [m]') -c = colorbar(ax1); -c.Label.String= 'n [m^{-3}]'; -Zx=(M.zgrid(Zindex)); -Rx=(M.rgrid(Rindex)); - -plot(ax1,Zx,Rx,'rx','Markersize',12); - -sgtitle(sprintf('r=%1.3g [m] z=%1.3g [m] \\Delta\\phi=%1.3g[kV] R=%1.3g',M.rgrid(Rindex),M.zgrid(Zindex),(M.potout-M.potinn)*M.phinorm,M.Rcurv)) - -M.savegraph(f,sprintf('%s/%s_phasespaceR%dZ%dpos',M.folder,M.name,Rindex,Zindex),[16,12]) - +h1=histogram(Vparend,'Binwidth',binwidth,'Normalization','count','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9)); +%h1=histfit(ax2,Vthetend); +set(h1,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(timesteppart)*1e9)); +ylabel('counts') +xlabel('v_{par} [m/s]') +legend(ax2,'location','northoutside','orientation','vertical') +grid on -%% -f=figure; -p(1)=scatter(Vpar,Vperp,'.','displayname','Initial'); +ax2=subplot(1,3,3); +if(length(Vpar)>1) +h1=histogram(sqrt(Vr.^2+Vthet.^2+Vz.^2),'Binwidth',binwidth,'Normalization','count','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9)); +set(h1,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(timeinit)*1e9)); +end hold on -p(2)=scatter(Vparend,Vperpend,'.','displayname','Final'); -xlabel('v_{par}/c') -ylabel('v_\perp/c') - +h1=histogram(sqrt(Vrend.^2+Vthetend.^2+Vzend.^2),'Binwidth',binwidth,'Normalization','count','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9)); +%h1=histfit(ax2,Vthetend); +set(h1,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(timesteppart)*1e9)); +ylabel('counts') +xlabel('|v|') +legend(ax2,'location','northoutside','orientation','vertical') +grid on -zleftlim=1; +% ax3=subplot(1,3,3); +% h1=histogram(Vz(:,1),'Binwidth',binwidth,'Normalization','probability','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(timeinit)*1e9)); +% if(length(Vz)>1) +% h2=histfit(ax3,sqrt(Vr.^2+Vthet.^2+Vz.^2),[],'Rayleigh'); +% set(h2,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9)); +% h2(1).FaceAlpha=0.6; +% h2(1).FaceColor=[0.00,0.45,0.74]; +% h2(2).Color='blue'; +% end +% hold on +% h1=histogram(Vzend,'Binwidth',binwidth,'Normalization','probability','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9)); +% h2=histfit(ax3,sqrt(Vrend.^2+Vthetend.^2+Vzend.^2),[],'Rayleigh'); +% set(h2,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(timesteppart)*1e9)); +% h2(1).FaceAlpha=0.6; +% h2(1).FaceColor=[0.85,0.33,0.10]; +% h2(2).Color='red'; +% ylabel('counts') +% xlabel('|v| [m/s]') +% ax3.Children=ax3.Children([1,3,2,4]); +% grid on -vpar=linspace(1,M.vlight,1000); -dN=1e13*M.weight; -dN=0; -rp=9.85e-3; -rm=7e-3; -V=2*pi*(rp^2-rm^2); -N=M.N(:,:,timestepNend); -levels=M.rAthet(Rindex,Zindex)*[1 1]; -rAthetposend=find(levels(1)>=M.rAthet(:,zleftlim),1,'last'); -rleft=sqrt(levels(1)/(M.rAthet(rAthetposend,1)/M.rgrid(rAthetposend)^2)); -rposleft=find(M.rgrid>=rleft,1,'first'); -n2=mean(mean(N(rposleft,[1 end],:),3)); -vpar2=(dN/V/n2)^2; -if n2==0 - vpar2=0; +f=gcf; +sgtitle(sprintf('R=%1.2e[m] Z=%1.2e[m] dt=%1.2e[ns]',M.rgrid(Rindex),M.zgrid(Zindex),M.dt*1e9)) +f.PaperOrientation='landscape'; +f.PaperUnits='centimeters'; +papsize=[16 10]; +f.PaperSize=papsize; +name=strcat(M.folder,'/',M.name); +if gcs + print(f,sprintf('%sParts_V_perparstarR%dZ%d',name,Rindex,Zindex),'-dpdf','-fillpage') + savefig(f,sprintf('%sParts_V_perparstarR%dZ%d',name,Rindex,Zindex)) +else + print(f,sprintf('%sParts_V_perparR%dZ%d',name,Rindex,Zindex),'-dpdf','-fillpage') + savefig(f,sprintf('%sParts_V_perparR%dZ%d',name,Rindex,Zindex)) end -b=M.rgrid(end); -a=M.rgrid(1); -vd1=((M.potout-M.potinn)*M.phinorm/M.rgrid(Rindex))/M.Bz(Zindex,Rindex)/log(b/a); -vd1=-M.Er(Rindex,Zindex,1)/M.Bz(Zindex,Rindex); -%vd1=-M.Er(Rindex,Zindex,timestepNend)/M.Bz(Zindex,Rindex); -vinit=sqrt(M.kb*10000/M.me); -%vinit=sqrt(M.qe*80/M.me) -%vinit=sqrt(120*M.qe/M.me) -%vd2=((M.potout-M.potinn)*M.phinorm/M.rgrid(25))/M.Bz(1,25)/log(b/a); -[Zmesh,Rmesh]=meshgrid(M.zgrid,M.rAthet(:,floor(M.nz/2)+1)); -Rcurv=griddata(Zmesh,M.rAthet,M.B',M.zgrid(zleftlim),M.rAthet(Rindex,Zindex),'natural')/M.B(Zindex,Rindex); -deltaphicomp=-0.5*M.me/M.qe*(vpar2-(vd1+vinit)^2*(1-M.Rcurv)) -ratioparper=vpar2/((vd1+vinit)^2*(1-M.Rcurv)) - -Ndistrib=M.N(:,:,timestepNend); -model=M.potentialwellmodel(timestepNend); -z=model.z; -r=model.r; -pot=model.pot; -rathet=model.rathet; -Zeval=[M.zgrid(zleftlim) M.zgrid(Zindex)]; -Psieval=M.rAthet(Rindex,Zindex)*[1 1]; -phis=griddata(Zmesh,M.rAthet,M.pot(:,:,end),Zeval,Psieval,'natural'); -deltaphi=-diff(phis) - -%deltaphi=-215; -R=griddata(Zmesh,M.rAthet,M.B',M.zgrid(zleftlim),M.rAthet(Rindex,Zindex),'natural')/M.B(Zindex,Rindex); -vper=sqrt(2*M.qe/M.me*deltaphi/(R-1)+vpar.^2/(R-1))/M.vlight; -vpar=vpar/M.vlight; -vpar=vpar(real(vper)~=0); -vper=vper(real(vper)~=0); -xlimits=xlim; -ylimits=ylim; -p(3)=plot(vpar,vper,'b-','displayname','Simulation'); -hold on -plot(-vpar,vper,'b-') - -vpar=linspace(1,M.vlight,1000); -vper=sqrt(-2*M.qe/M.me*deltaphicomp/(R-1)+vpar.^2/(R-1))/M.vlight; -vpar=vpar/M.vlight; -vpar=vpar(real(vper)~=0); -vper=vper(real(vper)~=0); -p(4)=plot(vpar,vper,'k--','displayname','Prediction'); -hold on -plot(-vpar,vper,'k--') -xlim(xlimits) -ylim(ylimits) -legend(p) -axis equal -legend('location','northeast') - -title(sprintf('r=%1.3g [m] z=%1.3g [m] \\Delta\\phi=%1.3g[kV] R=%1.3g',M.rgrid(Rindex),M.zgrid(Zindex),(M.potout-M.potinn)*M.phinorm,M.Rcurv)) -M.savegraph(f,sprintf('%s/%s_phasespaceR%dZ%d',M.folder,M.name,Rindex,Zindex),[12,8]) \ No newline at end of file +%% Show the phase space at begining and end of the simulation at position (Rindex,Zindex) as a probability density +f=figure(); +ax1=subplot(1,2,1); +ax2=subplot(1,2,2); +%[p,maxnb,c]=M.displayPhaseSpace(1:2,Rindex,Zindex,'',sprintf('t=%1.3g [s]',M.tpart(1)),ax1); +% ax1=subplot(1,2,1); +% hold(ax1,'off') +c=cell(2,1); +c{1}=linspace(-3e6,3e6,21); +c{2}=linspace(-0e6,1.3e7,51); +[p,maxnb,c]=M.displayPhaseSpace('parper',1:2,Rindex,Zindex,'',sprintf('t=%1.3g [s]',M.tpart(1)),ax1,-1,c); +M.displayPhaseSpace('parper',length(M.tpart)+(-1:0),Rindex,Zindex,'',sprintf('t=%1.3g [s]',M.tpart(end)),ax2,maxnb,c); +xlimits=xlim(ax1); +xlimits=[min([xlimits,xlim(ax2)]) max([xlimits,xlim(ax2)])]; +ylimits=ylim(ax1); +ylimits=[min([ylimits,ylim(ax2)]) max([ylimits,ylim(ax2)])]; +xlim([ax1,ax2],xlimits) +ylim([ax1,ax2],ylimits) +climitsmax=max([caxis(ax1),caxis(ax2)]); +caxis(ax1,[0 climitsmax]); +caxis(ax2,[0 climitsmax]); +xtickformat(ax1,'%.3g') +ytickformat(ax1,'%.3g') +ax1.YAxis.Exponent = 0; +ax1.XAxis.Exponent = 0; + +xtickformat(ax2,'%.3g') +ytickformat(ax2,'%.3g') +ax2.YAxis.Exponent = 0; +ax2.XAxis.Exponent = 0; + +sgtitle(sprintf('r=%1.2f [mm] z=%1.2f [mm] \\Delta\\phi=%1.1f[kV] R=%1.1f',M.rgrid(Rindex)*1e3,M.zgrid(Zindex)*1e3,(M.potout-M.potinn)*M.phinorm/1e3,M.Rcurv)) +M.savegraph(f,sprintf('%s/%s_phasespaceR%dZ%dbegend',M.folder,M.name,Rindex,Zindex),[15,10]); + + +%% Show the phase space at begining and end of the simulation at position (Rindex,Zindex) as a probability density +f=figure(); +ax1=subplot(1,2,1); +ax2=subplot(1,2,2); +%[p,maxnb,c]=M.displayPhaseSpace(1:2,Rindex,Zindex,'',sprintf('t=%1.3g [s]',M.tpart(1)),ax1); +% ax1=subplot(1,2,1); +% hold(ax1,'off') +c=cell(2,1); +c{1}=linspace(-1e7,1e7,21); +c{2}=linspace(-1.5e7,1.5e7,51); +[p,maxnb,c]=M.displayPhaseSpace('rthet',1:2,Rindex,Zindex,'',sprintf('t=%1.3g [s]',M.tpart(1)),ax1,-1,c); +M.displayPhaseSpace('rthet',length(M.tpart)+(-1:0),Rindex,Zindex,'',sprintf('t=%1.3g [s]',M.tpart(end)),ax2,maxnb,c); +xlimits=xlim(ax1); +xlimits=[min([xlimits,xlim(ax2)]) max([xlimits,xlim(ax2)])]; +ylimits=ylim(ax1); +ylimits=[min([ylimits,ylim(ax2)]) max([ylimits,ylim(ax2)])]; +xlim([ax1,ax2],xlimits) +ylim([ax1,ax2],ylimits) +climitsmax=max([caxis(ax1),caxis(ax2)]); +caxis(ax1,[0 climitsmax]); +caxis(ax2,[0 climitsmax]); +xtickformat(ax1,'%.3g') +ytickformat(ax1,'%.3g') +ax1.YAxis.Exponent = 0; +ax1.XAxis.Exponent = 0; + +xtickformat(ax2,'%.3g') +ytickformat(ax2,'%.3g') +ax2.YAxis.Exponent = 0; +ax2.XAxis.Exponent = 0; + +sgtitle(sprintf('r=%1.2f [mm] z=%1.2f [mm] \\Delta\\phi=%1.1f[kV] R=%1.1f',M.rgrid(Rindex)*1e3,M.zgrid(Zindex)*1e3,(M.potout-M.potinn)*M.phinorm/1e3,M.Rcurv)) +M.savegraph(f,sprintf('%s/%s_phasespaceR%dZ%dbegendrthet',M.folder,M.name,Rindex,Zindex),[15,10]); diff --git a/matlab/espic2dhdf5.m b/matlab/espic2dhdf5.m index 409f37a..c0da0cf 100644 --- a/matlab/espic2dhdf5.m +++ b/matlab/espic2dhdf5.m @@ -1,858 +1,1031 @@ classdef espic2dhdf5 %espic2dhdf5 General class used to treat hdf5 result files of espic2d code % A result file is loaded with a call to M=espic2dhdf5(filename) where filename is the relative or absolute file path - % after loading, several quantities and composite diagnostics such as moments of the distribution function or individual particles + % after loading, several quantities and composite diagnostics such as moments of the distribution function or individual particles % quantities can be accessed. properties filename name folder fullpath timestamp info t0d t1d t2d tpart it0 it1 it2 %% Physical constants vlight=299792458; qe=1.60217662E-19; me=9.109383E-31; eps_0=8.85418781762E-12; kb=1.38064852E-23; %% Run parameters dt % simulation time step nrun % number of time steps simulated nlres nlsave nlclassical % Was the equation of motion solved in the classical framework nlPhis % Was the self-consistent electric field computed nz % number of intervals in the z direction for the grid nnr % number of intervals in the r direction for the grid for each of the 3 mesh regions lz % physical axial dimension of the simulation space nplasma % Number of initial macro particles potinn % Normalized electric potential at the coaxial insert potout % Normalized electric potential at the cylinder surface B0 % Normalization for the magnetic field Rcurv % Magnetic mirror ratio width % Magnetic mirror length n0 % Initial particle density in case of old particle loading temp % Initial particle temperature in case of old particle loading femorder % finite element method order in z and r direction ngauss % Order of the Gauss integration method for the FEM plasmadim % initial dimensions of the plasma for the old particle loading system radii % Radial limits of the three mesh regions coarse,fine,coarse H0 % Initial particle Energy for Davidsons distribution function P0 % Initial particle Angular momentum for Davidsons distribution function normalized % Are the parts quantities normalized in the h5 file - nbspecies % Number of species simulated + nbspecies % Number of species simulated %% Frequencies omepe % Reference plasma frequency used for normalization omece % Reference cyclotronic frequency for normalization %% Normalizations tnorm % Time normalization rnorm % Dimension normalization bnorm % Magnetic field normalization enorm % Electric field normalization phinorm % Electric potential normalization vnorm % Velocity normalization %% Grid data rgrid % Radial grid position points zgrid % Axial grid position points dz % Axial grid step dr % Radial grid step for the three mesh regions CellVol % Volume of the cell used for density calculation %% Magnetic field Br % Radial magnetic field Bz % Axial magnetic field Athet % Azimuthal component of the Magnetic potential vector rAthet % r*Athet used for the representation of magnetic field lines B % Magnetic field amplitude %% Energies epot % Time evolution of the particles potential energy ekin % Time evolution of the particles kinetic energy etot % Time evolution of the particles total energy etot0 % Time evolution of the reference particle total energy eerr % Time evolution of the error on the energy conservation %% 2D time data evaluated on grid points N % main specie Density fluidUR % main specie radial fluid velocity fluidUZ % main specie axial fluid velocity fluidUTHET % main specie azimuthal fluid velocity pot % Electric potential Er % Radial electric field Ez % Axial electric field Presstens % Pressure tensor %% Splines knotsr % Spline radial knots knotsz % Spline axial knots %% Particle parameters weight % Macro particle numerical weight of the main specie qsim % Macro particle charge msim % Macro particle mass nbparts % Time evolution of the number of simulated particles partepot % Electric potential at the particles positions R % Particles radial position Z % Particles axial position Rindex % Particles radial grid index Zindex % Particles axial grid index partindex % Particles unique id for tracing trajectories VR % Particles radial velocity VZ % Particles axial velocity VTHET % Particles azimuthal velocity THET % Particles azimuthal position species % Array containing the other simulated species %% Celldiag celldiag % Array containing the cell diagnostic data nbcelldiag % Total number of cell diagnostics end methods function file=file(obj) - % returns the h5 file name + % returns the h5 file name file=obj.filename; end function obj = espic2dhdf5(filename,readparts) % Reads the new result file filename and read the parts data if readparts==true % Try catch are there for compatibility with older simulation files filedata=dir(filename); if (isempty(filedata)) error("File: ""%s"" doesn't exist",filename) end obj.folder=filedata.folder; obj.filename=filename; [~, obj.name, ~] = fileparts(obj.filename); obj.fullpath=[obj.folder,'/',obj.filename]; obj.timestamp=filedata.date; if nargin==1 readparts=true; end %obj.info=h5info(filename); %% Read the run parameters obj.dt = h5readatt(obj.fullpath,'/data/input.00/','dt'); obj.nrun = h5readatt(obj.fullpath,'/data/input.00/','nrun'); obj.nlres = strcmp(h5readatt(obj.fullpath,'/data/input.00/','nlres'),'y'); obj.nlsave = strcmp(h5readatt(obj.fullpath,'/data/input.00/','nlsave'),'y'); obj.nlclassical =strcmp(h5readatt(obj.fullpath,'/data/input.00/','nlclassical'),'y'); obj.nlPhis =strcmp(h5readatt(obj.fullpath,'/data/input.00/','nlPhis'),'y'); obj.nz = h5readatt(obj.fullpath,'/data/input.00/','nz'); obj.nnr = h5read(obj.fullpath,'/data/input.00/nnr'); obj.lz = h5read(obj.fullpath,'/data/input.00/lz'); obj.qsim = h5readatt(obj.fullpath,'/data/input.00/','qsim'); obj.msim = h5readatt(obj.fullpath,'/data/input.00/','msim'); try obj.weight=h5readatt(obj.fullpath,'/data/part/','weight'); catch obj.weight=obj.msim/obj.me; end obj.nplasma = h5readatt(obj.fullpath,'/data/input.00/','nplasma'); obj.potinn = h5readatt(obj.fullpath,'/data/input.00/','potinn'); obj.potout = h5readatt(obj.fullpath,'/data/input.00/','potout'); obj.B0 = h5readatt(obj.fullpath,'/data/input.00/','B0'); obj.Rcurv = h5readatt(obj.fullpath,'/data/input.00/','Rcurv'); obj.width = h5readatt(obj.fullpath,'/data/input.00/','width'); obj.n0 = h5readatt(obj.fullpath,'/data/input.00/','n0'); obj.temp = h5readatt(obj.fullpath,'/data/input.00/','temp'); try obj.it0 = h5readatt(obj.fullpath,'/data/input.00/','it0d'); obj.it1 = h5readatt(obj.fullpath,'/data/input.00/','it2d'); obj.it2 = h5readatt(obj.fullpath,'/data/input.00/','itparts'); catch obj.it0 = h5readatt(obj.fullpath,'/data/input.00/','it0'); obj.it1 = h5readatt(obj.fullpath,'/data/input.00/','it1'); obj.it1 = h5readatt(obj.fullpath,'/data/input.00/','it2'); end try obj.nbspecies=h5readatt(obj.fullpath,'/data/input.00/','nbspecies'); obj.normalized=strcmp(h5readatt(obj.fullpath,'/data/input.00/','rawparts'),'y'); catch obj.nbspecies=1; obj.normalized=false; end try obj.nbcelldiag=h5readatt(obj.fullpath,'/data/celldiag/','nbcelldiag'); catch obj.nbcelldiag=0; end obj.omepe=sqrt(abs(obj.n0)*obj.qe^2/(obj.me*obj.eps_0)); obj.omece=obj.qe*obj.B0/obj.me; obj.nbparts= h5read(obj.fullpath, '/data/var0d/nbparts'); try obj.H0 = h5read(obj.fullpath,'/data/input.00/H0'); obj.P0 = h5read(obj.fullpath,'/data/input.00/P0'); catch obj.H0=3.2e-14; obj.P0=8.66e-25; end % Normalizations obj.tnorm=1/obj.omepe; obj.rnorm=obj.vlight*obj.tnorm; obj.bnorm=obj.B0; obj.enorm=obj.vlight*obj.bnorm; obj.phinorm=obj.enorm*obj.rnorm; obj.vnorm=obj.vlight; % Grid data obj.rgrid= h5read(obj.fullpath, '/data/var1d/rgrid')*obj.rnorm; obj.zgrid= h5read(obj.fullpath, '/data/var1d/zgrid')*obj.rnorm; obj.dz=(obj.zgrid(end)-obj.zgrid(1))/double(obj.nz); rid=1; for i=1:length(obj.nnr) obj.dr(i)=(obj.rgrid(sum(obj.nnr(1:i))+1)-obj.rgrid(rid))/double(obj.nnr(i)); rid=rid+obj.nnr(i); end - + Br = h5read(obj.fullpath,'/data/fields/Br')*obj.bnorm; obj.Br= reshape(Br,length(obj.zgrid),length(obj.rgrid)); Bz = h5read(obj.fullpath,'/data/fields/Bz')*obj.bnorm; obj.Bz= reshape(Bz,length(obj.zgrid),length(obj.rgrid)); try Atheta = h5read(obj.fullpath,'/data/fields/Athet')*obj.bnorm; obj.Athet= reshape(Atheta,length(obj.zgrid),length(obj.rgrid)); [rmeshgrid,~]=meshgrid(obj.rgrid,obj.zgrid); obj.rAthet=(rmeshgrid.*obj.Athet)'; catch end obj.B=sqrt(obj.Bz.^2+obj.Br.^2); clear Br Bz try obj.t0d=h5read(obj.fullpath,'/data/var0d/time'); catch obj.t0d=obj.dt.*double(0:length(obj.epot)-1); end obj.femorder = h5read(obj.fullpath,'/data/input.00/femorder'); obj.ngauss = h5read(obj.fullpath,'/data/input.00/ngauss'); obj.plasmadim = h5read(obj.fullpath,'/data/input.00/plasmadim'); obj.radii = h5read(obj.fullpath,'/data/input.00/radii'); obj.epot = h5read(obj.fullpath,'/data/var0d/epot'); obj.ekin = h5read(obj.fullpath,'/data/var0d/ekin'); obj.etot = h5read(obj.fullpath,'/data/var0d/etot'); try obj.etot0 = h5read(obj.fullpath,'/data/var0d/etot0'); obj.eerr = obj.etot-obj.etot0; catch obj.eerr = obj.etot-obj.etot(2); end if(obj.normalized) obj.pot=gridquantity(obj.fullpath,'/data/fields/pot',sum(obj.nnr)+1, obj.nz+1,1); obj.Er=gridquantity(obj.fullpath,'/data/fields/Er',sum(obj.nnr)+1, obj.nz+1,1); obj.Ez=gridquantity(obj.fullpath,'/data/fields/Ez',sum(obj.nnr)+1, obj.nz+1,1); else obj.pot=gridquantity(obj.fullpath,'/data/fields/pot',sum(obj.nnr)+1, obj.nz+1,obj.phinorm); obj.Er=gridquantity(obj.fullpath,'/data/fields/Er',sum(obj.nnr)+1, obj.nz+1,obj.enorm); obj.Ez=gridquantity(obj.fullpath,'/data/fields/Ez',sum(obj.nnr)+1, obj.nz+1,obj.enorm); end try obj.t2d = h5read(obj.fullpath,'/data/fields/time'); catch info=h5info(obj.fullpath,'/data/fields/partdensity'); obj.t2d=obj.dt*(0:info.objspace.Size(2)-1)*double(obj.it1); end try info=h5info(obj.fullpath,'/data/fields/moments'); obj.femorder = h5read(obj.fullpath,'/data/input.00/femorder'); kr=obj.femorder(2)+1; obj.knotsr=augknt(obj.rgrid,kr); kz=obj.femorder(1)+1; obj.knotsz=augknt(obj.zgrid,kz); try obj.CellVol= reshape(h5read(obj.fullpath,'/data/fields/volume'),length(obj.knotsz)-kz,length(obj.knotsr)-kr); obj.CellVol=permute(obj.CellVol,[2,1,3])*obj.rnorm^3; catch zvol=fnder(spmak(obj.knotsz,ones(1,length(obj.knotsz)-kz)), -1 ); rvol=fnder(spmak(obj.knotsr,2*pi*[obj.rgrid' 2*obj.rgrid(end)-obj.rgrid(end-1)]), -1 ); ZVol=diff(fnval(zvol,obj.knotsz)); RVol=diff(fnval(rvol,obj.knotsr)); obj.CellVol=RVol(3:end-1)*ZVol(3:end-1)'; obj.CellVol=padarray(obj.CellVol,[1,1],'replicate','post'); end if(obj.normalized) obj.N=splinedensity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, 1, 1); else obj.N=splinedensity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, abs(obj.qsim/obj.qe), 1); end obj.fluidUR=splinevelocity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.vnorm, 2); obj.fluidUTHET=splinevelocity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.vnorm, 3); obj.fluidUZ=splinevelocity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.vnorm, 4); if(obj.normalized) obj.Presstens=splinepressure(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, obj.vnorm^2*obj.me); else obj.Presstens=splinepressure(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, obj.vnorm^2*obj.msim); end catch obj.CellVol=(obj.zgrid(2:end)-obj.zgrid(1:end-1))*((obj.rgrid(2:end).^2-obj.rgrid(1:end-1).^2)*pi)'; obj.CellVol=obj.CellVol'; obj.N=griddensity(obj.fullpath, '/data/fields/partdensity', sum(obj.nnr)+1, obj.nz+1, obj.CellVol, abs(obj.qsim/obj.qe), true); obj.fluidUR=gridquantity(obj.fullpath, '/data/fields/fluidur', sum(obj.nnr)+1, obj.nz+1, obj.vnorm, true); obj.fluidUTHET=gridquantity(obj.fullpath, '/data/fields/fluiduthet', sum(obj.nnr)+1, obj.nz+1, obj.vnorm, true); obj.fluidUZ=gridquantity(obj.fullpath, '/data/fields/fluiduz', sum(obj.nnr)+1, obj.nz+1, obj.vnorm, true); end if(readparts) if(obj.normalized) obj.R = h5partsquantity(obj.fullpath,'/data/part','R'); obj.Z = h5partsquantity(obj.fullpath,'/data/part','Z'); else obj.R = h5partsquantity(obj.fullpath,'/data/part','R',obj.rnorm); obj.Z = h5partsquantity(obj.fullpath,'/data/part','Z',obj.rnorm); end try obj.THET = h5partsquantity(obj.fullpath,'/data/part','THET'); catch clear obj.THET end try obj.Rindex=h5partsquantity(obj.fullpath,'/data/part','Rindex'); obj.Zindex=h5partsquantity(obj.fullpath,'/data/part','Zindex'); catch clear obj.Rindex obj.Zindex end obj.VR = h5partsquantity(obj.fullpath,'/data/part','UR',obj.vnorm); obj.VZ = h5partsquantity(obj.fullpath,'/data/part','UZ',obj.vnorm); obj.VTHET= h5partsquantity(obj.fullpath,'/data/part','UTHET',obj.vnorm); if(obj.normalized) obj.partepot = h5partsquantity(obj.fullpath,'/data/part','pot',sign(obj.qsim)*obj.qe); else obj.partepot = h5partsquantity(obj.fullpath,'/data/part','pot',sign(obj.qsim)*obj.qe*obj.phinorm); end try obj.partindex = h5partsquantity(obj.fullpath,'/data/part/','partindex'); catch end end try obj.tpart = h5read(obj.fullpath,'/data/part/time'); catch obj.tpart=obj.dt*(0:size(obj.R,2)-1)*double(obj.it2); end if(obj.nbspecies >1) obj.species=h5parts.empty; for i=2:obj.nbspecies nbparts=h5read(obj.fullpath,sprintf('%s/Nparts',sprintf('/data/part/%2d',i))); if (nbparts(1)>0) obj.species(i-1)=h5parts(obj.fullpath,sprintf('/data/part/%2d',i),obj); end end end if(obj.nbcelldiag > 0) obj.celldiag=h5parts.empty; for i=1:obj.nbcelldiag nbparts=h5read(obj.fullpath,sprintf('%s/Nparts',sprintf('/data/celldiag/%02d',i))); if (nbparts(1)>0) obj.celldiag(i)=h5parts(obj.fullpath,sprintf('/data/celldiag/%02d',i),obj); end end end end function Atheta=Atheta(obj,R,Z) halflz=(obj.zgrid(end)+obj.zgrid(1))/2; Atheta=0.5*obj.B0*(R-obj.width/pi*(obj.Rcurv-1)/(obj.Rcurv+1)... .*besseli(1,2*pi*R/obj.width).*cos(2*pi*(Z-halflz)/obj.width)); end function quantity=H(obj,indices) %% computes the total energy for the main specie particle indices{1} at time indices{2} if indices{1}==':' p=1:obj.VR.nparts; else p=indices{1}; end if indices{2}==':' t=1:length(obj.tpart); else t=indices{2}; end if size(indices,1)>2 track=indices{3}; else track=false; end quantity=0.5*obj.me*(obj.VR(p,t,track).^2+obj.VTHET(p,t,track).^2+obj.VZ(p,t,track).^2)+obj.partepot(p,t,track); end function quantity=P(obj,indices) %% computes the canonical angular momentum for the main specie particle indices{1} at time indices{2} if indices{1}==':' p=1:obj.R.nparts; else p=indices{1}; end if indices{2}==':' t=1:length(obj.tpart); else t=indices{2}; end - if size(indices,1)>2 + if size(indices,2)>2 track=indices{3}; else track=false; end quantity=obj.R(p,t,track).*(obj.VTHET(p,t,track)*obj.me+sign(obj.qsim)*obj.qe*obj.Atheta(obj.R(p,t,track),obj.Z(p,t,track))); end function quantity=Vpar(obj,varargin) %Vpar Computes the parallel velocity for the main specie particle indices{1} at time indices{2} if(~iscell(varargin)) indices=mat2cell(varargin); else indices=varargin; end if indices{1}==':' p=1:obj.R.nparts; else p=indices{1}; end if indices{2}==':' t=1:length(obj.tpart); else t=indices{2}; end - if size(indices,1)>2 + if size(indices,2)>2 track=indices{3}; else track=false; end - Zind=obj.Zindex(p,t,track)+1; - Rind=obj.Rindex(p,t,track)+1; - posind=sub2ind(size(obj.B),Zind,Rind); +% Zind=obj.Zindex(p,t,track)+1; +% Rind=obj.Rindex(p,t,track)+1; + rpos(1,1,:)=obj.rgrid; + zpos(1,1,:)=obj.zgrid; + [~,Rind]=min(abs(obj.R(p,t,track)-rpos),[],3); + [~,Zind]=min(abs(obj.Z(p,t,track)-zpos),[],3); + posind=sub2ind(size(obj.B),Zind,Rind); costhet=obj.Bz(posind)./obj.B(posind); sinthet=obj.Br(posind)./obj.B(posind); quantity=obj.VR(p,t,track).*sinthet+obj.VZ(p,t,track).*costhet; end function quantity=Vperp(obj,varargin) %Vperp Computes the perpendicular velocity in the guidind center reference frame, % for the main specie particle indices{1} at time indices{2} if(~iscell(varargin)) indices=mat2cell(varargin); else indices=varargin; end if indices{1}==':' p=1:obj.R.nparts; else p=indices{1}; end if indices{2}==':' t=1:length(obj.tpart); else t=indices{2}; end - if size(indices,1)>2 + if size(indices,2)>2 track=indices{3}; else track=false; end - Zind=obj.Zindex(p,t,track)+1; - Rind=obj.Rindex(p,t,track)+1; - posind=sub2ind(size(obj.B),Zind,Rind); + + if size(indices,2)>3 + gcs=indices{4}; + else + gcs=true; + end + %Zind=obj.Zindex(p,t,track)+1; + rpos(1,1,:)=obj.rgrid; + zpos(1,1,:)=obj.zgrid; + [~,Rind]=min(abs(obj.R(p,t,track)-rpos),[],3); + [~,Zind]=min(abs(obj.Z(p,t,track)-zpos),[],3); + posind=sub2ind(size(obj.B),Zind,Rind); costhet=obj.Bz(posind)./obj.B(posind); sinthet=obj.Br(posind)./obj.B(posind); Vdrift=zeros(size(Zind)); - for j=1:length(t) - tfield=find(obj.t2d==obj.tpart(t(j))); - timeEr=obj.Er(:,:,tfield); - timeEz=obj.Ez(:,:,tfield); - posindE=sub2ind(size(timeEr),Rind(:,j),Zind(:,j)); - Vdrift(:,j)=(timeEz(posindE).*obj.Br(posind(:,j))-timeEr(posindE).*obj.Bz(posind(:,j)))./obj.B(posind(:,j)).^2; + vr=obj.VR(p,t,track); + vz=obj.VZ(p,t,track); + vthet=obj.VTHET(p,t,track); + if gcs + for j=1:length(t) + [~,tfield]=min(abs(obj.t2d-obj.tpart(t(j)))); + timeEr=obj.Er(:,:,tfield); + timeEz=obj.Ez(:,:,tfield); + posindE=sub2ind(size(timeEr),Rind(:,j),Zind(:,j)); + VdriftE=(timeEz(posindE).*obj.Br(posind(:,j))-timeEr(posindE).*obj.Bz(posind(:,j)))./obj.B(posind(:,j)).^2; + dbr=(obj.B(posind(:,j)+1)-obj.B(posind(:,j)))./(obj.rgrid(Rind(:,j)+1)-obj.rgrid(Rind(:,j))); + dbz=(obj.B(posind(:,j)+1)-obj.B(posind(:,j)))./(obj.zgrid(Zind(:,j)+1)-obj.zgrid(Zind(:,j))); + VdriftB= sign(obj.weight)*obj.me/obj.qe*(vr(:,j).*sinthet(:,j)+vz(:,j).*costhet(:,j)).^2./ ... + obj.B(posind(:,j)).^3.*(obj.Bz(posind(:,j)).*dbr-obj.Br(posind(:,j)).*dbz); + Vdrift(:,j)= VdriftE+VdriftB; + end end - - quantity=sqrt((obj.VTHET(p,t,track)-Vdrift).^2+(obj.VR(p,t,track).*costhet-obj.VZ(p,t,track).*sinthet).^2); + quantity=sqrt((vthet-Vdrift).^2+(vr.*costhet-vz.*sinthet).^2); end function displaypsi(obj,deltat,Davidson) %% plot the initial and final radial profile at position z=0 and show the normalized enveloppe function Psi % relevant for Davidson annular distribution function f=figure('Name', sprintf('%s Psi',obj.name)); f.Name= sprintf('%s Psi',obj.name); zpos=floor(length(obj.zgrid)/2); tinit=1; tend=length(obj.t2d); deltat=cell2mat(deltat); if(obj.R.nt<2) H0=obj.H0; P0=obj.P0; else H0=mean(H(obj,{1:obj.VR.nparts,obj.VR.nt,false})); P0=mean(P(obj,{1:obj.VR.nparts,obj.VR.nt,false})); end lw=1.5; Mirrorratio=(obj.Rcurv-1)/(obj.Rcurv+1); locpot=mean(obj.pot(:,zpos,tend-deltat:tend),3); psi=1+obj.qe*locpot(:)/H0-1/(2*obj.me*H0)*(P0./obj.rgrid+obj.qe*0.5*obj.B0.*(obj.rgrid-obj.width/pi*Mirrorratio*cos(2*pi*obj.zgrid(zpos)/obj.width)*besseli(1,2*pi*obj.rgrid/obj.width))).^2; locdens=mean(obj.N(:,zpos,tend-deltat:tend),3); [maxn,In]=max(locdens);%M.N(:,zpos,tinit)); plot(obj.rgrid,obj.N(:,zpos,tinit),'bx-','DisplayName',sprintf('t=%1.2f[ns]',obj.t2d(tinit)*1e9),'linewidth',lw) hold on plot(obj.rgrid,locdens,'rx-','DisplayName',sprintf('t=[%1.2f-%1.2f] [ns] averaged',obj.t2d(tend-deltat)*1e9,obj.t2d(tend)*1e9),'linewidth',lw) plot(obj.rgrid(In-2:end),1./obj.rgrid(In-2:end)*maxn*obj.rgrid(In),'DisplayName','N=c*1/r','linewidth',lw) xlabel('r [m]') ylabel('n_e [m^{-3}]') I=find(psi>0); if (length(I)>1) I=[I(1)-2; I(1)-1; I; I(end)+1; I(end)+2]; else I=obj.nnr(1):length(psi); end rq=linspace(obj.rgrid(max(I(1),1)),obj.rgrid(I(end)),500); psiinterp=interp1(obj.rgrid(I),psi(I),rq,'pchip'); zeroindices=find(diff(psiinterp>=0),2); maxpsiinterp=max(psiinterp); plot(rq,maxn*psiinterp/abs(maxpsiinterp),'Displayname','normalized \Psi [a.u.]','linewidth',lw) ylim([0 inf]) for i=1:length(zeroindices) border=plot([rq(zeroindices(i)) rq(zeroindices(i))],[0 obj.rgrid(In)/obj.rgrid(In-2)*maxn],'k--','linewidth',lw); set(get(get(border,'Annotation'),'LegendInformation'),'IconDisplayStyle','off'); end legend xlim([0 0.02]) grid on title(sprintf('Radial density profile at z=%1.2e[m]',obj.zgrid(zpos))) obj.savegraph(f,sprintf('%sPsi',obj.name),[15 10]); end function f=displayrprofile(obj,deltat,zpos) - %% plot the initial and final radial profile at the middle of the simulation space + %% plot the initial and final radial profile at the middle of the simulation space f=figure('Name', sprintf('%s Prof',obj.name)); if nargin < 3 zpos=floor(length(obj.zgrid)/2); end tinit=1; tend=length(obj.t2d); if(iscell(deltat)) deltat=cell2mat(deltat); end lw=1.5; locdens=mean(obj.N(:,zpos,tend-deltat:tend),3); Rinv=1./obj.rgrid; Rinv(isnan(Rinv))=0; VTHET=mean(obj.fluidUTHET(:,zpos,tend-deltat:tend),3); omegare=(VTHET.*Rinv); plot(obj.rgrid,obj.N(:,zpos,tinit),'bx-','DisplayName',sprintf('t=%1.2f[ns]',obj.t2d(tinit)*1e9),'linewidth',lw) hold on plot(obj.rgrid,locdens,'rx-','DisplayName',sprintf('t=[%1.2f-%1.2f] [ns] averaged',obj.t2d(tend-deltat)*1e9,obj.t2d(tend)*1e9),'linewidth',lw) xlabel('r [m]') ylabel('n_e [m^{-3}]') legend xlim([obj.rgrid(1) obj.rgrid(sum(obj.nnr(1:2)))]) grid on yyaxis right plot(obj.rgrid,omegare,'DisplayName',sprintf('<\\omega_{re}> t=[%1.2f-%1.2f] [ns] averaged',obj.t2d(tend-deltat)*1e9,obj.t2d(tend)*1e9),'linewidth',lw) ylabel('\omega_{re} [1/s]') title(sprintf('Radial density profile at z=%1.2e[m]',obj.zgrid(zpos))) obj.savegraph(f,sprintf('%srProf',obj.name),[15 10]); - end + end function changed=ischanged(obj) %ischanged Check if the file has been changed and some data must be reloaded try filedata=dir(obj.fullpath); checkedtimestamp=filedata.date; if (max(checkedtimestamp > obj.timestamp) ) changed=true; return end changed=false; return catch changed=true; return end end function displayenergy(obj) %% Plot the time evolution of the system energy and number of simulated macro particles tmin=2; tmax=length(obj.ekin); f=figure('Name', sprintf('%s Energy',obj.name)); subplot(2,1,1) plot(obj.t0d(tmin:tmax),obj.ekin(tmin:tmax),'o-',... obj.t0d(tmin:tmax),obj.epot(tmin:tmax),'d-',... obj.t0d(tmin:tmax),obj.etot(tmin:tmax),'h-',... obj.t0d(tmin:tmax),obj.eerr(tmin:tmax),'x--') legend('ekin', 'epot', 'etot','eerr') xlabel('Time [s]') ylabel('Energies [J]') grid on xlimits=xlim(); subplot(2,1,2) try semilogy(obj.t0d(tmin:tmax),abs(obj.eerr(tmin:tmax)./obj.etot0(tmin:tmax)),'h-') catch semilogy(obj.t0d(tmin:tmax),abs(obj.eerr(tmin:tmax)/obj.etot(2)),'h-') end xlabel('t [s]') ylabel('E_{err}/E_{tot}') xlim(xlimits) grid on try yyaxis right plot(obj.t0d(tmin:tmax),abs(obj.nbparts(tmin:tmax)./obj.nbparts(1)*100),'d--') ylabel('Nparts %') %ylim([0 110]) catch end obj.savegraph(f,sprintf('%s/%sEnergy',obj.folder,obj.name)); end function SimParticles(obj) %% Plot the time evolution of the number of simulated physical particles in the main specie f=figure('Name', sprintf('%s Trapped particles',obj.name)); plot(obj.t0d,obj.nbparts*obj.weight) xlabel('t [s]') ylabel('N particles') obj.savegraph(f,sprintf('%s/%sntrapped',obj.folder,obj.name)); end function fighandle=savegraph(obj, fighandle, name, papsize) - %% Saves the given figure as a pdf and a .fig + %% Saves the given figure as a pdf and a .fig fighandle.PaperUnits='centimeters'; if (nargin < 4) papsize=[14 16]; end fighandle.PaperSize=papsize; print(fighandle,name,'-dpdf','-fillpage') savefig(fighandle,name) end function displayHP(obj,tstart) % Plot the histogramm of the total energy and canonical angular momentum at time tstart and % end time of the simulation over the full simulation space for the main specie if(iscell(tstart)) tstart=cell2mat(tstart); end if(obj.R.nt>=2) tstart=obj.R.nt; f=figure('Name', sprintf('%s HP',obj.name)); legtext=sprintf("t=%2.1f - %2.1f [ns]",obj.tpart(tstart)*1e9,obj.tpart(end)*1e9); subplot(1,2,1) partsmax=min(obj.nbparts(end),obj.R.nparts); Hloc=H(obj,{1:obj.nbparts(1),1,false}); h1=histogram(Hloc,20,'BinLimits',[min(Hloc(:)) max(Hloc(:))],'DisplayName',sprintf("t=%2.3d [ns]",obj.tpart(1)*1e9)); hold on Hloc=H(obj,{1:partsmax,obj.R.nt,false}); %,'Binwidth',h1.BinWidth h1=histogram(Hloc,20,'BinLimits',[min(Hloc(:)) max(Hloc(:))],'DisplayName',legtext); ylabel('counts') xlabel('H [J]') legend subplot(1,2,2) Ploc=P(obj,{1:obj.nbparts(1),1,false}); h2=histogram(Ploc,50,'BinLimits',[min(Ploc(:)) max(Ploc(:))],'DisplayName',sprintf("t=%2.3d [ns]",obj.tpart(1)*1e9)); hold on Ploc=P(obj,{1:partsmax,obj.R.nt,false}); histogram(Ploc,50,'BinLimits',[min(Ploc(:)) max(Ploc(:))],'DisplayName',legtext); ylabel('counts') xlabel('P [kg\cdotm^2\cdots^{-1}]') %clear P %clear H legend %xlim([0.95*h2.BinLimits(1) 1.05*h2.BinLimits(2)]) obj.savegraph(f,sprintf('%s/%sParts_HP',obj.folder,obj.name)); end - end + end function displayaveragetemp(obj) % Computes and show the particles average temperature as a function of time f=figure('Name',sprintf('%s potinn=%f part temperature',obj.name,obj.potinn*obj.phinorm)); vr2=obj.VR(:,:,false); vr2=mean(vr2.^2,1)-mean(vr2,1).^2; vz2=obj.VZ(:,:,false); vz2=mean(vz2.^2,1)-mean(vz2,1).^2; vthet2=obj.VTHET(:,:,false); vthet2=mean(vthet2.^2,1)-mean(vthet2,1).^2; plot(obj.tpart,0.5*obj.me*vr2/obj.qe,'displayname','T_r') hold on plot(obj.tpart,0.5*obj.me*vz2/obj.qe,'displayname','T_z') plot(obj.tpart,0.5*obj.me*vthet2/obj.qe,'displayname','T_{thet}') xlabel('time [s]') ylabel('T [eV]') title(sprintf('\\phi_a=%.1f kV \\phi_b=%.1f kV R=%.1f',obj.potinn*obj.phinorm/1e3,obj.potout*obj.phinorm/1e3,obj.Rcurv)) legend grid obj.savegraph(f,sprintf('%s/%s_partstemp',obj.folder,obj.name)); end function Gamma=Axialflux(obj,timestep,zpos) % Computes the axial particle flux n*Uz at timestep timestep and axial position zpos Gamma=obj.fluidUZ(:,zpos,timestep).*obj.N(:,zpos,timestep); end function I=OutCurrents(obj,timestep) % Computes the Outgoing currens at the simulation axial boundaries at timestep timestep % This is simply the surface integral of the axial flux I=zeros(2,length(timestep)); flux=obj.Axialflux(timestep,[1 obj.nz+1]); I=squeeze(trapz(obj.rgrid,flux.*obj.rgrid)*2*pi*obj.qsim/obj.weight); end function displayCurrentsevol(obj,timesteps) % Computes and display the time evolution of the outgoing currents at timesteps timesteps if nargin<2 timesteps=1:length(obj.t2d); end currents=obj.OutCurrents(timesteps); f=figure('Name',sprintf('%s Currents',obj.name)); plot(obj.t2d(timesteps),currents(1,:),'Displayname',sprintf('z=%.2f cm',obj.zgrid(1)*100)); hold on plot(obj.t2d(timesteps),-currents(2,:),'Displayname',sprintf('z=%.2f cm',obj.zgrid(end)*100)); + plot(obj.t2d(timesteps),(currents(1,:)-currents(2,:))/2,'Displayname','average'); legend('location','east') xlabel('time [s]') ylabel('I [A]') grid on title(sprintf('\\phi_b-\\phi_a=%.2g kV, R=%.1f',(obj.potout-obj.potinn)*obj.phinorm/1e3,obj.Rcurv)) obj.savegraph(f,sprintf('%s/%s_outCurrents',obj.folder,obj.name)); end function model=potentialwellmodel(obj,timestep) - % Computes the potential well at the given timestep and return the model to be able to + % Computes the potential well at the given timestep and return the model to be able to % interpolate either using grid coordinates or magnetic field line coordinates if iscell(timestep) timestep=cell2mat(timestep); end Phi=-obj.pot(:,:,timestep); contpoints=contourc(obj.zgrid,obj.rgrid,obj.rAthet,obj.rAthet(:,1)'); [x,y,zcont]=C2xyz(contpoints); k=1; zdiff=[diff(zcont),0]; potfinal=zeros(numel(cell2mat(x)),length(timestep)); pot=cell(1,size(x,2)); rmin=obj.rgrid(1); rmax=obj.rgrid(end); rathet=x; for i=1:length(timestep) locPhi=Phi(:,:,i); - for j=1:size(x,2) - %for i=1:size(pot,1) - xloc=x{j}; - yloc=y{j}; -% xloc=xloc(ylocrmin); -% yloc=yloc(ylocrmin); -% x{j}=xloc; -% y{j}=yloc; - rathet{j}=zcont(j)*ones(1,length(xloc)); - if(length(xloc)>=1) - pot{j}=interp2(obj.zgrid,obj.rgrid,locPhi,xloc,yloc,'makima'); - pot{j}=pot{j}-max(pot{j}); + for j=1:size(x,2) + %for i=1:size(pot,1) + xloc=x{j}; + yloc=y{j}; + % xloc=xloc(ylocrmin); + % yloc=yloc(ylocrmin); + % x{j}=xloc; + % y{j}=yloc; + rathet{j}=zcont(j)*ones(1,length(xloc)); + if(length(xloc)>=1) + pot{j}=interp2(obj.zgrid,obj.rgrid,locPhi,xloc,yloc,'spline'); + pot{j}=pot{j}-max(pot{j}); + end + k=k+(zdiff(j)~=0); end - k=k+(zdiff(j)~=0); - end - potfinal(:,i)=cell2mat(pot); + potfinal(:,i)=cell2mat(pot); end model.z=cell2mat(x); model.r=cell2mat(y); model.pot=potfinal; model.rathet=cell2mat(rathet); end function [pot] = PotentialWell(obj,timestep) % PotentialWell Computes the potential well at the given timestep on the grid points model=obj.potentialwellmodel(timestep); z=model.z; r=model.r; modpot=model.pot; rathet=model.rathet; [Zmesh,Rmesh]=meshgrid(obj.zgrid,obj.rgrid); pot=zeros(length(obj.zgrid),length(obj.rgrid),length(fieldstep)); for i=1:length(fieldstep) pot(:,:,i)=griddata(z,r,modpot(i,:),Zmesh,Rmesh); end end function displaypotentialwell(obj,timestep,rcoord) % Display the potential well at timestep timestep % if rcoord is true, the potential is evaluated at grid points in r,z coordinates % if false, the potential is evaluated at grid points in magnetic field line coordinates if iscell(timestep) timestep=cell2mat(timestep); end if nargin <3 rcoord=true; end f=figure('Name',sprintf('%f Potential well',obj.name)); model=obj.potentialwellmodel(timestep); z=model.z; r=model.r; pot=model.pot; rathet=model.rathet; title(sprintf('Potential well t=%1.2f [ns]',obj.t2d(timestep)*1e9)) if rcoord [Zmesh,Rmesh]=meshgrid(obj.zgrid,obj.rgrid); pot=griddata(z,r,pot,Zmesh,Rmesh); surface(obj.zgrid(1:end),obj.rgrid(1:end),pot(1:end,1:end),'edgecolor','none') xlabel('z [m]') ylabel('r [m]') xlim([obj.zgrid(1) obj.zgrid(end)]) ylim([obj.rgrid(1) obj.rgrid(end)]) else [Zmesh,Rmesh]=meshgrid(obj.zgrid,obj.rAthet(:,1)); pot=griddata(z,rathet,pot,Zmesh,Rmesh); surface(obj.zgrid(1:end),obj.rAthet(:,1),pot(1:end,1:end),'edgecolor','none') ylabel('rA_\theta [Tm^2]') xlabel('z [m]') xlim([obj.zgrid(1) obj.zgrid(end)]) ylim([obj.rAthet(1,1) obj.rAthet(sum(obj.nnr(1:2)),1)]) end - c=colorbar; - colormap('jet'); - c.Label.String= 'depth [eV]'; - obj.savegraph(f,sprintf('%s/%s_well',obj.folder,obj.name)); + c=colorbar; + colormap('jet'); + c.Label.String= 'depth [eV]'; + obj.savegraph(f,sprintf('%s/%s_well',obj.folder,obj.name)); end function Epar = Epar(obj,fieldstep) % Computes the electric field component parallel to the magnetic field line Epar=obj.Er(:,:,fieldstep).*(obj.Br./obj.B)' + (obj.Bz./obj.B)'.*obj.Ez(:,:,fieldstep); end function Eperp = Eperp(obj,fieldstep) % Computes the electric field component perpendicular to the magnetic field line Eperp=obj.Er(:,:,fieldstep).*(obj.Bz./obj.B)' - (obj.Br./obj.B)'.*obj.Ez(:,:,fieldstep); end + + function charge=totcharge(obj,fieldstep) + % Integrates the density profile over the full volume to obtain + % the total number of electrons in the volume + N=splinedensity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder,ones(size(obj.CellVol)), 1, 1); + charge=sum(sum(N(:,:,end))); + end + + function [p, maxnb, c]=displayPhaseSpace(obj,type,partsstep, Rindex, Zindex,legendtext, figtitle, f, maxnb, c, gcs) + if nargin<8 + f=figure; + f=gca; + end + if nargin<7 + figtitle=sprintf('r=%1.2f [mm] z=%1.2f [mm] \\Delta\\phi=%1.1f[kV] R=%1.1f',obj.rgrid(Rindex)*1e3,obj.zgrid(Zindex)*1e3,(obj.potout-obj.potinn)*obj.phinorm/1e3,obj.Rcurv); + end + if nargin <6 + legendtext=sprintf('t=%1.3g [s]',obj.tpart(partsstep)); + end + fieldstep=find(obj.tpart(partsstep(end))==obj.t2d,1); + + zleftlim=1; + if nargin>=10 + ctemp=c; + n=zeros(length(c{1}),length(c{2})); + else + nbins=15; + n=zeros(nbins); + end + if nargin <11 + gcs=true; + end + + + for i=1:length(partsstep) + odstep=find(obj.tpart(partsstep(i))==obj.t0d); + Rp=obj.R(1:obj.nbparts(odstep),partsstep(i),false); + Zp=obj.Z(1:obj.nbparts(odstep),partsstep(i),false); + deltar=obj.dr(2)/2; + deltarm=obj.rgrid(Rindex)-sqrt(obj.rgrid(Rindex)^2-deltar^2-2*obj.rgrid(Rindex)*deltar); + deltaz=obj.dz/2; + Indices=Rp>=obj.rgrid(Rindex)-deltarm & Rp=obj.zgrid(Zindex)-deltaz & Zp2 track=indices{3}; else track=false; end quantity=0.5*obj.m*(obj.VR(p,t,track).^2+obj.VTHET(p,t,track).^2+obj.VZ(p,t,track).^2)+obj.partepot(p,t,track); end function quantity=P(obj,indices) if indices{1}==':' p=1:obj.R.nparts; else p=indices{1}; end if indices{2}==':' t=1:length(obj.tpart); else t=indices{2}; end if size(indices,1)>2 track=indices{3}; else track=false; end quantity=obj.R(p,t,track).*(obj.VTHET(p,t,track)*obj.m+obj.q*obj.parent.Atheta(obj.R(p,t,track),obj.Z(p,t,track))); end function quantity=Vpar(obj,indices) %Vpar Computes the parallel velocity for the particle indices{1} at time indices{2} if indices{1}==':' p=1:obj.R.nparts; else p=indices{1}; end if indices{2}==':' t=1:length(obj.tpart); else t=indices{2}; end - if size(indices,1)>2 + if size(indices,2)>2 track=indices{3}; else track=false; end - Zind=obj.Zindex(p,t,track); - Rind=obj.Rindex(p,t,track); + %Zind=obj.Zindex(p,t,track); + %Rind=obj.Rindex(p,t,track); + rpos(1,1,:)=obj.parent.rgrid; + zpos(1,1,:)=obj.parent.zgrid; + [~,Rind]=min(abs(obj.R(p,t,track)-rpos),[],3); + [~,Zind]=min(abs(obj.Z(p,t,track)-zpos),[],3); posind=sub2ind(size(obj.parent.B),Zind,Rind); costhet=obj.parent.Bz(posind)./obj.parent.B(posind); sinthet=obj.parent.Br(posind)./obj.parent.B(posind); quantity=obj.VR(p,t,track).*sinthet+obj.VZ(p,t,track).*costhet; end function quantity=Vperp(obj,indices) %Vperp Computes the perpendicular velocity in the guidind center reference frame, % for the main specie particle indices{1} at time indices{2} if indices{1}==':' p=1:obj.R.nparts; else p=indices{1}; end if indices{2}==':' t=1:length(obj.tpart); else t=indices{2}; end - if size(indices,1)>2 + if size(indices,2)>2 track=indices{3}; else track=false; end - Zind=obj.Zindex(p,t,track); - Rind=obj.Rindex(p,t,track); + + if size(indices,2)>3 + gcs=indices{4}; + else + gcs=true; + end + +% Zind=obj.Zindex(p,t,track)+1; +% Rind=obj.Rindex(p,t,track)+1; + rpos(1,1,:)=obj.parent.rgrid; + zpos(1,1,:)=obj.parent.zgrid; + [~,Rind]=min(abs(obj.R(p,t,track)-rpos),[],3); + [~,Zind]=min(abs(obj.Z(p,t,track)-zpos),[],3); posind=sub2ind(size(obj.parent.B),Zind,Rind); costhet=obj.parent.Bz(posind)./obj.parent.B(posind); sinthet=obj.parent.Br(posind)./obj.parent.B(posind); Vdrift=zeros(size(Zind)); - for j=1:length(t) - tfield=find(obj.parent.t2d==obj.tpart(t(j))); - timeEr=obj.parent.Er(:,:,tfield); - timeEz=obj.parent.Ez(:,:,tfield); - posindE=sub2ind(size(timeEr),Rind(:,j),Zind(:,j)); - Vdrift(:,j)=(timeEz(posindE).*obj.Br(posind(:,j))-timeEr(posindE).*obj.Bz(posind(:,j)))./obj.B(posind(:,j)).^2; + vr=obj.VR(p,t,track); + vz=obj.VZ(p,t,track); + vthet=obj.VTHET(p,t,track); + if gcs + for j=1:length(t) + [~,tfield]=min(abs(obj.parent.t2d-obj.tpart(t(j)))); + timeEr=obj.parent.Er(:,:,tfield); + timeEz=obj.parent.Ez(:,:,tfield); + posindE=sub2ind(size(timeEr),Rind(:,j),Zind(:,j)); + VdriftE=(timeEz(posindE).*obj.parent.Br(posind(:,j))-timeEr(posindE).*obj.parent.Bz(posind(:,j)))./obj.parent.B(posind(:,j)).^2; + dbr=(obj.parent.B(posind(:,j)+1)-obj.parent.B(posind(:,j)))./(obj.parent.rgrid(Rind(:,j)+1)-obj.parent.rgrid(Rind(:,j))); + dbz=(obj.parent.B(posind(:,j)+1)-obj.parent.B(posind(:,j)))./(obj.parent.zgrid(Zind(:,j)+1)-obj.parent.zgrid(Zind(:,j))); + VdriftB= +sign(obj.weight)*obj.m/obj.q*(vr(:,j).*sinthet(:,j)+vz(:,j).*costhet(:,j)).^2./ ... + obj.parent.B(posind(:,j)).^3.*(obj.parent.Bz(posind(:,j)).*dbr-obj.parent.Br(posind(:,j)).*dbz); + Vdrift(:,j)= VdriftE+VdriftB; + end end - - quantity=sqrt((obj.VTHET(p,t,track)-Vdrift).^2+(obj.VR(p,t,track).*costhet-obj.VZ(p,t,track).*sinthet).^2); + quantity=sqrt((vthet-Vdrift).^2+(vr.*costhet-vz.*sinthet).^2); end function sref = subsref(obj,s) % obj(i) is equivalent to obj.Data(i) switch s(1).type case '.' if(strcmp(s(1).subs,'H')) sref=H(obj,s(2).subs); elseif(strcmp(s(1).subs,'P')) sref=P(obj,s(2).subs); else sref=builtin('subsref',obj,s); end case '()' sref=builtin('subsref',obj,s); case '{}' error('MYDataClass:subsref',... 'Not a supported subscripted reference') end end end end \ No newline at end of file diff --git a/matlab/timescales.m b/matlab/timescales.m index a44527f..7e06a1f 100644 --- a/matlab/timescales.m +++ b/matlab/timescales.m @@ -1,93 +1,93 @@ -Pn=1e-8;%mbar +Pn=1e-6;%mbar kb=1.38064852e-23; T=300; estimneutraldensity=Pn*100/(kb*T); nb=estimneutraldensity; qe=1.60217662E-19; me=9.109383E-31; eps_0=8.85418781762e-12; vlight=299792458; Ee=[10]; %eV gridsize1=ceil(sqrt(length(Ee))); gridsize2=ceil(length(Ee)/gridsize1); fighandle=figure; for i=1:length(Ee) subplot(gridsize2,gridsize1,i) ve=sqrt(Ee(i)*qe/me); %m/s gamma=1/(1-ve^2/vlight^2); ne=5*logspace(13,17,3000); B0=0.21; sigion=sigmabeb(15.58,54.91,2,1,Ee(i)); sigion=sigmabeb(41.72,71.13,2,1,Ee(i))+sigion; sigion=sigmabeb(21,63.18,2,1,Ee(i))+sigion; sigion=sigmabeb(17.07,44.30,4,1,Ee(i))+sigion; Coulomblog=log(sqrt(eps_0*Ee(i)/qe*ne)/qe^2*2*pi*eps_0*me*ve^2); ne_0=1; tauion=1./(nb*sigion*ve); tauloss=tauion*log(ne/ne_0); taucycl=2*pi*me*gamma/qe/B0; taupe=2*pi*sqrt(eps_0*me*gamma/qe^2./ne); tauthermal=1./(ve*ne*qe^4.*Coulomblog./(gamma^2*me^2*ve^4*2*pi*eps_0^2)); ldw=1.5; loglog(ne,tauion*ones(size(ne)),'displayname','\tau_{ion}','linewidth',ldw) hold on loglog(ne,tauloss,'displayname','t_{accum}','linewidth',ldw) loglog(ne,tauthermal,'displayname','\tau_{thermal}','linewidth',ldw) loglog(ne,taucycl*ones(size(ne)),'displayname','\tau_{ce}','linewidth',ldw) loglog(ne,taupe,'displayname','\tau_{pe}','linewidth',ldw) legend('location','eastoutside') %title(sprintf('E_e=%#.3g eV',Te(i))) xlabel('n_e [m^{-3}]') ylabel('\tau [s]') grid on end sgtitle(sprintf('E_e=%3.2g [eV], P_n=%3.2g [mbar]',Ee(1),Pn)) %% Saves the given figure as a pdf and a .fig fighandle.PaperUnits='centimeters'; papsize=[12 10]; fighandle.PaperSize=papsize; name=sprintf('timescaleE%2.0eP%2.0embar',Ee(1),Pn); print(fighandle,name,'-dpdf','-fillpage') savefig(fighandle,name) % Te=500; % ve=sqrt(Te*qe/me); %m/s % gamma=1/(1-ve^2/vlight^2); % time=logspace(-5,1,500); % dens=exp(time/tauion); % Coulomblog=log(sqrt(eps_0*Te/qe*dens)/qe^2*2*pi*eps_0*me*ve^2); % tauthermal=1./(ve*dens*qe^4.*Coulomblog./(gamma^2*me^2*ve^4*2*pi*eps_0^2)); % figure % yyaxis left % loglog(time,dens) % yyaxis right % loglog(time,tauthermal) % Sigma BEB see (https://physics.nist.gov/PhysRefData/Ionization/intro.html) % B binding energy % U average kinetic energy % N electron occupation number % Q dipole constant % T incident energy function sig=sigmabeb(B,U,N,Q,T) t=T/B; %incident energy u=U/B; % Average kinetic energy % B: a0=0.52918e-10; R=13.6057; S=4*pi*a0^2*N*(R/B)^2; n=1; sig=S/(t+(u+1)/n)*(Q*log(t)*0.5*(1-t^-2)+(2-Q)*(1-t^-1-log(t)/(t+1))); end