diff --git a/matlab/PlotParticleTrajectory.m b/matlab/PlotParticleTrajectory.m
index ffef0e1..5b632c6 100644
--- a/matlab/PlotParticleTrajectory.m
+++ b/matlab/PlotParticleTrajectory.m
@@ -1,65 +1,72 @@
function PlotParticleTrajectory(M,partnumber,t,text)
%UNTITLED3 Summary of this function goes here
% Detailed explanation goes here
if (nargin<4)
text='';
else
text=['_',text];
end
f=figure();
%time=((t(1):t(end))-1)*M.dt*double(M.it2);
time=M.tpart(t(1):t(end));
%sgtitle(sprintf('Particles %d',partnumber));
ax1=subplot(1,2,1);
hold(ax1);
ax2=subplot(1,2,2);
hold(ax2);
Msize=8;
i=1;
for part=partnumber
p1(i)=plot(ax1,M.Z(part,t),M.R(part,t),'x-.','Linewidth',1.1,...
'Displayname',sprintf('part=%d',part),'MarkerIndices',1,...
'Markersize',Msize);
plot(ax1,M.Z(part,t(end)),M.R(part,t(end)),'^','Linewidth',1.1,...
'Displayname',sprintf('part=%d',part),'MarkerIndices',1,...
'Markersize',Msize,'color',p1(i).Color)
x=M.R(part,t).*cos(M.THET(part,t));
y=M.R(part,t).*sin(M.THET(part,t));
%pp=spline(x,y);
p2=plot(ax2,x,y,'x-.','Linewidth',1.1,...
'Displayname',sprintf('part=%d',part),'MarkerIndices',1,...
'Markersize',Msize);
plot(ax2,M.R(part,t(end)).*cos(M.THET(part,t(end))),M.R(part,t(end)).*sin(M.THET(part,t(end))),'^',...
'Linewidth',1.1,'Displayname',sprintf('part=%d',part),'MarkerIndices',1,...
'Markersize',Msize,'color',p2.Color)
i=i+1;
end
xlabel(ax1,'z [m]')
ylabel(ax1,'r [m]')
sgtitle(sprintf('Position between t=[%3.2f-%3.2f] ns',time(1)*1e9,time(end)*1e9))
-%xlim(ax1,[M.zgrid(1) M.zgrid(end)])
-xlim(ax1,[-.17 .17])
-%ylim(ax1,[M.rgrid(1) M.rgrid(end)])
-ylim(ax1,[0.005 0.025])
+xlim(ax1,[M.zgrid(1) M.zgrid(end)])
+%xlim(ax1,[-.17 .17])
+ylim(ax1,[M.rgrid(1) M.rgrid(end)])
+%ylim(ax1,[0.005 0.025])
legend(ax1,p1(1:end))
grid(ax1,'on');
xlabel(ax2,'x [m]')
ylabel(ax2,'y [m]')
+t=2*pi*linspace(0,1,100);
+xin=M.rgrid(1)*cos(t);
+yin=M.rgrid(1)*sin(t);
+plot(ax2,xin,yin,'k-')
+xout=M.rgrid(end)*cos(t);
+yout=M.rgrid(end)*sin(t);
+plot(ax2,xout,yout,'k-')
axis(ax2,'equal')
%legend(ax2)
grid(ax2,'on');
% ax=subplot(1,2,2);
% xlabel(ax,'t [s]')
% ylabel(ax,'VR [m/s]')
% title(ax,'Radial velocity')
% grid on;
f.PaperOrientation='landscape';
-[~, name, ~] = fileparts(M.file);
+[~, name, ~] = fileparts(M.fullpath);
f.PaperUnits='centimeters';
f.PaperSize=[16 9];
print(f,sprintf('%sTrajectories%s',name,text),'-dpdf','-fillpage')
end
diff --git a/matlab/analyse_espicfields.m b/matlab/analyse_espicfields.m
index 8251528..fe2014f 100644
--- a/matlab/analyse_espicfields.m
+++ b/matlab/analyse_espicfields.m
@@ -1,369 +1,372 @@
% file='teststablegauss_13_traject2.h5';
% file='unigauss.h5';
% file='Dav_5e14fine_wr+.h5';
% file='stablegauss_8e19fine.h5'
% % %file='teststable_Dav1.h5';
% % %close all hidden;
% if (~exist('M','var'))
% M.file=file;
% end
% M=load_espic2d(file,M,'all');
t2dlength=size(M.t2d,1)
-fieldstart=max(1,t2dlength-800);%floor(0.95*t2dlength);
+fieldstart=max(1,t2dlength-100);%floor(0.95*t2dlength);
deltat=t2dlength-fieldstart;
%[~, rgridend]=min(abs(M.rgrid-0.045));
rgridend=sum(M.nnr(1:2));
[~, name, ~] = fileparts(M.file);
M.plotEnergy
%fieldstart=2300; fieldend=2500;
dens=mean(M.N(:,:,fieldstart:end),3);
pot=mean(M.pot(:,:,fieldstart:end),3);
brillratio=2*dens*M.me./(M.eps_0*M.Bz'.^2);
[R,Z]=meshgrid(M.rgrid,M.zgrid);
Rinv=1./R;
Rinv(:,1)=0;
VTHET=mean(M.fluidUTHET(:,:,fieldstart:end),3);
omegare=(VTHET.*Rinv');
maxdens=max(max(dens));
% f=figure();
% ax3=gca;
% surface(ax3,M.zgrid(1:end-1),M.rgrid(1:end-1),(squeeze(mean(M.Presstens(4,:,:,fieldstart:end),4)))*M.vlight)
% xlabel(ax3,'z [m]')
% ylabel(ax3,'r [m]')
% xlim(ax3,[M.zgrid(1) M.zgrid(end)])
% ylim(ax3,[M.rgrid(1) M.rgrid(50)])
% colormap(ax3,'jet')
% c = colorbar(ax3);
% c.Label.String= 'thermal v_\theta [m/s]';
% %c.Limits=[min(M.fluidUTHET(:)) max(M.fluidUTHET(:))];
% %caxis(ax3,[min(M.fluidUTHET(:)) max(M.fluidUTHET(:))])
% title(ax3,'Azimuthal velocity')
% %set(ax3,'colorscale','log')
% grid(ax3, 'on')
% view(ax3,2)
% f.PaperOrientation='landscape';
% [~, name, ~] = fileparts(M.file);
% f.PaperUnits='centimeters';
% papsize=[16 14];
% f.PaperSize=papsize;
% print(f,sprintf('%sfluid_thermtheta',name),'-dpdf','-fillpage')
%% Psi function
-displaypsi(M,deltat)
+try
+ displaypsi(M,deltat)
+catch
+end
%% Density
f=figure('Name', sprintf('%sfluid_dens',name));
ax1=gca;
h=surface(ax1,M.zgrid,M.rgrid,dens);
+%set(h,'edgecolor','none');
ylim(ax1,[M.rgrid(1) M.rgrid(rgridend)])
xlim(ax1,[M.zgrid(1) M.zgrid(end)])
xlabel(ax1,'z [m]')
ylabel(ax1,'r [m]')
title(ax1,sprintf('Density t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(fieldstart),M.t2d(end),double(maxdens)))
c = colorbar(ax1);
c.Label.String= 'n[m^{-3}]';
view(ax1,2)
-set(h,'edgecolor','none');
grid on;
+hold on;
f.PaperOrientation='landscape';
[~, name, ~] = fileparts(M.file);
f.PaperUnits='centimeters';
papsize=[16 14];
f.PaperSize=papsize;
print(f,sprintf('%sfluid_dens',name),'-dpdf','-fillpage')
savefig(f,sprintf('%sfluid_dens',name))
%% Fieldswide
f=figure('Name', sprintf('%s fields',name));
ax1=gca;
h=contourf(ax1,M.zgrid,M.rgrid,dens,'Displayname','n_e [m^{-3}]');
hold on;
%[c1,h]=contour(ax1,M.zgrid,M.rgrid,pot./1e3,5,'m--','ShowText','off','linewidth',1.5,'Displayname','Equipotentials [kV]');
%clabel(c1,h,'Color','white')
Blines=zeros(size(M.Athet));
for i=1:length(M.zgrid)
Blines(i,:)=M.Athet(i,:).*M.rgrid';
end
zindex=floor(length(M.zgrid/2));
levels=logspace( log10(min(min(Blines(:,2:end)))), log10(max(max(Blines(:,:)))),8);
contour(ax1,M.zgrid,M.rgrid,Blines',levels,'r-.','linewidth',1.5,'Displayname','Magnetic field lines');
xlim(ax1,[M.zgrid(1) M.zgrid(end)])
ylim(ax1,[M.rgrid(1) M.rgrid(end)])
legend
xlabel(ax1,'z [m]')
ylabel(ax1,'r [m]')
title(ax1,sprintf('Density t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(fieldstart),M.t2d(end),double(maxdens)))
c = colorbar(ax1);
c.Label.String= 'n[m^{-3}]';
view(ax1,2)
%set(h,'edgecolor','none');
grid on;
%rectangle('position',[M.zgrid(5) M.rgrid(3) M.zgrid(end-5)-M.zgrid(5) (M.rgrid(rgridend)-M.rgrid(1))],'Edgecolor','white','linestyle','--','linewidth',2)
f.PaperOrientation='landscape';
[~, name, ~] = fileparts(M.file);
f.PaperUnits='centimeters';
papsize=[18 10];
f.PaperSize=papsize;
print(f,sprintf('%sFieldswide',name),'-dpdf','-fillpage')
savefig(f,sprintf('%sFieldswide',name))
%% Fields
f=figure('Name', sprintf('%s fields',name));
ax1=gca;
h=contourf(ax1,M.zgrid,M.rgrid,dens,'Displayname','n_e [m^{-3}]');
hold on;
[c1,h]=contour(ax1,M.zgrid,M.rgrid,pot./1e3,'m--','ShowText','on','linewidth',1.5,'Displayname','Equipotentials [kV]');
clabel(c1,h,'Color','white')
Blines=zeros(size(M.Athet));
for i=1:length(M.zgrid)
Blines(i,:)=M.Athet(i,:).*M.rgrid';
end
zindex=floor(length(M.zgrid/2));
levels=logspace( log10(min(min(Blines(:,2:end)))), log10(max(max(Blines(:,:)))),10);
contour(ax1,M.zgrid,M.rgrid,Blines',levels,'r-.','linewidth',1.5,'Displayname','Magnetic field lines');
xlim(ax1,[M.zgrid(1) M.zgrid(end)])
ylim(ax1,[M.rgrid(1) M.rgrid(rgridend)])
legend
xlabel(ax1,'z [m]')
ylabel(ax1,'r [m]')
title(ax1,sprintf('Density t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(fieldstart),M.t2d(end),double(maxdens)))
c = colorbar(ax1);
c.Label.String= 'n[m^{-3}]';
view(ax1,2)
%set(h,'edgecolor','none');
grid on;
f.PaperOrientation='landscape';
[~, name, ~] = fileparts(M.file);
f.PaperUnits='centimeters';
papsize=[18 10];
f.PaperSize=papsize;
print(f,sprintf('%sFields',name),'-dpdf','-fillpage')
savefig(f,sprintf('%sFields',name))
%% Wells
-f=figure('Name', sprintf('%s Wells',name));
-ax1=gca;
-rpos=zeros(size(pot));
-[r,~]=meshgrid(M.rgrid,M.zgrid);
-rathet=(M.Athet.*r)';
-poteff=zeros(size(pot));
-for j=2 :size(rpos,2)-1
- for i=1:size(rpos,1)
- if(isempty(find(rathet(i,j)>rathet(:,1),1,'last')))
- rpos(i,j)=1;
- rinterp=0;
- else
- rpos(i,j)=find(rathet(i,j)>rathet(:,1),1,'last');
- rinterp=(rathet(i,j)-rathet(rpos(i,j),1))./(rathet(rpos(i,j)+1,1)-rathet(rpos(i,j),1));
- end
- bmin=M.B(1,rpos(i,j))+rinterp*(M.B(1,rpos(i,j)+1)-M.B(1,rpos(i,j)));
- potmin=pot(rpos(i,j),1)+rinterp*(pot(rpos(i,j)+1,1)-pot(rpos(i,j),1));
- poteff(i,j)=-(pot(i,j)-potmin)*bmin/(M.B(j,i)-bmin);
- end
-end
-h=surface(ax1,M.zgrid,M.rgrid,poteff,'edgecolor','none','facecolor','interp');
-Blines=rathet';
-zindex=floor(length(M.zgrid/2));
-levels=logspace( log10(min(min(Blines(:,2:end)))), log10(max(max(Blines(:,:)))),8);
-hold on
-%contour(ax1,M.zgrid,M.rgrid,Blines',levels,'r-.','linewidth',1.5,'Displayname','Magnetic field lines');
-xlim(ax1,[M.zgrid(1) M.zgrid(end)])
-ylim(ax1,[M.rgrid(1) M.rgrid(end)])
-legend
-xlabel(ax1,'z [m]')
-ylabel(ax1,'r [m]')
-title(ax1,sprintf('Ion Well t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(fieldstart),M.t2d(end),double(maxdens)))
-c = colorbar(ax1);
-c.Label.String= 'Potential [eV]';
-view(ax1,2)
-%set(h,'edgecolor','none');
-grid on;
-f.PaperOrientation='landscape';
-[~, name, ~] = fileparts(M.file);
-f.PaperUnits='centimeters';
-papsize=[16 14];
-f.PaperSize=papsize;
-print(f,sprintf('%sWell',name),'-dpdf','-fillpage')
-savefig(f,sprintf('%sWell',name))
+Wells(M,pot,[fieldstart t2dlength],maxdens)
+Wells(M,M.pot(:,:,1),1, max(max(M.N(:,:,1))))
%% Brillouin
f=figure('Name', sprintf('%s Brillouin',name));
ax1=gca;
h=surface(ax1,M.zgrid,M.rgrid,brillratio);
set(h,'edgecolor','none');
xlim(ax1,[M.zgrid(1) M.zgrid(end)])
ylim(ax1,[M.rgrid(1) M.rgrid(rgridend)])
xlabel(ax1,'z [m]')
ylabel(ax1,'r [m]')
title(ax1,'Position')
c = colorbar(ax1);
c.Label.String= 'Brillouin Ratio';
view(ax1,2)
caxis([0 1.2])
title(ax1,sprintf('Brillouin Ratio t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(fieldstart),M.t2d(end),double(maxdens)))
f.PaperOrientation='landscape';
[~, name, ~] = fileparts(M.file);
f.PaperUnits='centimeters';
papsize=[16 14];
f.PaperSize=papsize;
print(f,sprintf('%sfluid_Brillouin',name),'-dpdf','-fillpage')
savefig(f,sprintf('%sfluid_Brillouin',name))
%% Omegar
f=figure('Name', sprintf('%s Omegar',name));
ax3=gca;
surface(ax3,M.zgrid,M.rgrid,omegare,'edgecolor','None')
xlabel(ax3,'z [m]')
ylabel(ax3,'r [m]')
xlim(ax3,[M.zgrid(1) M.zgrid(end)])
ylim(ax3,[M.rgrid(1) M.rgrid(rgridend)])
colormap(ax3,'jet')
c = colorbar(ax3);
c.Label.String= '|\omega_\theta| [1/s]';
%c.Limits=[min(M.fluidUTHET(:)) max(M.fluidUTHET(:))];
%caxis(ax3,[min(M.fluidUTHET(:)) max(M.fluidUTHET(:))])
title(ax3,sprintf('Azimuthal frequency t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(fieldstart),M.t2d(end),double(maxdens)))
%set(ax3,'colorscale','log')
grid(ax3, 'on')
view(ax3,2)
f.PaperOrientation='landscape';
[~, name, ~] = fileparts(M.file);
f.PaperUnits='centimeters';
papsize=[16 14];
f.PaperSize=papsize;
print(f,sprintf('%sfluid_omegar',name),'-dpdf','-fillpage')
savefig(f,sprintf('%sfluid_omegar',name))
%% HP
if(size(M.H,2)>=2)
-f=figure('Name', sprintf('%s HP',name));
-tstart=max(floor(0.8*size(M.tpart,1)),2);
-legtext=sprintf("t=%2.1f - %2.1f [ns]",M.tpart(tstart)*1e9,M.tpart(end)*1e9);
-subplot(1,2,1)
-partsmax=min(M.nbparts(end),size(M.H,1));
-H=M.H(1:M.nbparts(1),1);
-h1=histogram(H,10,'BinLimits',[min(H(:)) max(H(:))],'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9));
-hold on
-H=mean(M.H(1:partsmax,tstart:end),2);
-%,'Binwidth',h1.BinWidth
-h1=histogram(H,20,'BinLimits',[min(H(:)) max(H(:))],'DisplayName',legtext);
-ylabel('counts')
-xlabel('H [J]')
-legend
+taveragestart=max(floor(0.8*size(M.tpart,1)),2);
+M.displayHP(taveragestart);
+end
-subplot(1,2,2)
-P=M.P(1:M.nbparts(1),1);
-h2=histogram(P,10,'BinLimits',[min(P(:)) max(P(:))],'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9));
-hold on
-P=mean(M.P(1:partsmax,tstart:end),2);
-histogram(P,50,'BinLimits',[min(P(:)) max(P(:))],'DisplayName',legtext);
-ylabel('counts')
-xlabel('P [kg\cdotm^2\cdots^{-1}]')
-%clear P
-%clear H
-legend
+%% 0D diagnostics
+BrillouinRatio=2*M.omepe^2/M.omece^2
+NpartsTrapped=mean(M.nbparts(end-10:end))*M.msim/M.me
+dblengthiter=fieldstart:length(M.t2d);
+
+Pr=squeeze(mean(M.Presstens(1,:,:,dblengthiter),4));
+Pz=squeeze(mean(M.Presstens(6,:,:,dblengthiter),4));
+enddens=dens;
+Debye_length=sqrt(abs(min(min(min(Pr(Pr>0)./enddens(Pr>0).^2)),min(min(Pz(Pz>0)./enddens(Pz>0).^2))))*M.eps_0/M.qe^2)
+
+
+%% Debye length
+f=figure('Name', sprintf('%s Debye',name));
+surface(M.zgrid,M.rgrid,min(Pr./enddens.^2,Pz./enddens.^2)*M.eps_0/M.qe^2)
+colorbar
f.PaperOrientation='landscape';
[~, name, ~] = fileparts(M.file);
-%xlim([0.95*h2.BinLimits(1) 1.05*h2.BinLimits(2)])
-print(f,sprintf('%sParts_HP',name),'-dpdf','-fillpage')
-savefig(f,sprintf('%sParts_HP',name))
-end
+f.PaperUnits='centimeters';
+papsize=[16 14];
+f.PaperSize=papsize;
+print(f,sprintf('%s_dblength',name),'-dpdf','-fillpage')
+savefig(f,sprintf('%s_dblength',name))
-% %% 0D diagnostics
-% BrillouinRatio=2*M.omepe^2/M.omece^2
-% NpartsTrapped=mean(M.nbparts(end-10:end))*M.msim/M.me
-% dblengthiter=fieldstart:length(M.t2d);
-% Pr=squeeze(mean(M.Presstens(1,:,:,dblengthiter),4));
-% Pz=squeeze(mean(M.Presstens(6,:,:,dblengthiter),4));
-% enddens=dens;
-% Debye_length=sqrt(abs(min(min(min(Pr(Pr>0)./enddens(Pr>0).^2)),min(min(Pz(Pz>0)./enddens(Pz>0).^2))))*M.eps_0/M.qe^2)
-%
-%
-% %% Debye length
-% figure
-% surface(M.zgrid,M.rgrid,min(Pr./enddens.^2,Pz./enddens.^2)*M.eps_0/M.qe^2)
-% colorbar
-%
-% %% Temperature
-% figure
-% ax1=subplot(2,1,1);
-% title('Radial temperature')
-% surface(M.zgrid,M.rgrid,Pr./enddens/M.qe,'edgecolor','none')
-% colormap('jet')
-% c = colorbar;
-% c.Label.String= 'T_er [eV]';
-% templimits1=caxis;
-%
-% ax2=subplot(2,1,2);
-% title('Axial tempreature')
-% surface(M.zgrid,M.rgrid,Pz./enddens/M.qe,'edgecolor','none')
-% colormap('jet')
-% c = colorbar;
-% c.Label.String= 'T_ez [eV]';
-% maxtemp=max([templimits1,caxis]);
-% caxis(ax1,[0 maxtemp])
-% caxis(ax2,[0 maxtemp])
+%% Temperature
+f=figure('Name', sprintf('%s Temperature',name));
+ax1=subplot(2,1,1);
+title('Radial temperature')
+surface(M.zgrid,M.rgrid,Pr./enddens/M.qe,'edgecolor','none')
+colormap('jet')
+c = colorbar;
+c.Label.String= 'T_er [eV]';
+templimits1=caxis;
+
+ax2=subplot(2,1,2);
+title('Axial tempreature')
+surface(M.zgrid,M.rgrid,Pz./enddens/M.qe,'edgecolor','none')
+colormap('jet')
+c = colorbar;
+c.Label.String= 'T_ez [eV]';
+maxtemp=max([templimits1,caxis]);
+caxis(ax1,[0 maxtemp])
+caxis(ax2,[0 maxtemp])
+f.PaperOrientation='landscape';
+[~, name, ~] = fileparts(M.file);
+f.PaperUnits='centimeters';
+papsize=[16 14];
+f.PaperSize=papsize;
+print(f,sprintf('%s_temp',name),'-dpdf','-fillpage')
+savefig(f,sprintf('%s_temp',name))
function displaypsi(M,deltat)
%% Psi function
[~, name, ~] = fileparts(M.file);
f=figure('Name', sprintf('%s Psi',name));
zpos=floor(length(M.zgrid)/2);
tinit=1;
tend=length(M.t2d);
if(size(M.H,2)<2)
H0=M.H0;
P0=M.P0;
else
H0=mean(M.H(:,end));
P0=mean(M.P(:,end));
end
lw=1.5;
Mirrorratio=(M.Rcurv-1)/(M.Rcurv+1);
locpot=mean(M.pot(:,zpos,tend-deltat:tend),3);
psi=1+M.qe*locpot(:)/H0-1/(2*M.me*H0)*(P0./M.rgrid+M.qe*0.5*M.B0.*(M.rgrid-M.width/pi*Mirrorratio*cos(2*pi*M.zgrid(zpos)/M.width)*besseli(1,2*pi*M.rgrid/M.width))).^2;
locdens=mean(M.N(:,zpos,tend-deltat:tend),3);
[maxn,In]=max(locdens);%M.N(:,zpos,tinit));
plot(M.rgrid,M.N(:,zpos,tinit),'bx-','DisplayName',sprintf('t=%1.2f[ns]',M.t2d(tinit)*1e9),'linewidth',lw)
hold on
plot(M.rgrid,locdens,'rx-','DisplayName',sprintf('t=[%1.2f-%1.2f] [ns] averaged',M.t2d(tend-deltat)*1e9,M.t2d(tend)*1e9),'linewidth',lw)
plot(M.rgrid(In-2:end),1./M.rgrid(In-2:end)*maxn*M.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=M.nnr(1):length(psi);
end
-rq=linspace(M.rgrid(I(1)),M.rgrid(I(end)),500);
+rq=linspace(M.rgrid(max(I(1),1)),M.rgrid(I(end)),500);
psiinterp=interp1(M.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])
ylimits=ylim;
for i=1:length(zeroindices)
border=plot([rq(zeroindices(i)) rq(zeroindices(i))],[0 M.rgrid(In)/M.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]',M.zgrid(zpos)))
f.PaperOrientation='landscape';
[~, name, ~] = fileparts(M.file);
f.PaperUnits='centimeters';
papsize=[15 10];
f.PaperSize=papsize;
print(f,sprintf('%sPsi',name),'-dpdf','-fillpage')
savefig(f,sprintf('%sPsi',name))
end
+
+
+function Wells(M,pot,fieldtime,maxdens)
+
+[~, name, ~] = fileparts(M.file);
+f=figure('Name', sprintf('%s Wells',name));
+ax1=gca;
+rpos=zeros(size(pot));
+[r,~]=meshgrid(M.rgrid,M.zgrid);
+rathet=(M.Athet.*r)';
+poteff=zeros(size(pot));
+for j=2 :size(rpos,2)-1
+ for i=1:size(rpos,1)
+ if(isempty(find(rathet(i,j)>rathet(:,1),1,'last')))
+ rpos(i,j)=1;
+ rinterp=0;
+ else
+ rpos(i,j)=find(rathet(i,j)>rathet(:,1),1,'last');
+ rinterp=(rathet(i,j)-rathet(rpos(i,j),1))./(rathet(rpos(i,j)+1,1)-rathet(rpos(i,j),1));
+ end
+ bmin=M.B(1,rpos(i,j))+rinterp*(M.B(1,rpos(i,j)+1)-M.B(1,rpos(i,j)));
+ potmin=pot(rpos(i,j),1)+rinterp*(pot(rpos(i,j)+1,1)-pot(rpos(i,j),1));
+ poteff(i,j)=-(pot(i,j)-potmin)*bmin/(M.B(j,i)-bmin);
+ end
+end
+h=surface(ax1,M.zgrid,M.rgrid,poteff,'edgecolor','none','facecolor','interp');
+Blines=rathet';
+zindex=floor(length(M.zgrid/2));
+levels=logspace( log10(min(min(Blines(:,2:end)))), log10(max(max(Blines(:,:)))),8);
+hold on
+%contour(ax1,M.zgrid,M.rgrid,Blines',levels,'r-.','linewidth',1.5,'Displayname','Magnetic field lines');
+xlim(ax1,[M.zgrid(1) M.zgrid(end)])
+ylim(ax1,[M.rgrid(1) M.rgrid(end)])
+legend
+xlabel(ax1,'z [m]')
+ylabel(ax1,'r [m]')
+if size(fieldtime)<2
+ title(ax1,sprintf('Ion Well t=%1.2g s n_e=%1.2gm^{-3}',M.t2d(fieldtime),double(maxdens)))
+ filename=sprintf('%sWell%d',name,fieldtime);
+else
+ title(ax1,sprintf('Ion Well t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(fieldtime(1)),M.t2d(fieldtime(2)),double(maxdens)))
+ filename=sprintf('%sWell%d_%d',name,fieldtime);
+end
+c = colorbar(ax1);
+c.Label.String= 'Potential [eV]';
+view(ax1,2)
+%set(h,'edgecolor','none');
+grid on;
+f.PaperOrientation='landscape';
+f.PaperUnits='centimeters';
+papsize=[16 14];
+f.PaperSize=papsize;
+print(f,filename,'-dpdf','-fillpage')
+savefig(f,filename)
+end
\ No newline at end of file
diff --git a/matlab/coaxial_ion_loss_criteria.m b/matlab/coaxial_ion_loss_criteria.m
new file mode 100644
index 0000000..0a8e4b6
--- /dev/null
+++ b/matlab/coaxial_ion_loss_criteria.m
@@ -0,0 +1,230 @@
+%% wr+
+H0=3.2e-14;
+P0=8.66e-25;
+r0=0.005;
+
+qe=1.60217662E-19;
+me=9.109383E-31;
+eps_0=8.85418781762e-12;
+vlight=299792458;
+
+ne=5e15;
+ne=5*logspace(10,19,5000);
+
+kb=1.38e-23;
+T=300;
+P=1e-9*100000;
+estimneutraldensity=P/(kb*T)
+
+I=40;
+V=90;
+alpha=1.3;
+
+
+B0=0.21;
+L=0.48;
+z=0.0;
+
+% Davidson
+prefix='Davidson';
+phia=-60000;
+phib=0;
+R=1.5;
+B0=0.21;
+width=0.48;
+b=0.16;
+a=0.001;
+r=0.0075;
+z=0.22;
+deltar=(R-1)/(R+1)*besseli(1,2*pi*r/width)*(width/2/pi)*(cos(2*pi*z/width)+1);
+rb_=0.007;
+rbplus=0.009872;
+Te=400; %eV
+%1rbplus=0.020;
+
+% rb_=0.01;
+% rbplus=rb_+0.1;
+% r=(rb_+rbplus)/2;
+% r=0.015
+% deltar=(R-1)/(R+1)*r;
+
+
+
+% ITER
+% prefix='ITER'
+% R=1.0001;
+% phia=-10000;
+% phib=0;
+% b=0.09;
+% a=0.068;
+% r=0.085;
+% deltar=(R-1)/(R+1)*r;
+% rb_=0.075;
+% rbplus=0.088;
+
+spcharge=@(ne)-qe*ne/2/eps_0;
+%
+% Apart=@(ne)spcharge(ne)/2*(2*r*deltar+deltar^2);
+%
+% Bpart=@(ne)log(1+deltar/r)/log(b/a)*(phib-phia+spcharge(ne)*rb_^2*log(rbplus/rb_));
+%
+% Cpart=@(ne)-log(1+deltar/r)/log(rbplus/rb_)*(phib+spcharge(ne)*rb_^2*log(rbplus/rb_));
+%
+% Dpart=@(ne)spcharge(ne)*(rbplus^2-rb_^2)/(log(b/a)*log(rbplus/rb_))*log(1+deltar/r)*...
+% ( log(b/rbplus)*(log(rbplus/a)+1/2) - log(rb_/a)*(log(b/rbplus)-1/2) );
+
+deltaphicloud=@(ne) Phi(r+deltar,rb_,rbplus,a,b,phia,phib, ne)-Phi(r,rb_,rbplus,a,b,phia,phib, ne);
+deltaphivacuum=(phib-phia)/log(b/a)*log(2*R/(R+1));
+
+Emaxcloud=1/(R-1)*deltaphicloud(ne);
+Emaxvacuum=1/(R-1)*deltaphivacuum
+
+figure('name',sprintf('%s Phi_a=%3.2f kV Phi_b=%3.2f kV r=%3.2f mm R=%3.2f',prefix,phia/1e3,phib/1e3,r*1e3, R));
+loglog(ne,abs(Emaxcloud),'displayname','flat top density cloud')
+hold on
+loglog([ne(1) ne(end)], Emaxvacuum*[1 1],'displayname','Vacuum')
+xlabel('electron density [m^{-3}]')
+ylabel('|E_{min}| [eV]')
+title(sprintf('%s \\Phi_a=%3.2f kV \\Phi_b=%3.2f kV r=%3.2f mm R=%3.2f',prefix,phia/1e3,phib/1e3,r*1e3, R))
+savefig(sprintf('%s_Phia%3.2f_Phib%3.2f_r%3.2f_R%3.2f.fig',prefix,phia/1e3,phib/1e3,r*1e3, R))
+
+
+figure
+subplot(2,2,1)
+sgtitle(sprintf('%s \\Phi_a=%3.2f kV \\Phi_b=%3.2f kV r_0=%3.2f mm z_0=%3.2f mm R=%3.2f',prefix,phia/1e3,phib/1e3,(r+deltar)*1e3,z*1e3, R))
+neinter=5e13;
+vpar=linspace(1,vlight,1000);
+vper=sqrt(-2*qe/me*deltaphicloud(neinter)/(R-1)+vpar.^2/(R-1))/vlight;
+vpar=vpar/vlight;
+vpar=vpar(real(vper)~=0);
+vper=vper(real(vper)~=0);
+
+%Davidson vperp
+phi0=Phi(r,rb_,rbplus,a,b,phia,phib, neinter);
+vrz=sqrt(2/me*H0+2/me*qe*phi0-1/me^2*(P0/r+qe*Athet(r,z,B0,width,R))^2);
+
+theta=2*pi*linspace(0,1,1000);
+vthet=(P0/r+qe*Athet(r,z,B0,width,R))/me;
+vperdav=sqrt((vrz^2*cos(theta).^2+vthet^2))/vlight;
+vpardav=vrz*sin(theta)/vlight;
+
+vpargauss=sqrt(Te*qe/me)*linspace(-1,1,200);
+vpergauss=sqrt(2*Te*qe/me)*sqrt(1-vpargauss.^2*me/Te/qe);
+vpargauss=[flip(vpargauss(2:end)) vpargauss]/vlight;
+vpergauss=[-flip(vpergauss(2:end)) vpergauss]/vlight;
+
+plot(vpar,vper,'b-')
+hold on
+plot(vpar,-vper,'b-')
+plot(-vpar,vper,'b-')
+plot(-vpar,-vper,'b-')
+xlabel('\beta_z [m/s]')
+ylabel('\beta_\perp [m/s]')
+title(sprintf('electrons n_e=%3.2e m^{-3}',neinter))
+grid
+plot(vpardav,vperdav,'r--')
+plot(vpargauss,vpergauss,'g--')
+xlim([-1 1])
+ylim([-1 1])
+axis equal
+
+subplot(2,2,3)
+neinter=5e17;
+vpar=linspace(1,vlight,1000);
+vper=sqrt(-2*qe/me*deltaphicloud(neinter)/(R-1)+vpar.^2/(R-1))/vlight;
+vpar=vpar/vlight;
+vpar=vpar(real(vper)~=0);
+vper=vper(real(vper)~=0);
+
+%Davidson vperp
+phi0=Phi(r,rb_,rbplus,a,b,phia,phib, neinter);
+vrz=sqrt(2/me*H0+2/me*qe*phi0-1/me^2*(P0/r+qe*Athet(r,z,B0,width,R))^2);
+theta=2*pi*linspace(0,1,1000);
+vthet=(P0/r+qe*Athet(r,z,B0,width,R))/me;
+vperdav=sqrt((vrz^2*cos(theta).^2+vthet^2))/vlight;
+vpardav=vrz*sin(theta)/vlight;
+
+
+vpargauss=sqrt(Te*qe/me)*linspace(-1,1,200);
+vpergauss=sqrt(2*Te*qe/me)*sqrt(1-vpargauss.^2*me/Te/qe);
+vpargauss=[flip(vpargauss(2:end)) vpargauss]/vlight;
+vpergauss=[-flip(vpergauss(2:end)) vpergauss]/vlight;
+
+plot(vpar,vper,'b-')
+hold on
+plot(vpar,-vper,'b-')
+plot(-vpar,vper,'b-')
+plot(-vpar,-vper,'b-')
+xlabel('\beta_z [m/s]')
+ylabel('\beta_\perp [m/s]')
+title(sprintf('electrons n_e=%3.2e m^{-3}',neinter))
+grid
+plot(vpardav,vperdav,'r--')
+plot(vpargauss,vpergauss,'g--')
+xlim([-1 1])
+ylim([-1 1])
+axis equal
+
+subplot(2,2,2)
+neinter=5e13;
+vpar=linspace(1,vlight,1000);
+m=28.0134/1000*6.02214076e23;
+vper=sqrt(2*qe/m*deltaphicloud(neinter)/(R-1)+vpar.^2/(R-1))/vlight;
+vpar=vpar/vlight;
+vpar=vpar(real(vper)~=0);
+vper=vper(real(vper)~=0);
+plot(vpar,vper,'b-')
+hold on
+plot(vpar,-vper,'b-')
+plot(-vpar,vper,'b-')
+plot(-vpar,-vper,'b-')
+xlabel('\beta_z [m/s]')
+ylabel('\beta_\perp [m/s]')
+title(sprintf('N_2^+ ions n_e=%3.2e m^{-3}',neinter))
+grid
+xlim([-1 1])
+ylim([-1 1])
+
+subplot(2,2,4)
+neinter=5e17;
+vpar=linspace(1,vlight,1000);
+vper=sqrt(2*qe/m*deltaphicloud(neinter)/(R-1)+vpar.^2/(R-1))/vlight;
+vpar=vpar/vlight;
+vpar=vpar(real(vper)~=0);
+vper=vper(real(vper)~=0);
+plot(vpar,vper,'b-')
+hold on
+plot(vpar,-vper,'b-')
+plot(-vpar,vper,'b-')
+plot(-vpar,-vper,'b-')
+xlabel('\beta_z [m/s]')
+ylabel('\beta_\perp [m/s]')
+title(sprintf('N_2^+ ions n_e=%3.2e m^{-3}',neinter))
+grid
+savefig(sprintf('%s_Phia%3.2f_Phib%3.2f_r%3.2f_z%3.2f_R%3.2f_cone.fig',prefix,phia/1e3,phib/1e3,(r+deltar)*1e3,z*1e3, R))
+xlim([-1 1])
+ylim([-1 1])
+
+function Atheta=Athet(r,z,B0,width,R)
+ Atheta=0.5*B0*(r-width/pi*(R-1)/(R+1)...
+ .*besseli(1,2*pi*r/width).*cos(2*pi*z/width));
+end
+
+function phi=Phi(r,rbm,rbp,a,b,phia,phib, ne)
+qe=1.60217662E-19;
+me=9.109383E-31;
+eps_0=8.85418781762e-12;
+
+phibp=1/log(b/a)*(ne/2*qe/eps_0*rbm^2*log(b/rbp)*log(rbp/rbm)+(-qe*ne/2/eps_0*(rbp^2-rbm^2)*(log(rbp/a)+0.5)+phia-phib)*log(b/rbp));
+phibm=1/log(b/a)*(-qe/eps_0*ne/2*rbm^2*log(rbm/a)*log(rbp/rbm)+(-qe/eps_0*ne/2*(rbp^2-rbm^2)*(log(b/rbp)-0.5)+phib)*log(rbm/a)+phia*log(b/rbm));
+if(r<rbm && r>a)
+ phi=((phibm-phia)*log(r)+phia*log(rbm)-phibm*log(a))/log(rbm/a);
+elseif(r>=rbm && r<=rbp)
+ phi=-qe/eps_0*ne/4*(r^2-rbp^2)+phibp+(phibp-phibm+qe/eps-0*ne/4*(rbp^2-rbm^2))*log(r/rbp)/log(rbp/rbm);
+elseif(r>rbp && r< b)
+ phi=((phib-phibp)*log(r)+phibp*log(b)-phib*log(rbp))/log(b/rbp);
+else
+ error('invalid radial position')
+end
+
+end
\ No newline at end of file
diff --git a/matlab/dispespicFields.m b/matlab/dispespicFields.m
index 5d12b0c..5d481d3 100644
--- a/matlab/dispespicFields.m
+++ b/matlab/dispespicFields.m
@@ -1,192 +1,192 @@
function dispespicFields(M,logdensity,showgrid,fixed)
%UNTITLED Summary of this function goes here
% Detailed explanation goes here
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','Grid data');
+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;
for i=edt.Value:sld.Limits(2)
edt.Value=i;
sld.Value=i;
updatesubplotsdata(i,mf);
pause(0.01)
if stop
stop=false;
break;
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));
%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.Er(:,:,fieldstep))
xlabel(ax3,'z [m]')
ylabel(ax3,'r [m]')
c = colorbar(ax3);
c.Label.String= 'E_r [V/m]';
title(ax3,'Radial Electric field')
grid(ax3, 'on')
ax4=subplot(2,2,4,'Parent',fig);
contourf(ax4,M.zgrid,M.rgrid,M.Ez(:,:,fieldstep))
xlabel(ax4,'z [m]')
ylabel(ax4,'r [m]')
c = colorbar(ax4);
c.Label.String= 'E_z [V/m]';
title(ax4,'Axial Electric field')
grid(ax4, 'on')
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,double((fieldstep-1)*M.it1)*M.dt))
%% 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.Er(:,:,fieldstep);
%% update Axial electric field
fig.Children(end-6).Children(1).ZData=M.Ez(:,:,fieldstep);
drawnow limitrate
end
end
diff --git a/matlab/dispespicParts.m b/matlab/dispespicParts.m
index 4b41db5..245f13f 100644
--- a/matlab/dispespicParts.m
+++ b/matlab/dispespicParts.m
@@ -1,138 +1,138 @@
function dispespicParts(M)
%UNTITLED2 Summary of this function goes here
% Detailed explanation goes here
fieldstep=1;
f2=uifigure('Name','Particles data');
mf2=uipanel(f2,'Position',[5 50 f2.Position(3)-10 f2.Position(4)-55]);
mf2.AutoResizeChildren='off';
m2=uipanel(f2,'Position',[5 5 f2.Position(3)-10 40]);
sgtitle(mf2,sprintf('step=%d t=%0.5e s',fieldstep*M.it2,M.tpart(fieldstep)))
sld2 = uislider(m2,'Position',[10 30 0.6*m2.Position(3) 3]);
sld2.Value=fieldstep;
sld2.Limits=[1 size(M.VR,2)];
edt2 = uieditfield(m2,'numeric','Limits',[1 size(M.VR,2)],'Value',1);
edt2.Position=[sld2.Position(1)+sld2.Position(3)+25 5 40 20];
edt2.RoundFractionalValues='on';
Printbt2=uibutton(m2,'Position',[edt2.Position(1)+edt2.Position(3)+10 5 40 20],'Text', 'Save');
sld2.ValueChangingFcn={@updatefigpartdata,edt2,mf2};
edt2.ValueChangedFcn={@updatefigpartdata,sld2,mf2};
Printbt2.ButtonPushedFcn={@plotPartButtonPushed};
PlotEspic2dpartdata(mf2,M,fieldstep);
function PlotEspic2dpartdata(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.it2,M.tpart(fieldstep)))
ax1=subplot(3,2,1,'Parent',fig);
plot(ax1,M.Z(:,fieldstep),M.R(:,fieldstep),'.');
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')
grid(ax1, 'on')
ax2=subplot(3,2,2,'Parent',fig);
plot(ax2,M.VZ(:,fieldstep),M.VR(:,fieldstep),'.');
xlabel(ax2,'VZ [m/s]')
ylabel(ax2,'VR [m/s]')
axis(ax2, 'equal')
grid(ax2, 'on')
ax3=subplot(3,2,3,'Parent',fig);
plot(ax3,M.Z(:,fieldstep),M.VZ(:,fieldstep),'.');
xlim(ax3,[M.zgrid(1) M.zgrid(end)])
xlabel(ax3,'Z [m]')
ylabel(ax3,'VZ [m/s]')
grid(ax3, 'on')
ax4=subplot(3,2,4,'Parent',fig);
plot(ax4,M.VZ(:,fieldstep),M.R(:,fieldstep),'.')
ylabel(ax4,'R [m]')
xlabel(ax4,'VZ [m/s]')
ylim(ax4,[M.rgrid(1) M.rgrid(end)])
grid(ax4, 'on')
ax5=subplot(3,2,5,'Parent',fig);
plot(ax5,M.Z(:,fieldstep),M.VR(:,fieldstep),'.');
xlim(ax5,[M.zgrid(1) M.zgrid(end)])
xlabel(ax5,'Z [m]')
ylabel(ax5,'VR [m/s]')
grid(ax5, 'on')
ax6=subplot(3,2,6,'Parent',fig);
plot(ax6,M.VR(:,fieldstep),M.R(:,fieldstep),'.');
ylim(ax6,[M.rgrid(1) M.rgrid(end)])
ylabel(ax6,'R [m]')
xlabel(ax6,'VR [m/s]')
grid(ax6, 'on')
end
function updatefigpartdata(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;
- sgtitle(fig,sprintf('step=%d t=%0.5e s',(fieldstep-1)*M.it2,double((fieldstep-1)*M.it2)*M.dt))
+ sgtitle(fig,sprintf('t=%0.5e s',M.tpart(fieldstep-1)))
%% update Position plot
ax1=fig.Children(end);
ax1.Children.XData=M.Z(:,fieldstep);
ax1.Children.YData=M.R(:,fieldstep);
view(ax1,2)
%% update VZ VPERP
fig.Children(end-1).Children(1).XData=M.VZ(:,fieldstep);
fig.Children(end-1).Children(1).YData=M.VR(:,fieldstep);
axis(fig.Children(end-1),'equal')
%% update Z VZ
fig.Children(end-2).Children(1).XData=M.Z(:,fieldstep);
fig.Children(end-2).Children(1).YData=M.VZ(:,fieldstep);
%% update VZ VR
fig.Children(end-3).Children(1).XData=M.VZ(:,fieldstep);
fig.Children(end-3).Children(1).YData=M.R(:,fieldstep);
%% update R VR
fig.Children(end-4).Children(1).XData=M.Z(:,fieldstep);
fig.Children(end-4).Children(1).YData=M.VR(:,fieldstep);
%% update Z VR
fig.Children(end-5).Children(1).XData=M.VR(:,fieldstep);
fig.Children(end-5).Children(1).YData=M.R(:,fieldstep);
drawnow limitrate
end
function plotPartButtonPushed(btn,ax)
%UNTITLED2 Summary of this function goes here
% Detailed explanation goes here
f=figure();
PlotEspic2dpartdata(f,M,sld.Value);
f.PaperOrientation='portrait';
[~, name, ~] = fileparts(M.file);
print(f,sprintf('%sParts%d',name,sld.Value),'-dpdf','-fillpage')
end
end
diff --git a/matlab/distribution_function.m b/matlab/distribution_function.m
index 7732200..9cdc897 100644
--- a/matlab/distribution_function.m
+++ b/matlab/distribution_function.m
@@ -1,84 +1,83 @@
fig=figure;
ax1=gca;
Ndistrib=M.N(:,:,end);
-pcolor(ax1,M.zgrid,M.rgrid,Ndistrib);
+h=pcolor(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')
stepmin=2000;
stepend=size(M.t2d,1);
papsize=[16 14];
[x,y]=ginput(1);
Zindex=find(x>M.zgrid,1,'last');
Rindex=find(y>M.rgrid,1,'last');
Z=(M.zgrid(Zindex));
R=(M.rgrid(Rindex));
plot(Z,R,'rx','Markersize',12);
% Rindex=16; Zindex=64;
% Rindex=28; Zindex=floor(size(M.zgrid,1)/2)+1;
Indices=M.Rindex(:,1)==Rindex & M.Zindex(:,1)==Zindex;
Vr=M.VR(Indices,1)/M.vlight;
Vz=M.VZ(Indices,1)/M.vlight;
Vthet=M.VTHET(Indices,1)/M.vlight;
Indicesend=M.Rindex(:,end)==Rindex & M.Zindex(:,end)==Zindex;
Vrend=M.VR(Indicesend,end)/M.vlight;
Vzend=M.VZ(Indicesend,end)/M.vlight;
Vthetend=M.VTHET(Indicesend,end)/M.vlight;
-binwidth=0.05;
+binwidth=abs(max(Vrend)-min(Vrend))/10;
f=figure('Name',sprintf("%s v distrib",M.file));
tstudied=0;
legtext=sprintf("t=%2.1f [ns]",M.tpart(end)*1e9);
subplot(1,3,1)
h1=histogram(Vr,'Binwidth',binwidth,'Normalization','probability','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9));
hold on
h1=histogram(Vrend,'Binwidth',binwidth,'Normalization','probability','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9));
ylabel('probability')
xlabel('\beta_r')
grid on
ax2=subplot(1,3,2);
-binwidth=0.01;
h1=histogram(Vthet(:,1),'Binwidth',binwidth,'Normalization','probability','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9));
hold on
h1=histogram(Vthetend,'Binwidth',binwidth,'Normalization','probability','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9));
ylabel('probability')
xlabel('\beta_\theta')
legend(ax2,'location','northoutside','orientation','vertical')
grid on
subplot(1,3,3)
-binwidth=0.05;
h1=histogram(Vz(:,1),'Binwidth',binwidth,'Normalization','probability','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9));
hold on
h1=histogram(Vzend,'Binwidth',binwidth,'Normalization','probability','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9));
ylabel('probability')
xlabel('\beta_z')
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))
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))
\ No newline at end of file
diff --git a/matlab/espic2dhdf5.m b/matlab/espic2dhdf5.m
index 773d635..ffdf7a8 100644
--- a/matlab/espic2dhdf5.m
+++ b/matlab/espic2dhdf5.m
@@ -1,409 +1,499 @@
classdef espic2dhdf5
%UNTITLED Summary of this class goes here
% Detailed explanation goes here
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
nrun
nlres
nlsave
nlclassical
nlPhis
nz
nnr
lz
- qsim
- msim
nplasma
potinn
potout
B0
Rcurv
width
n0
temp
femorder
ngauss
plasmadim
radii
H0
P0
+ normalized
+ nbspecies
%% Frequencies
omepe
omece
%% Normalizations
tnorm
rnorm
bnorm
enorm
phinorm
vnorm
%% Grid data
rgrid
zgrid
dz
dr
Vcell
%% Magnetic field
Br
Bz
Athet
B
%% Energies
epot
ekin
etot
etot0
eerr
%% 2D time data
N
fluidUR
fluidUZ
fluidUTHET
pot
Er
Ez
CellVol
Presstens
%% Splines
knotsr
knotsz
%% Particle parameters
+ weight
+ qsim
+ msim
nbparts
H
P
partepot
R
Z
Rindex
Zindex
partindex
VR
VZ
VTHET
THET
+ species
end
methods
function file=file(obj)
file=obj.filename;
end
function obj = espic2dhdf5(filename,readparts)
filedata=dir(filename);
- obj.filename=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);
- obj.dt = h5readatt(obj.filename,'/data/input.00/','dt');
- obj.nrun = h5readatt(obj.filename,'/data/input.00/','nrun');
- obj.nlres = strcmp(h5readatt(obj.filename,'/data/input.00/','nlres'),'y');
- obj.nlsave = strcmp(h5readatt(obj.filename,'/data/input.00/','nlsave'),'y');
- obj.nlclassical =strcmp(h5readatt(obj.filename,'/data/input.00/','nlclassical'),'y');
- obj.nlPhis =strcmp(h5readatt(obj.filename,'/data/input.00/','nlPhis'),'y');
- obj.nz = h5readatt(obj.filename,'/data/input.00/','nz');
- obj.nnr = h5read(obj.filename,'/data/input.00/nnr');
- obj.lz = h5read(obj.filename,'/data/input.00/lz');
- obj.qsim = h5readatt(obj.filename,'/data/input.00/','qsim');
- obj.msim = h5readatt(obj.filename,'/data/input.00/','msim');
- obj.nplasma = h5readatt(obj.filename,'/data/input.00/','nplasma');
- obj.potinn = h5readatt(obj.filename,'/data/input.00/','potinn');
- obj.potout = h5readatt(obj.filename,'/data/input.00/','potout');
- obj.B0 = h5readatt(obj.filename,'/data/input.00/','B0');
- obj.Rcurv = h5readatt(obj.filename,'/data/input.00/','Rcurv');
- obj.width = h5readatt(obj.filename,'/data/input.00/','width');
- obj.n0 = h5readatt(obj.filename,'/data/input.00/','n0');
- obj.temp = h5readatt(obj.filename,'/data/input.00/','temp');
+ 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.it0 = h5readatt(obj.filename,'/data/input.00/','it0d');
- obj.it1 = h5readatt(obj.filename,'/data/input.00/','it2d');
- obj.it2 = h5readatt(obj.filename,'/data/input.00/','itparts');
+ obj.nbspecies=h5readatt(obj.fullpath,'/data/input.00/','nbspecies');
+ obj.normalized=strcmp(h5readatt(obj.fullpath,'/data/input.00/','rawparts'),'y');
catch
- obj.it0 = h5readatt(obj.filename,'/data/input.00/','it0');
- obj.it1 = h5readatt(obj.filename,'/data/input.00/','it1');
- obj.it1 = h5readatt(obj.filename,'/data/input.00/','it2');
+ obj.nbspecies=1;
+ obj.normalized=false;
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.filename, '/data/var0d/nbparts');
+ obj.nbparts= h5read(obj.fullpath, '/data/var0d/nbparts');
try
- obj.H0 = h5read(obj.filename,'/data/input.00/H0');
- obj.P0 = h5read(obj.filename,'/data/input.00/P0');
+ 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.filename, '/data/var1d/rgrid')*obj.rnorm;
- obj.zgrid= h5read(obj.filename, '/data/var1d/zgrid')*obj.rnorm;
+ 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
obj.Vcell=(obj.zgrid(2:end)-obj.zgrid(1:end-1))*((obj.rgrid(2:end).^2-obj.rgrid(1:end-1).^2)*pi)';
- Br = h5read(obj.filename,'/data/fields/Br')*obj.bnorm;
+ Br = h5read(obj.fullpath,'/data/fields/Br')*obj.bnorm;
obj.Br= reshape(Br,length(obj.zgrid),length(obj.rgrid));
- Bz = h5read(obj.filename,'/data/fields/Bz')*obj.bnorm;
+ Bz = h5read(obj.fullpath,'/data/fields/Bz')*obj.bnorm;
obj.Bz= reshape(Bz,length(obj.zgrid),length(obj.rgrid));
try
- Atheta = h5read(obj.filename,'/data/fields/Athet')*obj.bnorm;
+ Atheta = h5read(obj.fullpath,'/data/fields/Athet')*obj.bnorm;
obj.Athet= reshape(Atheta,length(obj.zgrid),length(obj.rgrid));
catch
end
obj.B=sqrt(obj.Bz.^2+obj.Br.^2);
clear Br Bz
try
- obj.t0d=h5read(obj.filename,'/data/var0d/time');
+ obj.t0d=h5read(obj.fullpath,'/data/var0d/time');
catch
obj.t0d=obj.dt.*double(0:length(obj.epot)-1);
end
- obj.femorder = h5read(obj.filename,'/data/input.00/femorder');
- obj.ngauss = h5read(obj.filename,'/data/input.00/ngauss');
- obj.plasmadim = h5read(obj.filename,'/data/input.00/plasmadim');
- obj.radii = h5read(obj.filename,'/data/input.00/radii');
+ 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.filename,'/data/var0d/epot');
- obj.ekin = h5read(obj.filename,'/data/var0d/ekin');
- obj.etot = h5read(obj.filename,'/data/var0d/etot');
+ 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(filename,'/data/var0d/etot0');
+ obj.etot0 = h5read(obj.fullpath,'/data/var0d/etot0');
obj.eerr = obj.etot-obj.etot0;
catch
obj.eerr = obj.etot-obj.etot(2);
end
- obj.pot=gridquantity(obj.filename,'/data/fields/pot',sum(obj.nnr)+1, obj.nz+1,obj.phinorm);
- obj.Er=gridquantity(obj.filename,'/data/fields/Er',sum(obj.nnr)+1, obj.nz+1,obj.enorm);
- obj.Ez=gridquantity(obj.filename,'/data/fields/Ez',sum(obj.nnr)+1, obj.nz+1,obj.enorm);
+ 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.filename,'/data/fields/time');
+ obj.t2d = h5read(obj.fullpath,'/data/fields/time');
catch
- info=h5info(obj.filename,'/data/fields/partdensity');
+ 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.filename,'/data/fields/moments');
- obj.femorder = h5read(obj.filename,'/data/input.00/femorder');
+ 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.filename,'/data/fields/volume'),length(obj.knotsz)-kz,length(obj.knotsr)-kr);
+ 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
-
- obj.N=splinedensity(obj.filename, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, abs(obj.qsim/obj.qe), 1);
- obj.fluidUR=splinevelocity(obj.filename, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.vnorm, 2);
- obj.fluidUTHET=splinevelocity(obj.filename, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.vnorm, 3);
- obj.fluidUZ=splinevelocity(obj.filename, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.vnorm, 4);
- obj.Presstens=splinepressure(obj.filename, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, obj.vnorm^2*obj.msim);
+ 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.filename, '/data/fields/partdensity', sum(obj.nnr)+1, obj.nz+1, obj.CellVol, abs(obj.qsim/obj.qe), true);
- obj.fluidUR=gridquantity(obj.filename, '/data/fields/fluidur', sum(obj.nnr)+1, obj.nz+1, obj.vnorm, true);
- obj.fluidUTHET=gridquantity(obj.filename, '/data/fields/fluiduthet', sum(obj.nnr)+1, obj.nz+1, obj.vnorm, true);
- obj.fluidUZ=gridquantity(obj.filename, '/data/fields/fluiduz', sum(obj.nnr)+1, obj.nz+1, obj.vnorm, true);
+ 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)
-
- obj.R = h5read(filename,'/data/part/R')*obj.rnorm;
- obj.Z = h5read(filename,'/data/part/Z')*obj.rnorm;
- try
- obj.THET = h5read(filename,'/data/part/THET');
- catch
- clear obj.THET
- end
- try
- obj.Rindex=h5read(filename,'/data/part/Rindex');
- obj.Zindex=h5read(filename,'/data/part/Zindex');
- catch
- clear obj.Rindex obj.Zindex
- end
- UR = h5read(filename,'/data/part/UR');
- UZ = h5read(filename,'/data/part/UZ');
- UTHET= h5read(filename,'/data/part/UTHET');
- if(obj.nlclassical)
- obj.VR = UR*obj.vnorm;
- clear UR
- obj.VZ = UZ*obj.vnorm;
- clear UZ
- obj.VTHET = UTHET*obj.vnorm;
- clear UTHET
- else
- GAMM = sqrt(1+UR.^2+UZ.^2+UTHET.^2);
- obj.VR = UR./GAMM*obj.vnorm;
- clear UR
- obj.VZ = UZ./GAMM*obj.vnorm;
- clear UZ
- obj.VTHET = UTHET./GAMM*obj.vnorm;
- clear UTHET
- clear GAMM
- end
- %obj.VPERP = sqrt(obj.VR.*obj.VR+obj.VTHET.*obj.VTHET);
- obj.partepot = sign(obj.qsim)*h5read(filename,'/data/part/pot')*obj.phinorm*obj.qe;
- obj.H=0.5*obj.me*(obj.VR.^2+obj.VTHET.^2+obj.VZ.^2)+obj.partepot;
- %clear partepot
- obj.P=obj.R.*(obj.VTHET*obj.me+sign(obj.qsim)*obj.qe*Athet(obj.R,obj.Z));
- try
- obj.partindex = h5read(filename,'/data/part/partindex');
- partindex=obj.partindex;
- [~,Indices]=sort(partindex,'descend');
- for i=1:size(Indices,2)
- obj.H(:,i)=obj.H(Indices(:,i),i);
- obj.P(:,i)=obj.P(Indices(:,i),i);
- obj.R(:,i)=obj.R(Indices(:,i),i);
- obj.Z(:,i)=obj.Z(Indices(:,i),i);
- obj.Rindex(:,i)=obj.Rindex(Indices(:,i),i);
- obj.Zindex(:,i)=obj.Zindex(Indices(:,i),i);
- obj.VR(:,i)=obj.VR(Indices(:,i),i);
- obj.VZ(:,i)=obj.VZ(Indices(:,i),i);
- obj.VTHET(:,i)=obj.VTHET(Indices(:,i),i);
- obj.THET(:,i)=obj.THET(Indices(:,i),i);
+
+ if(obj.normalized)
+ obj.R = h5read(obj.fullpath,'/data/part/R');
+ obj.Z = h5read(obj.fullpath,'/data/part/Z');
+ else
+ obj.R = h5read(obj.fullpath,'/data/part/R')*obj.rnorm;
+ obj.Z = h5read(obj.fullpath,'/data/part/Z')*obj.rnorm;
end
- catch
+ try
+ obj.THET = h5read(obj.fullpath,'/data/part/THET');
+ catch
+ clear obj.THET
+ end
+ try
+ obj.Rindex=h5read(obj.fullpath,'/data/part/Rindex');
+ obj.Zindex=h5read(obj.fullpath,'/data/part/Zindex');
+ catch
+ clear obj.Rindex obj.Zindex
+ end
+ UR = h5read(obj.fullpath,'/data/part/UR');
+ UZ = h5read(obj.fullpath,'/data/part/UZ');
+ UTHET= h5read(obj.fullpath,'/data/part/UTHET');
+ if(obj.nlclassical)
+ obj.VR = UR*obj.vnorm;
+ clear UR
+ obj.VZ = UZ*obj.vnorm;
+ clear UZ
+ obj.VTHET = UTHET*obj.vnorm;
+ clear UTHET
+ else
+ GAMM = sqrt(1+UR.^2+UZ.^2+UTHET.^2);
+ obj.VR = UR./GAMM*obj.vnorm;
+ clear UR
+ obj.VZ = UZ./GAMM*obj.vnorm;
+ clear UZ
+ obj.VTHET = UTHET./GAMM*obj.vnorm;
+ clear UTHET
+ clear GAMM
+ end
+ %obj.VPERP = sqrt(obj.VR.*obj.VR+obj.VTHET.*obj.VTHET);
+ if(obj.normalized)
+ obj.partepot = sign(obj.qsim)*h5read(filename,'/data/part/pot')*obj.qe;
+ else
+ obj.partepot = sign(obj.qsim)*h5read(filename,'/data/part/pot')*obj.phinorm*obj.qe;
+ end
+ obj.H=0.5*obj.me*(obj.VR.^2+obj.VTHET.^2+obj.VZ.^2)+obj.partepot;
+ %clear partepot
+ obj.P=obj.R.*(obj.VTHET*obj.me+sign(obj.qsim)*obj.qe*Athet(obj.R,obj.Z));
+ try
+ obj.partindex = h5read(obj.fullpath,'/data/part/partindex');
+ partindex=obj.partindex;
+ [~,Indices]=sort(partindex,'descend');
+ for i=1:size(Indices,2)
+ obj.H(:,i)=obj.H(Indices(:,i),i);
+ obj.P(:,i)=obj.P(Indices(:,i),i);
+ obj.R(:,i)=obj.R(Indices(:,i),i);
+ obj.Z(:,i)=obj.Z(Indices(:,i),i);
+ obj.Rindex(:,i)=obj.Rindex(Indices(:,i),i);
+ obj.Zindex(:,i)=obj.Zindex(Indices(:,i),i);
+ obj.VR(:,i)=obj.VR(Indices(:,i),i);
+ obj.VZ(:,i)=obj.VZ(Indices(:,i),i);
+ obj.VTHET(:,i)=obj.VTHET(Indices(:,i),i);
+ obj.THET(:,i)=obj.THET(Indices(:,i),i);
+ end
+ catch
+ end
+ clear partindex;
end
- clear partindex;
- end
try
- obj.tpart = h5read(filename,'/data/part/time');
+ obj.tpart = h5read(obj.fullpath,'/data/part/time');
catch
obj.tpart=obj.dt*(0:size(obj.R,2)-1)*double(obj.it2);
end
function Atheta=Athet(R,Z)
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/obj.width));
+ .*besseli(1,2*pi*R/obj.width).*cos(2*pi*Z/obj.width));
+ end
+ if(obj.nbspecies >1)
+ obj.species=h5parts.empty;
+ for i=2:obj.nbspecies
+ obj.species(i-1)=h5parts(obj.fullpath,sprintf('/data/part/%2d',i),obj);
+ end
end
end
function changed=ischanged(obj)
%ischanged Check if the file has been changed and some data
%must be reloaded
try
- filedata=dir(obj.filename);
+ 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 plotEnergy(obj)
tmin=2;
tmax=length(obj.ekin);
- [~, name, ~] = fileparts(obj.file);
- f=figure('Name', sprintf('%s Energy',name));
+ 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--')
+ 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('%sEnergy',name));
+ obj.savegraph(f,sprintf('%s/%sEnergy',obj.folder,obj.name));
+ end
+
+ function trappedParticles(obj)
+ 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 savegraph(obj, figure, name, papsize)
figure.PaperUnits='centimeters';
if (nargin < 4)
papsize=[14 16];
end
figure.PaperSize=papsize;
print(figure,name,'-dpdf','-fillpage')
savefig(figure,name)
end
+
+ function displayHP(obj,tstart)
+ if(size(obj.H,2)>=2)
+ 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),size(obj.H,1));
+ H=obj.H(1:obj.nbparts(1),1);
+ h1=histogram(H,10,'BinLimits',[min(H(:)) max(H(:))],'DisplayName',sprintf("t=%2.3d [ns]",obj.tpart(1)*1e9));
+ hold on
+ H=mean(obj.H(1:partsmax,tstart:end),2);
+ %,'Binwidth',h1.BinWidth
+ h1=histogram(H,20,'BinLimits',[min(H(:)) max(H(:))],'DisplayName',legtext);
+ ylabel('counts')
+ xlabel('H [J]')
+ legend
+
+ subplot(1,2,2)
+ P=obj.P(1:obj.nbparts(1),1);
+ h2=histogram(P,10,'BinLimits',[min(P(:)) max(P(:))],'DisplayName',sprintf("t=%2.3d [ns]",obj.tpart(1)*1e9));
+ hold on
+ P=mean(obj.P(1:partsmax,tstart:end),2);
+ histogram(P,50,'BinLimits',[min(P(:)) max(P(:))],'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
end
diff --git a/matlab/h5parts.m b/matlab/h5parts.m
new file mode 100644
index 0000000..323e3bf
--- /dev/null
+++ b/matlab/h5parts.m
@@ -0,0 +1,98 @@
+classdef h5parts
+ properties
+ fullpath
+ R
+ Z
+ THET
+ VR
+ VZ
+ VTHET
+ Rindex
+ Zindex
+ partindex
+ H
+ P
+ partepot
+ nlclassical
+ q
+ m
+ weight
+ tpart
+ it2
+ nbparts
+
+ rgrid
+ zgrid
+ vnorm
+
+ end
+ methods
+ function obj=h5parts(filename,hdf5group,parent)
+ obj.fullpath=filename;
+ obj.nlclassical=parent.nlclassical;
+ obj.it2=parent.it2;
+ obj.rgrid=parent.rgrid;
+ obj.zgrid=parent.zgrid;
+ obj.vnorm=parent.vnorm;
+ obj.R = h5read(obj.fullpath,sprintf('%s/R',hdf5group));
+ obj.Z = h5read(obj.fullpath,sprintf('%s/Z',hdf5group));
+ obj.THET = h5read(obj.fullpath,sprintf('%s/THET',hdf5group));
+ obj.Rindex=h5read(obj.fullpath,sprintf('%s/Rindex',hdf5group));
+ obj.Zindex=h5read(obj.fullpath,sprintf('%s/Zindex',hdf5group));
+ UR = h5read(obj.fullpath,sprintf('%s/UR',hdf5group));
+ UZ = h5read(obj.fullpath,sprintf('%s/UZ',hdf5group));
+ UTHET= h5read(obj.fullpath,sprintf('%s/UTHET',hdf5group));
+ obj.q= h5readatt(obj.fullpath,hdf5group,'q');
+ obj.m= h5readatt(obj.fullpath,hdf5group,'m');
+ obj.weight= h5readatt(obj.fullpath,hdf5group,'weight');
+
+ if(obj.nlclassical)
+ obj.VR = UR*obj.vnorm;
+ clear UR
+ obj.VZ = UZ*obj.vnorm;
+ clear UZ
+ obj.VTHET = UTHET*obj.vnorm;
+ clear UTHET
+ else
+ GAMM = sqrt(1+UR.^2+UZ.^2+UTHET.^2);
+ obj.VR = UR./GAMM*obj.vnorm;
+ clear UR
+ obj.VZ = UZ./GAMM*obj.vnorm;
+ clear UZ
+ obj.VTHET = UTHET./GAMM*obj.vnorm;
+ clear UTHET
+ clear GAMM
+ end
+ obj.partepot = h5read(filename,sprintf('%s/pot',hdf5group))*obj.q;
+ obj.H=0.5*obj.m*(obj.VR.^2+obj.VTHET.^2+obj.VZ.^2)+obj.partepot;
+ %clear partepot
+ obj.P=obj.R.*(obj.VTHET*obj.m+obj.q*Athet(obj.R,obj.Z));
+ try
+ obj.partindex = h5read(obj.fullpath,sprintf('%s/partindex',hdf5group));
+ partindex=obj.partindex;
+ [~,Indices]=sort(partindex,'descend');
+ for i=1:size(Indices,2)
+ obj.H(:,i)=obj.H(Indices(:,i),i);
+ obj.P(:,i)=obj.P(Indices(:,i),i);
+ obj.R(:,i)=obj.R(Indices(:,i),i);
+ obj.Z(:,i)=obj.Z(Indices(:,i),i);
+ obj.Rindex(:,i)=obj.Rindex(Indices(:,i),i);
+ obj.Zindex(:,i)=obj.Zindex(Indices(:,i),i);
+ obj.VR(:,i)=obj.VR(Indices(:,i),i);
+ obj.VZ(:,i)=obj.VZ(Indices(:,i),i);
+ obj.VTHET(:,i)=obj.VTHET(Indices(:,i),i);
+ obj.THET(:,i)=obj.THET(Indices(:,i),i);
+ end
+ catch
+ end
+ clear partindex;
+ obj.tpart = h5read(obj.fullpath,sprintf('%s/time',hdf5group));
+ obj.nbparts = h5read(obj.fullpath,sprintf('%s/Nparts',hdf5group));
+
+ function Atheta=Athet(R,Z)
+ Atheta=0.5*parent.B0*(R-parent.width/pi*(parent.Rcurv-1)/(parent.Rcurv+1)...
+ .*besseli(1,2*pi*R/parent.width).*cos(2*pi*Z/parent.width));
+ end
+ end
+ end
+end
\ No newline at end of file
diff --git a/matlab/rbfinder.m b/matlab/rbfinder.m
new file mode 100644
index 0000000..ae36ec7
--- /dev/null
+++ b/matlab/rbfinder.m
@@ -0,0 +1,150 @@
+%% wr+
+H0=3.2e-14;
+P0=8.66e-25;
+r0=0.005;
+consideredne=5e16;
+
+%% wr-
+% H0=2.0e-16;
+% P0=-5.7e-25;
+% r0=0.0025;
+
+%% simulated
+% 5e14 5e15 1e16 5e16 1e17 2e17 3e17 4e17 1e18 3e18 5e18 1e19
+r_sim= [0.00497, 0.0053, 0.0056, 0.0059, 0.00628, 0.0050, 0.0053, 0.0052, 0.0054, 0.0052, 0.0054];%, 0.0063];
+rplussim=[0.0124, 0.0118, 0.0118, 0.0118, 0.0121, 0.0137, 0.0140, 0.016, 0.0163, 0.020, 0.025];%, 0.042];
+nsim= [4.8e14, 4.12e15, 7.5e15, 2.54e16, 4.19e16, 6.42e16, 8.55e16, 9.79e16, 1.34e17, 2.408e17, 2.416e17];%, 3.49e17];
+
+
+B0=0.21;
+Rcurv=1.5;
+L=0.48;
+z=0.0;
+
+qe=1.60217662E-19;
+me=9.109383E-31;
+eps0=8.85418781762e-12;
+
+ne=5*logspace(14,19,800);
+rb_=zeros(size(ne));
+rbplus=zeros(size(ne));
+i=1;
+minus=@(r) 1-(P0/r+qe*0.5*B0*(r-L/pi*(Rcurv-1)/(Rcurv+1)*besseli(1,2*pi*r/L)*cos(2*pi*z/L)))^2/(2*me*H0);
+%minus=@(r) 1+qe*((phib-phia)*log(r)+phia*log(b)-phib*log(a))/(log(b/a))/H0-(P0/r+qe*0.5*B0*(r-L/pi*(Rcurv-1)/(Rcurv+1)*besseli(1,2*pi*r/L)*cos(2*pi*z/L)))^2/(2*me*H0);
+format long;
+rb_(i)=fzero(minus,r0);
+plus=@(r) qe^2*ne(i)*rb_(i)/(eps0*H0)*(r-rb_(i)-rb_(i)*log(r/rb_(i)))+1-1/(2*me*H0)*(P0/r+qe*0.5*B0*(r-L/pi*(Rcurv-1)/(Rcurv+1)*besseli(1,2*pi*r/L)*cos(2*pi*z/L)))^2;
+%plus=@(r) 1+qe*((phib-phia)*log(r)+phia*log(b)-phib*log(a))/(log(b/a))/H0-(P0/r+qe*0.5*B0*(r-L/pi*(Rcurv-1)/(Rcurv+1)*besseli(1,2*pi*r/L)*cos(2*pi*z/L)))^2/(2*me*H0);
+rbplus(i)=fzero(plus,5*rb_(i));
+remainder=ones(size(ne));
+remainder(i)=plus(rbplus(i));
+for i=2:length(ne)
+ try
+ minus=@(r) 1-(P0/r+qe*0.5*B0*(r-L/pi*(Rcurv-1)/(Rcurv+1)*besseli(1,2*pi*r/L)*cos(2*pi*z/L)))^2/(2*me*H0);
+ %minus=@(r) 1+qe*((phib-phia)*log(r)+phia*log(b)-phib*log(a))/(log(b/a))/H0-(P0/r+qe*0.5*B0*(r-L/pi*(Rcurv-1)/(Rcurv+1)*besseli(1,2*pi*r/L)*cos(2*pi*z/L)))^2/(2*me*H0);
+ rb_(i)=fzero(minus,0.004);
+ rb_t=rb_(i);
+ plus=@(r) qe^2*ne(i)*rb_t/(eps0*H0)*(r-rb_t-rb_t*log(r/rb_t))+1-1/(2*me*H0)*(P0/r+qe*0.5*B0*(r-L/pi*(Rcurv-1)/(Rcurv+1)*besseli(1,2*pi*r/L)*cos(2*pi*z/L)))^2;
+ %plus=@(r) 1+qe*((phib-phia)*log(r)+phia*log(b)-phib*log(a))/(log(b/a))/H0-(P0/r+qe*0.5*B0*(r-L/pi*(Rcurv-1)/(Rcurv+1)*besseli(1,2*pi*r/L)*cos(2*pi*z/L)))^2/(2*me*H0);
+ rbplus(i)=fzero(plus,1.5*rbplus(i-1));
+ remainder(i)=plus(rbplus(i));
+ catch
+ end
+end
+
+
+f=figure;
+%subplot(2,1,1)
+semilogx(ne(remainder~=1),rbplus(remainder~=1)*1e3,'-','displayname','r_b^+')
+hold on
+semilogx(ne,rb_*1000,'-','displayname','r_b^-')
+semilogx(nsim,rplussim*1e3,'+--','displayname','r_b^+ simulated','linewidth',1.5)
+semilogx(nsim,r_sim*1e3,'+--','displayname','r_b^- simulated','linewidth',1.5)
+legend
+ylabel('r_b [mm]')
+ylim([0 160])
+xlabel('n_e(r_b^-) [m^{-3}]')
+% subplot(2,1,2)
+% semilogx(ne(remainder~=1),remainder(remainder~=1),'x--')
+% ylabel('remainder [a.u.]')
+% xlabel('n_e(r_b^-) [m^{-3}]')
+
+f.PaperOrientation='landscape';
+f.PaperUnits='centimeters';
+papsize=[10 10];
+f.PaperSize=papsize;
+print(f,'rbne','-dpdf','-fillpage')
+
+
+i=find(consideredne>=ne,1,'last');
+z=linspace(0,L/2,1000);
+rb_z=zeros(size(z));
+rbplusz=rb_z;
+minus=@(r) 1-(P0/r+qe*0.5*B0*(r-L/pi*(Rcurv-1)/(Rcurv+1)*besseli(1,2*pi*r/L)*cos(2*pi*z(1)/L)))^2/(2*me*H0);
+rb_z(1)=fzero(minus,rb_(i));
+rb_t=rb_z(1);
+z_t=z(1);
+plus=@(r) qe^2*ne(i)*rb_t/(eps0*H0)*(r-rb_t-rb_t*log(r/rb_t))+1-1/(2*me*H0)*(P0/r+qe*0.5*B0*(r-L/pi*(Rcurv-1)/(Rcurv+1)*besseli(1,2*pi*r/L)*cos(2*pi*z_t/L)))^2;
+rbplusz(1)=fzero(plus,rbplus(i));
+remainder(1)=plus(rbplusz(1));
+for j=2:length(z)
+minus=@(r) 1-(P0/r+qe*0.5*B0*(r-L/pi*(Rcurv-1)/(Rcurv+1)*besseli(1,2*pi*r/L)*cos(2*pi*z(j)/L)))^2/(2*me*H0);
+format long;
+try
+rb_z(j)=fzero(minus,rb_z(j-1));
+rb_t=rb_z(j);
+z_t=z(j);
+plus=@(r) qe^2*ne(i)*rb_t/(eps0*H0)*(r-rb_t-rb_t*log(r/rb_t))+1-1/(2*me*H0)*(P0/r+qe*0.5*B0*(r-L/pi*(Rcurv-1)/(Rcurv+1)*besseli(1,2*pi*r/L)*cos(2*pi*z_t/L)))^2;
+rbplusz(j)=fzero(plus,rbplusz(j-1));
+remainder(j)=plus(rbplusz(j));
+catch
+end
+end
+figure()
+z=[-flip(z(2:end)) z];
+rb_z=[flip(rb_z(2:end)) rb_z];
+rbplusz=[flip(rbplusz(2:end)) rbplusz];
+plot(z(rb_z~=0),rb_z(rb_z~=0),'b')
+hold on
+plot(z(rbplusz~=0),rbplusz(rbplusz~=0),'r')
+ylabel('r_b [m]')
+ylim([0 0.16])
+xlim([-L/2 L/2])
+xlabel('z [m]')
+title(sprintf('n_e=%7.2e [m^{-3}]',ne(i)))
+
+ind=rbplusz~=0 & rb_z~=0& ~isnan(rbplusz) & ~isnan(rb_z);
+ra=rb_z(ind);
+rb=rbplusz(ind);
+zvol=z(ind);
+ra=(ra(1:end-1)+ra(2:end))/2;
+rb=(rb(1:end-1)+rb(2:end))/2;
+VOL=2*pi*diff(zvol).*min(ra).*(rb-ra);
+
+nplasma=2116800
+Nbesl=ne(i)*VOL;
+Nbe=sum(Nbesl)
+Nmacrosl=floor(Nbesl/Nbe*nplasma);
+Nmacro=sum(Nmacrosl)
+remain=nplasma-Nmacro
+Nmacrosl(ceil(length(Nmacrosl)/2-remain/2):floor(length(Nmacrosl)/2+remain/2))=Nmacrosl(ceil(length(Nmacrosl)/2-remain/2):floor(length(Nmacrosl)/2+remain/2))+1;
+Nmacrosl(floor(length(Nmacrosl)/2))=Nmacrosl(floor(length(Nmacrosl)/2))+nplasma-sum(Nmacrosl);
+Nmacro=sum(Nmacrosl)
+remain=nplasma-Nmacro
+
+
+msim=me*Nbe/nplasma
+qsim=qe*Nbe/nplasma
+
+datas.ra=ra;
+datas.rb=rb;
+datas.z=zvol;
+datas.npartsslice=Nmacrosl;
+datas.m=me;
+datas.weight=sum(Nbesl)/nplasma;
+datas.q=-qe;
+datas.H0=H0;
+datas.P0=P0;
+datas.temperature=10000;
+datas.velocitytype=2;
+datas.radialtype=1;
\ No newline at end of file
diff --git a/matlab/rbfindercoax.m b/matlab/rbfindercoax.m
new file mode 100644
index 0000000..b913f36
--- /dev/null
+++ b/matlab/rbfindercoax.m
@@ -0,0 +1,147 @@
+%% wr+
+H0=3.2e-14;
+P0=8.66e-25;
+r0=0.005;
+ne=5e14;
+
+B0=0.21;
+Rcurv=1.5;
+L=0.48;
+z=0.0;
+consideredphia=-60000;
+
+qe=1.60217662E-19;
+me=9.109383E-31;
+eps0=8.85418781762e-12;
+phia=-linspace(0,90000,80);
+phib=0;
+b=0.16;
+a=0.001;
+
+%% simulated
+
+r_sim= [5.8e-3, 5.4e-3];%, 0.0063];
+rplussim=[1.11e-2, 1.20e-2];%, 0.042];
+phiasim= [-45000, -30000];%, 3.49e17];
+
+
+
+rb_phi=zeros(size(phia));
+rbplusphi=zeros(size(phia));
+i=1;
+minus=@(r) 1+qe*((phib-phia(i))*log(r)+phia(i)*log(b)-phib*log(a))/(log(b/a))/H0-(P0/r+qe*0.5*B0*(r-L/pi*(Rcurv-1)/(Rcurv+1)*besseli(1,2*pi*r/L)*cos(2*pi*z/L)))^2/(2*me*H0);
+format long;
+rb_phi(i)=fzero(minus,r0);
+%plus=@(r) qe^2*ne(i)*rb_(i)/(eps0*H0)*(r-rb_(i)-rb_(i)*log(r/rb_(i)))+1-1/(2*me*H0)*(P0/r+qe*0.5*B0*(r-L/pi*(Rcurv-1)/(Rcurv+1)*besseli(1,2*pi*r/L)*cos(2*pi*z/L)))^2;
+plus=@(r) 1+qe*((phib-phia(i))*log(r)+phia(i)*log(b)-phib*log(a))/(log(b/a))/H0-(P0/r+qe*0.5*B0*(r-L/pi*(Rcurv-1)/(Rcurv+1)*besseli(1,2*pi*r/L)*cos(2*pi*z/L)))^2/(2*me*H0);
+rbplusphi(i)=fzero(plus,5*rb_phi(i));
+remainderphi=ones(size(phia));
+remainderphi(i)=plus(rbplusphi(i));
+for i=2:length(phia)
+ try
+ minus=@(r) 1+qe*((phib-phia(i))*log(r)+phia(i)*log(b)-phib*log(a))/(log(b/a))/H0-(P0/r+qe*0.5*B0*(r-L/pi*(Rcurv-1)/(Rcurv+1)*besseli(1,2*pi*r/L)*cos(2*pi*z/L)))^2/(2*me*H0);
+ rb_phi(i)=fzero(minus,rb_phi(i-1));
+ rb_t=rb_phi(i);
+ plus=@(r) 1+qe*((phib-phia(i))*log(r)+phia(i)*log(b)-phib*log(a))/(log(b/a))/H0-(P0/r+qe*0.5*B0*(r-L/pi*(Rcurv-1)/(Rcurv+1)*besseli(1,2*pi*r/L)*cos(2*pi*z/L)))^2/(2*me*H0);
+ rbplusphi(i)=fzero(plus,rbplusphi(i-1));
+ remainderphi(i)=plus(rbplusphi(i));
+ catch
+ end
+end
+
+f=figure;
+%subplot(2,1,1)
+plot(phia(remainderphi~=1)/1e3,rb_phi(remainderphi~=1)*1e3,'x-','displayname','r_b^-')
+hold on
+plot(phia(remainderphi~=1)/1e3,rbplusphi(remainderphi~=1)*1e3,'x-','displayname','r_b^+')
+plot(phiasim/1e3,rplussim*1e3,'+--','displayname','r_b^+ simulated','linewidth',1.5)
+plot(phiasim/1e3,r_sim*1e3,'+--','displayname','r_b^- simulated','linewidth',1.5)
+legend
+ylabel('r_b [mm]')
+ylim([0 160])
+xlabel('\Phi_a [kV]')
+% subplot(2,1,2)
+% semilogx(ne(remainder~=1),remainder(remainder~=1),'x--')
+% ylabel('remainder [a.u.]')
+% xlabel('n_e(r_b^-) [m^{-3}]')
+
+f.PaperOrientation='landscape';
+f.PaperUnits='centimeters';
+papsize=[10 10];
+f.PaperSize=papsize;
+print(f,'rbne','-dpdf','-fillpage')
+
+
+
+i=find(consideredphia>=phia,1,'first');
+z=linspace(0,L/2,6000);
+rb_z=zeros(size(z));
+rbplusz=rb_z;
+
+minus=@(r) 1+qe*((phib-phia(i))*log(r)+phia(i)*log(b)-phib*log(a))/(log(b/a))/H0-(P0/r+qe*0.5*B0*(r-L/pi*(Rcurv-1)/(Rcurv+1)*besseli(1,2*pi*r/L)*cos(2*pi*z(1)/L)))^2/(2*me*H0);
+rb_z(1)=fzero(minus,rb_phi(i));
+rb_t=rb_z(1);
+z_t=z(1);
+plus=@(r) 1+qe*((phib-phia(i))*log(r)+phia(i)*log(b)-phib*log(a))/(log(b/a))/H0-(P0/r+qe*0.5*B0*(r-L/pi*(Rcurv-1)/(Rcurv+1)*besseli(1,2*pi*r/L)*cos(2*pi*z(1)/L)))^2/(2*me*H0);
+rbplusz(1)=fzero(plus,rbplusphi(i));
+remainder(1)=plus(rbplusz(1));
+for j=2:length(z)
+minus=@(r) 1+qe*((phib-phia(i))*log(r)+phia(i)*log(b)-phib*log(a))/(log(b/a))/H0-(P0/r+qe*0.5*B0*(r-L/pi*(Rcurv-1)/(Rcurv+1)*besseli(1,2*pi*r/L)*cos(2*pi*z(j)/L)))^2/(2*me*H0);
+format long;
+try
+rb_z(j)=fzero(minus,rb_z(j-1));
+z_t=z(j);
+plus=@(r) 1+qe*((phib-phia(i))*log(r)+phia(i)*log(b)-phib*log(a))/(log(b/a))/H0-(P0/r+qe*0.5*B0*(r-L/pi*(Rcurv-1)/(Rcurv+1)*besseli(1,2*pi*r/L)*cos(2*pi*z_t/L)))^2/(2*me*H0);
+rbplusz(j)=fzero(plus,rbplusz(j-1));
+remainder(j)=plus(rbplusz(j));
+catch
+end
+end
+figure()
+z=[-flip(z(2:end)) z];
+rb_z=[flip(rb_z(2:end)) rb_z];
+rbplusz=[flip(rbplusz(2:end)) rbplusz];
+plot(z(rb_z~=0),rb_z(rb_z~=0),'b')
+hold on
+plot(z(rb_z~=0),rbplusz(rb_z~=0),'r')
+ylabel('r_b [m]')
+ylim([0 0.16])
+xlim([-L/2 L/2])
+xlabel('z [m]')
+title(sprintf('\\phi_a=%5.2f kV',phia(i)/1e3))
+
+ind=rbplusz~=0 & rb_z~=0& ~isnan(rbplusz) & ~isnan(rb_z);
+ra=rb_z(ind);
+rb=rbplusz(ind);
+zvol=z(ind);
+ra=(ra(1:end-1)+ra(2:end))/2;
+rb=(rb(1:end-1)+rb(2:end))/2;
+VOL=2*pi*diff(zvol).*min(ra).*(rb-ra);
+
+nplasma=2116800
+Nbesl=ne*VOL;
+Nbe=sum(Nbesl)
+Nmacrosl=floor(Nbesl/Nbe*nplasma);
+Nmacro=sum(Nmacrosl)
+remain=nplasma-Nmacro
+Nmacrosl(ceil(length(Nmacrosl)/2-remain/2):floor(length(Nmacrosl)/2+remain/2))=Nmacrosl(ceil(length(Nmacrosl)/2-remain/2):floor(length(Nmacrosl)/2+remain/2))+1;
+Nmacrosl(floor(length(Nmacrosl)/2))=Nmacrosl(floor(length(Nmacrosl)/2))+nplasma-sum(Nmacrosl);
+Nmacro=sum(Nmacrosl)
+remain=nplasma-Nmacro
+
+
+msim=me*Nbe/nplasma
+qsim=-qe*Nbe/nplasma
+
+datas.ra=ra;
+datas.rb=rb;
+datas.z=zvol;
+datas.npartsslice=Nmacrosl;
+datas.m=me;
+datas.weight=Nbe/nplasma;
+datas.q=-qe;
+datas.H0=H0;
+datas.P0=P0;
+datas.temperature=10000;
+datas.velocitytype=2;
+datas.radialtype=1;
\ No newline at end of file
diff --git a/matlab/timescales.m b/matlab/timescales.m
new file mode 100644
index 0000000..d64fe6f
--- /dev/null
+++ b/matlab/timescales.m
@@ -0,0 +1,59 @@
+Pn=1e-8;%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;
+
+
+Te=1; %eV
+ve=sqrt(Te*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,Te);
+sigion=sigmabeb(41.72,71.13,2,1,Te)+sigion;
+sigion=sigmabeb(21,63.18,2,1,Te)+sigion;
+sigion=sigmabeb(17.07,44.30,4,1,Te)+sigion;
+
+Coulomblog=log(sqrt(eps_0*Te/qe*ne)/qe^2*2*pi*eps_0*me*ve^2);
+
+tauion=1./(nb*sigion*ve);
+tauloss=tauion*log(ne);
+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));
+
+figure
+loglog(ne,tauion*ones(size(ne)),'displayname','\tau_{ion}')
+hold on
+loglog(ne,tauloss,'displayname','\tau_{accum}')
+loglog(ne,tauthermal,'displayname','\tau_{thermal}')
+loglog(ne,taucycl*ones(size(ne)),'displayname','\tau_{ce}')
+loglog(ne,taupe,'displayname','\tau_{pe}')
+legend('location','southeast')
+title(sprintf('E_{in}=%.2f eV, n_b=%3.2e [m^{-3}]',Te,nb))
+xlabel('n_e [m^{-3}]')
+ylabel('\tau [s]')
+
+% 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