diff --git a/.gitignore b/.gitignore
index 2b10b2a..644c56a 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,13 +1,14 @@
 *.h5
 *.h5_old
 *_genmod.f90
 *.mod
 *.o
 *.out
 *.m~
 *.trace
 *.log
 dataanalysis/**/*
 src/release
 src/debug
 wk/advi
+f2matlab
diff --git a/matlab/CheckAxialConfinement.m b/matlab/CheckAxialConfinement.m
index 04f864d..3168910 100644
--- a/matlab/CheckAxialConfinement.m
+++ b/matlab/CheckAxialConfinement.m
@@ -1,61 +1,62 @@
 t2dlength=size(M.t2d,1);
 tstep=t2dlength;
-rpos=160;
-zpos=355;
+rpos=96;
+zpos=320;
 rAthet_pos=M.rAthet(rpos,zpos);
 vpar0=sqrt(M.kb*22000/M.msim*M.weight);
 
 
 model=M.potentialwellmodel(tstep);
 N=M.N(:,:,tstep);
 Psieval=rAthet_pos*ones(M.nz+1,1);
 Zeval=M.zgrid;
 
 %rdisp=sort(unique([M.rAthet(:,end);M.rAthet(end,:)]));
 [Zmesh,Rmesh]=meshgrid(M.zgrid,M.rgrid);
 pot=griddata(model.z,model.rathet,model.pot,Zeval,Psieval);
 [Zinit,~]=meshgrid(M.zgrid,M.rAthet(:,1));
 N0=griddata(Zinit,M.rAthet,N,Zeval,Psieval);
 pot=pot+min(pot);
 
 % Magnetic field mirror ratio at each grid position
 % compared to local position
 R=griddata(Zmesh,M.rAthet,M.B',Zeval,Psieval,'natural')/M.B(zpos,rpos);
 
 % Electrostatic potential on magnetic field line
 % coordinates
 phis=griddata(Zmesh,M.rAthet,M.pot(:,:,tstep),Zeval,Psieval,'natural')-M.pot(rpos,zpos,tstep);
 
 rposline=griddata(Zmesh,M.rAthet,Rmesh,Zeval,Psieval,'natural');
 
 %%
 figure
 subplot(3,1,1)
 contourf(M.zgrid,M.rgrid,N,'edgecolor','none')
 hold on
 plot(M.zgrid(zpos),M.rgrid(rpos),'rx')
 plot(Zeval,rposline,'r--')
 contour(M.zgrid,M.rgrid,M.geomweight(:,:,1),[0 0],'w-.','linewidth',1.5,'Displayname','Vessel Boundaries');
+ylim([1e-3 12e-3])
 
 subplot(3,1,2)
 plot(M.zgrid,pot)
 ylabel('-e\Delta\phi [eV]')
 yyaxis right
 plot(M.zgrid,N0)
 ylabel('n [m^{-3}]')
 xlabel('z [m]')
 
 subplot(3,1,3)
 EparB=0.5*M.msim/M.weight*((M.Ez(rpos,zpos,tstep)*M.Br(zpos,rpos)-M.Er(rpos,zpos,tstep)*M.Bz(zpos,rpos))/M.B(zpos,rpos)^2)^2*(1-R)/M.qe;
 Eparphi=-M.qsim/M.weight*phis/M.qe;
 Epar=EparB+Eparphi+0.5*M.msim/M.weight/M.qe*vpar0^2;
 plot(M.zgrid,Epar,'displayname','E_{par,tot}')
 hold on
 plot(M.zgrid,EparB,'displayname','E_{par,mag}')
 plot(M.zgrid,Eparphi,'displayname','E_{par,elec}')
 plot(M.zgrid,0.5*M.msim/M.weight*(vpar0^2)/M.qe*ones(size(M.zgrid)),'displayname','E_{par,0}')
 grid on
 xlabel('z [m]')
 ylabel('E_{par} [eV]')
 legend('location','south')
-ylim(1000*[-1 1])
\ No newline at end of file
+%ylim(1000*[-1 1])
\ No newline at end of file
diff --git a/matlab/CheckDiocStability_post.m b/matlab/CheckDiocStability_post.m
index cd0c731..6d524f7 100644
--- a/matlab/CheckDiocStability_post.m
+++ b/matlab/CheckDiocStability_post.m
@@ -1,126 +1,126 @@
 nbzid=15;
 zindices=linspace(length(M.zgrid)/6,length(M.zgrid)*5/6,nbzid);
 zconsid=1:nbzid;
 
 for zid=zconsid%1:length(zindices)%;%floor(length(M.zgrid)/2);
     zindex=floor(zindices(zid));
     zfolder=sprintf('diocz_%03d',zindex);
     if ~isfolder(zfolder) || ~isfile(sprintf('%s/results.mat',zfolder))
         fprintf('Warning: result %s doesn''t exist\n',zfolder)
         continue
     end
     result=load(sprintf('%s/results.mat',zfolder));
     if ~isfield(result,'geomweight')
         result.geomweight=M.geomweight(:,:,1);
     end
     rindex=M.geomweight(:,zindex,1)>=0;
     
     for l=result.lsteps
         
         r=result.r(:,l);
         phin=result.phinum(:,:,l);
         w=result.w(:,l);
         n=result.n;
         omegar=result.omegadrift;
         rgrid=result.rgrid;
         omegar=omegar(rindex,zindex);
         dn=(n(3:end)-n(1:end-2))./(rgrid(3:end)-rgrid(1:end-2));
         a=min(rgrid);
         b=max(rgrid);
         idmin=find(r==a);
         idmax=find(r==b);
         wd=(M.qsim/M.weight)^2*n/M.eps_0/(M.msim/M.weight)/2/result.wce;
         f=figure;
         subplot(4,1,1)
         hold on
         xlabel('R [m]')
         ylabel('n [m^{-3}]')
         plots=gobjects(2,1);
         plots(1)=plot(rgrid,n,'displayname','n');
         yyaxis right
         plots(2)=plot(rgrid,omegar,'displayname','\omega_{re}');
         hold on
         %plots(3)=plot(rgrid,wd,'displayname','\omega_{d}');
         legend(plots,{'n','\omega_r','\omega_d'},'location','northwest')
         ylabel('\omega_{re} [rad\cdots^{-1}]')
         grid minor
         
         ax=gca;
         ax.LineStyleOrder={'-','--','-.',':'};
         subplot(4,1,[2,3])
         hold on
         p=gobjects();
         j=1;
         for i=1:size(phin,1)
             if imag(w(i))/abs(real(w(i)))>1e-4
                 freq=w(i);
                 disphinum=phin(i,:);
                 disphinum=disphinum/max(abs(disphinum));
                 yyaxis left
                 p(j)=plot(r(idmin:idmax),[0 real(disphinum(idmin:idmax-2)) 0],'displayname',sprintf('\\omega=%1.3g%+1.3gi',real(freq),imag(freq)));
                 j=j+1;
                 yyaxis right
                 plot(r(idmin:idmax),[0 imag(disphinum(idmin:idmax-2)) 0],'displayname',sprintf('\\omega=%1.3g%+1.3gi',real(freq),imag(freq)))
             end
         end
         subplot(4,1,1)
         if (abs(imag(w(1))/real(w(1)))>1e-10)
-            title(sprintf('Unstable case l=%d, \\omega_d=%1.3g, z=%1.3f [mm]',l,max(wpe2/(2*wce)),M.zgrid(zindex)*1e3))
+            title(sprintf('Unstable case l=%d, \\omega_d=%1.3g, z=%1.3f [mm]',l,max(wd),M.zgrid(zindex)*1e3))
         else
-            title(sprintf('Stable case l=%d, \\omega_d=%1.3g z=%1.3f [mm]',l,max(wpe2/(2*wce)),M.zgrid(zindex)*1e3))
+            title(sprintf('Stable case l=%d, \\omega_d=%1.3g z=%1.3f [mm]',l,max(wd),M.zgrid(zindex)*1e3))
         end
         subplot(4,1,[2,3])
         yyaxis left
         ylim([-1 1])
         xlabel('R [m]')
         ylabel('real(\delta\phi) [a.u.]')
         grid on
         yyaxis right
         ylim([-1 1])
         grid minor
         ylabel('imag(\delta\phi) [a.u.]')
         f.PaperUnits='centimeters';
         f.PaperSize=[16,20];
         legend('location','northwest')
         if length(p)>=1 && all(ishghandle(p))
             legend(p)
         end
         subplot(4,1,4)
         plot(real(w),imag(w),'x')
         xlabel('\omega_r [rad\cdots^{-1}]')
         ylabel('\omega_i [rad\cdots^{-1}]')
         subplot(4,1,1)
         ax1=gca;
         subplot(4,1,[2,3])
         ax2=gca;
         drawnow
         posit1=get(ax1,'position');
         posit2=get(ax2,'position');
         posit2(3)=posit1(3);
         set(ax2,'position',posit2);
         print(f,sprintf('%s/%s_dioc_l%d',zfolder,M.name,l),'-dpdf','-fillpage')
         savefig(f,sprintf('%s/%s_dioc_l%d',zfolder,M.name,l))
         close(f) 
         f=figure;
         subplot(2,1,1)
         hold on
         plot(rgrid,n/max(n),'displayname','n')
         plot(rgrid(2:end-1),dn/max(abs(dn)),'displayname','dn')
         plot(rgrid,(real(w(1))-l*omegar)/max(abs((real(w(1))-l*omegar))),'displayname','w-lwre')
         legend
         grid on
         xlabel('r [m]')
         ylabel('Normalized quantity a.u.')
         subplot(2,1,2)
         disphinum=phin(1,:);
         disphinum=disphinum/max(abs(disphinum));
         p(j)=plot(r(idmin:idmax),[0 real(disphinum(idmin:idmax-2)) 0],'displayname',sprintf('\\omega=%1.3g%+1.3gi',real(freq),imag(freq)));
         xlabel('r [m]')
         ylabel('Re(\delta\phi) a.u.')
         f.PaperUnits='centimeters';
         f.PaperSize=[12,6];
         print(f,sprintf('%s/%s_maxw_l%d',zfolder,M.name,l),'-dpdf','-fillpage')
         savefig(f,sprintf('%s/%s_maxw_l%d',zfolder,M.name,l))
         close(f) 
     end
 end
\ No newline at end of file
diff --git a/matlab/CheckDiocStability_post_brd.m b/matlab/CheckDiocStability_post_brd.m
index fa9e474..67a9570 100644
--- a/matlab/CheckDiocStability_post_brd.m
+++ b/matlab/CheckDiocStability_post_brd.m
@@ -1,92 +1,92 @@
 nbzid=15;
 zindices=linspace(length(M.zgrid)/6,length(M.zgrid)*5/6,nbzid);
-zconsid=2:nbzid;
+zconsid=2:2:nbzid-1;
 
+
+dens=mean(M.N(:,:,fieldstart:fieldend),3);
+wd=M.qe*dens./(2*M.eps_0*M.B');
+
+wd=max(wd(:));
 f=figure('Name', sprintf('%s dioc_stab',M.name));
 fieldstart=max(1,length(M.t2d)-500);
 fieldend=length(M.t2d);
-ax1=subplot(3,1,1);
+ax1=subplot(2,1,1);
 geomw=M.geomweight(:,:,1);
 dens=mean(M.N(:,:,fieldstart:fieldend),3);
 dens(geomw<=0)=0;
-[C,h]=contourf(ax1,M.zgrid,M.rgrid,dens,'LineColor','none');
+[C,h]=contourf(ax1,M.zgrid*1e3,M.rgrid*1e3,dens,'LineColor','none');
 set(h,'LineColor','none')
 hold on
-contour(ax1,M.zgrid,M.rgrid,M.geomweight(:,:,1),[0 0],'w-.','linewidth',1.5,'Displayname','Boundaries');
+contour(ax1,M.zgrid*1e3,M.rgrid*1e3,M.geomweight(:,:,1),[0 0],'w-.','linewidth',1.5,'Displayname','Boundaries');
 if(M.conformgeom)
-    ylim(ax1,[M.rgrid(1) M.rgrid(sum(M.nnr(1:2)))])
+    ylim(ax1,[M.rgrid(1) M.rgrid(sum(M.nnr(1:2)))]*1e3)
 else
-    ylim(ax1,[M.rgrid(1) M.rgrid(end)])
+    ylim(ax1,[M.rgrid(1) M.rgrid(end)]*1e3)
 end
-xlim(ax1,[M.zgrid(1) M.zgrid(end)])
-xlabel(ax1,'z [m]')
-ylabel(ax1,'r [m]')
+xlim(ax1,[M.zgrid(1) M.zgrid(end)]*1e3)
+xlabel(ax1,'z [mm]')
+ylabel(ax1,'r [mm]')
 c = colorbar(ax1);
 c.Label.String= 'n[m^{-3}]';
 view(ax1,2)
 grid on;
 hold on;
 
 zpos=zeros(size(zconsid));
 wrslt=zeros(size(zconsid));
 lmax=zeros(size(zconsid));
 for zid=1:length(zconsid)%1:length(zindices)%;%floor(length(M.zgrid)/2);
     zindex=floor(zindices(zconsid(zid)));
-    plot(ax1,M.zgrid(zindex)*[1 1],[M.rgrid(1) M.rgrid(end)],'r--')
-    zpos(zid)=M.zgrid(zindex);
+    plot(ax1,M.zgrid(zindex)*[1 1]*1e3,[M.rgrid(1) M.rgrid(end)]*1e3,'r--')
+    zpos(zid)=M.zgrid(zindex)*1e3;
     zfolder=sprintf('diocz_%03d',zindex);
     if ~isfolder(zfolder) || ~isfile(sprintf('%s/results.mat',zfolder))
         fprintf('Warning: result %s doesn''t exist\n',zfolder)
         continue
     end
     result=load(sprintf('%s/results.mat',zfolder));
     if ~isfield(result,'geomweight')
         result.geomweight=M.geomweight(:,:,1);
     end
     rindex=M.geomweight(:,zindex,1)>=0;
     
     for l=result.lsteps
         
         r=result.r(:,l);
         phin=result.phinum(:,:,l);
         w=result.w(:,l);
         [~,idwmax]=max(imag(w));
         if(imag(w(idwmax))>imag(wrslt(zid)))
             wrslt(zid)=w(idwmax);
             lmax(zid)=l;
         end
     end
 end
-ax2=subplot(3,1,2);
-plot(zpos,real(wrslt),'x')
-xlabel('z [m]')
-xlim([M.zgrid(1) M.zgrid(end)])
-ylabel('real(\omega)')
+ax2=subplot(2,1,2);
+semilogy(zpos,real(wrslt)/wd,'bx','displayname','Real')
+xlabel('z [mm]')
+hold on
+semilogy(zpos,imag(wrslt)/wd,'b+','displayname','Imag')
+xlim([M.zgrid(1) M.zgrid(end)]*1e3)
+ylabel('\omega/\omega_d')
+ylim([1e-2 2])
 yyaxis right
-plot(zpos,imag(wrslt),'+')
-ylabel('imag(\omega)')
+plot(zpos,lmax,'d','displayname','l')
+ylabel('most unstable l')
 grid on
+title('Most unstable modes')
+legend('location','east')
 
-ax3=subplot(3,1,3);
-plot(zpos,lmax,'x')
-xlabel('z [m]')
-xlim([M.zgrid(1) M.zgrid(end)])
-ylabel('Most unstable l')
-ylim([0 max(result.lsteps)])
 
-grid on
 drawnow
 posit1=get(ax1,'position');
 posit2=get(ax2,'position');
 posit2(3)=posit1(3);
 set(ax2,'position',posit2);
-posit2=get(ax3,'position');
-posit2(3)=posit1(3);
-set(ax3,'position',posit2);
 
 f.PaperOrientation='landscape';
 f.PaperUnits='centimeters';
-papsize=[16 14];
+papsize=[10 12];
 f.PaperSize=papsize;
 print(f,sprintf('%sdioc_stab',M.name),'-dpdf','-fillpage')
 savefig(f,sprintf('%sdioc_stab',M.name))
\ No newline at end of file
diff --git a/matlab/analyse_espicfields.m b/matlab/analyse_espicfields.m
index 03a8eed..9831748 100644
--- a/matlab/analyse_espicfields.m
+++ b/matlab/analyse_espicfields.m
@@ -1,453 +1,391 @@
 %% time extent of plotted variables
 % 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-1000);%floor(0.95*t2dlength);
+fieldstart=max(1,t2dlength-500);%floor(0.95*t2dlength);
 fieldend=t2dlength;%23;
 deltat=t2dlength-fieldstart;
 %[~, rgridend]=min(abs(M.rgrid-0.045));
 rgridend=sum(M.nnr(1:2));
 [~, name, ~] = fileparts(M.file);
 
 M.displayenergy
 
 
 
 %% Radial profile 
 try
     M.displayrprofile([1,fieldstart:fieldend],[],true)
 catch
 end
 
 %%
 %fieldstart=2300; fieldend=2500;
 dens=mean(M.N(:,:,fieldstart:fieldend),3);
-<<<<<<< HEAD
 geomw=M.geomweight(:,:,1);
 geomw(geomw<0)=0;
 geomw(geomw>0)=max(max(dens));
 dispdens=dens;
-dispdens(geomw==0)=NaN;
+dispdens(geomw<=0)=NaN;
 pot=mean(M.pot(:,:,fieldstart:fieldend),3);
-brillratio=2*dens*M.me./(M.eps_0*M.Bz'.^2);
+brillratio=2*dens*M.me./(M.eps_0*M.B'.^2);
 [R,Z]=meshgrid(M.rgrid,M.zgrid);
 Rinv=1./R;
 Rinv(isnan(Rinv))=0;
 VTHET=mean(M.fluidUTHET(:,:,fieldstart:fieldend),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')
 
 %% Density
 M.displaydensity(fieldstart,fieldend);
 
 %% Fieldswide
 f=figure('Name', sprintf('%s fields',name));
 ax1=gca;
 h=contourf(ax1,M.zgrid,M.rgrid,dispdens,'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',real(levels),'r-.','linewidth',1.5,'Displayname','Magnetic field lines');
 [~, hContour]=contourf(ax1,M.zgrid,M.rgrid,geomw,2);
 drawnow;
 xlim(ax1,[M.zgrid(1) M.zgrid(end)])
 ylim(ax1,[M.rgrid(1) M.rgrid(end)])
 legend({'n_e [m^{-3}]','Magnetic field lines'},'location','southwest')
 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(fieldend),double(maxdens)))
 c = colorbar(ax1);
 c.Label.String= 'n[m^{-3}]';
 view(ax1,2)
 %set(h,'edgecolor','none');
 grid on;
 hFills=hContour.FacePrims;
 [hFills.ColorType] = deal('truecoloralpha');  % default = 'truecolor'
 hFills(1).ColorData = uint8([255;255;255;255]);
 for idx = 2 : numel(hFills)
    hFills(idx).ColorData(4) = 0;   % default=255
 end
 %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*1000,M.rgrid*1000,dispdens,'Displayname','n_e [m^{-3}]');
 hold on;
 pot(M.geomweight(:,:,1)<0)=0;
 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);
+%levels=logspace( log10(min(min(Blines(:,2:end)))), log10(max(max(Blines(:,:)))),10);
+levels=linspace(min(min(Blines(:,2:end))),max(max(Blines(:,:))),10);
 [~,h1]=contour(ax1,M.zgrid*1000,M.rgrid*1000,Blines',real(levels),'r-.','linewidth',1.5,'Displayname','Magnetic field lines');
-
+max(abs(pot(:)))/1e3
 [c1,h2]=contour(ax1,M.zgrid*1000,M.rgrid*1000,pot./1e3,'m--','ShowText','on','linewidth',1.2,'Displayname','Equipotentials [kV]');
 clabel(c1,h2,'Color','white')
 
 %contour(ax1,M.zgrid*1000,M.rgrid*1000,M.geomweight(:,:,1),[0 0],'w-','linewidth',1.5);
 
 [c1,hContour]=contourf(ax1,M.zgrid*1000,M.rgrid*1000,geomw);
 drawnow;
 xlim(ax1,[M.zgrid(1)*1000 M.zgrid(end)*1000])
 if(M.conformgeom)
     ylim(ax1,[M.rgrid(1)*1000 M.rgrid(rgridend)*1000])
 else
     ylim(ax1,[M.rgrid(1)*1000 M.rgrid(end)*1000])
 end
 legend([h1,h2],{'Magnetic field lines','Equipotentials [kV]'},'location','southwest')
 xlabel(ax1,'z [mm]')
 ylabel(ax1,'r [mm]')
 %title(ax1,sprintf('Density t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(fieldstart),M.t2d(fieldend),double(maxdens)))
 c = colorbar(ax1);
 c.Label.String= 'n[m^{-3}]';
 view(ax1,2)
 %set(h,'edgecolor','none');
 grid on;
 hFills=hContour.FacePrims;
 [hFills.ColorType] = deal('truecoloralpha');  % default = 'truecolor'
 try
 hFills(1).ColorData = uint8([150;150;150;255]);
 for idx = 2 : numel(hFills)
    hFills(idx).ColorData(4) = 0;   % default=255
 end
 catch
 end
 f.PaperOrientation='landscape';
 [~, name, ~] = fileparts(M.file);
+if M.maxwellsrce.present
+    rlen=diff(M.maxwellsrce.rlim);
+    zlen=diff(M.maxwellsrce.zlim);
+rectangle('Position',[M.maxwellsrce.zlim(1) M.maxwellsrce.rlim(1) zlen rlen]*1000,'Edgecolor','g','Linewidth',3,'Linestyle','--')
+end
+%rectangle('Position',[-10 73 44 8],'Edgecolor','g','Linewidth',3,'Linestyle','--')
 f.PaperUnits='centimeters';
-papsize=[20 16];
+papsize=[12 8];
 f.PaperSize=papsize;
 rectangle('Position',[M.zgrid(1) M.rgrid(1) M.zgrid(end)-M.zgrid(1) 0.064-M.rgrid(1)]*1000,'FaceColor',[150 150 150]/255)
 print(f,sprintf('%sFields',name),'-dpdf','-fillpage')
 savefig(f,sprintf('%sFields',name))
 
 %% Brillouin
 f=figure('Name', sprintf('%s Brillouin',name));
 ax1=gca;
+brillratio(geomw<=0)=NaN;
 h=surface(ax1,M.zgrid,M.rgrid,brillratio);
 set(h,'edgecolor','none');
 xlim(ax1,[M.zgrid(1) M.zgrid(end)])
 if(M.conformgeom)
     ylim(ax1,[M.rgrid(1) M.rgrid(rgridend)])
 else
     ylim(ax1,[M.rgrid(1) M.rgrid(end)])
 end
 xlabel(ax1,'z [m]')
 ylabel(ax1,'r [m]')
 title(ax1,'Position')
 c = colorbar(ax1);
 c.Label.String= 'Brillouin Ratio';
 view(ax1,2)
 caxis([0 max(1.2,max(brillratio(:)))])
 title(ax1,sprintf('Brillouin Ratio t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(fieldstart),M.t2d(fieldend),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;
+omegare(geomw<=0)=NaN;
 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)])
 if(M.conformgeom)
     ylim(ax3,[M.rgrid(1) M.rgrid(rgridend)])
 else
     ylim(ax3,[M.rgrid(1) M.rgrid(end)])
 end
 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(fieldend),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))
 
 %% VTHET
 f=figure('Name', sprintf('%s Vthet',name));
 ax3=gca;
+VTHET(geomw<=0)=NaN;
 surface(ax3,M.zgrid,M.rgrid,VTHET,'edgecolor','None')
 xlabel(ax3,'z [m]')
 ylabel(ax3,'r [m]')
 xlim(ax3,[M.zgrid(1) M.zgrid(end)])
 if(M.conformgeom)
     ylim(ax3,[M.rgrid(1) M.rgrid(rgridend)])
 else
     ylim(ax3,[M.rgrid(1) M.rgrid(end)])
 end
 colormap(ax3,'jet')
 c = colorbar(ax3);
 c.Label.String= '|V_\theta| [m/s]';
 %c.Limits=[min(M.fluidUTHET(:)) max(M.fluidUTHET(:))];
 %caxis(ax3,[min(M.fluidUTHET(:)) max(M.fluidUTHET(:))])
 title(ax3,sprintf('Azimuthal velocity t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(fieldstart),M.t2d(fieldend),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_vthet',name),'-dpdf','-fillpage')
 savefig(f,sprintf('%sfluid_vthet',name))
 
 %% vdrift
 f=figure('Name', sprintf('%s V drift',name));
 ax3=gca;
 vdrift=(mean(M.Ez(:,:,fieldstart:end),3).*M.Br'-mean(M.Er(:,:,fieldstart:end),3).*M.Bz')./(M.B.^2)';
+vdrift(geomw<=0)=NaN;
 v=vdrift;
 surface(ax3,M.zgrid,M.rgrid,v,'edgecolor','None')
 xlabel(ax3,'z [m]')
 ylabel(ax3,'r [m]')
 xlim(ax3,[M.zgrid(1) M.zgrid(end)])
 if(M.conformgeom)
     ylim(ax3,[M.rgrid(1) M.rgrid(rgridend)])
 else
     ylim(ax3,[M.rgrid(1) M.rgrid(end)])
 end
 colormap(ax3,'jet')
 c = colorbar(ax3);
 set(ax3,'zscale','log')
 set(ax3,'colorscale','log')
 c.Label.String= 'V_{ExB}';
 %c.Limits=[min(M.fluidUTHET(:)) max(M.fluidUTHET(:))];
 %caxis(ax3,[min(M.fluidUTHET(:)) max(M.fluidUTHET(:))])
 title(ax3,sprintf('V_{ExB} velocity t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(fieldstart),M.t2d(fieldend),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_EB',name),'-dpdf','-fillpage')
 savefig(f,sprintf('%sfluid_EB',name))
 
 
 %% HP
 if(M.R.nt>=2)
 taveragestart=max(floor(0.8*size(M.tpart,1)),2);
 M.displayHP(taveragestart);
 end
 
 %% VR, VTHET, VZ
-if(M.R.nt>=2)
-timesteppart=length(M.tpart);
-
-fig=figure;
-ax1=gca;
-nbp=min(M.nbparts(1),M.R.nparts);
-Vr=M.VR(1:nbp,1,false)/M.vlight;
-Vz=M.VZ(1:nbp,1,false)/M.vlight;
-Vthet=M.VTHET(1:nbp,1,false)/M.vlight;
-Vperp=sqrt(Vthet.^2+Vr.^2);
-
-nbp=min(M.nbparts(end),M.R.nparts);
-Vrend=M.VR(1:nbp,timesteppart,false)/M.vlight;
-Vzend=M.VZ(1:nbp,timesteppart,false)/M.vlight;
-Vthetend=M.VTHET(1:nbp,timesteppart,false)/M.vlight;
-Vperpend=sqrt(Vthetend.^2+Vrend.^2);
-
-binwidth=abs(max(Vrend)-min(Vrend))/sqrt(length(Vrend));
-f=figure('Name',sprintf("%s v distrib",M.file));
-tstudied=0;
-legtext=sprintf("t=%2.1f [ns]",M.tpart(end)*1e9);
-ax1=subplot(1,3,1);
-if(length(Vr)>1)
-h1=histogram(Vr,'Binwidth',binwidth,'Normalization','count','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*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));
-ylabel('counts')
-xlabel('\beta_r')
-
-grid on
-
-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));
-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));
-ylabel('counts')
-xlabel('\beta_\theta')
-legend(ax2,'location','northoutside','orientation','vertical')
-grid on
-
-ax3=subplot(1,3,3);
-%h1=histogram(Vz(:,1),'Binwidth',binwidth,'Normalization','probability','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9));
-if(length(Vz)>1)
-h1=histogram(ax3,Vz(:,1));
-set(h1,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9));
-end
-hold on
-%h1=histogram(Vzend,'Binwidth',binwidth,'Normalization','probability','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9));
-h1=histogram(ax3,Vzend);
-set(h1,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9));
-ylabel('counts')
-xlabel('\beta_z')
-grid on
-f=gcf;
-sgtitle(sprintf('dt=%1.2e[ns]',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))
+M.displayVdistribRThetZ()
 
-end
+%% Vperp Vpar
+M.displayVdistribParPer()
 
 
 %% 0D diagnostics
 BrillouinRatioth=2*M.omepe^2/M.omece^2
 BrillouinMax=max(max(brillratio))
 NpartsTrapped=mean(M.npart(end-10:end))*M.msim/M.me
 dblengthiter=fieldstart:fieldend;
 
 
 Pr=squeeze(M.Presstens(1,:,:,dblengthiter));
 Pz=squeeze(M.Presstens(6,:,:,dblengthiter));
 enddens=M.N(:,:,dblengthiter);
 meanPr=mean(Pr,3);
 meanPz=mean(Pz,3);
 Tr=Pr./enddens;
 Tz=Pz./enddens;
 Debye_length=sqrt(abs(min(min(min(mean(Tr(Tr>0)./enddens(Tr>0),3))),min(min(mean(Tz(Tz>0)./enddens(Tz>0),3)))))*M.eps_0/M.qe^2)
 
 
 %% Debye length
 f=figure('Name', sprintf('%s Debye',name));
 subplot(2,1,1)
 surface(M.zgrid,M.rgrid,sqrt(mean(abs(Tr./enddens*M.eps_0/M.qe^2),3)),'edgecolor','none')
 colorbar
 f.PaperOrientation='landscape';
 xlabel('z [m]')
 ylabel('r [m]')
 c = colorbar;
 c.Label.String= '\lambda_r [m]';
 
 subplot(2,1,2)
 surface(M.zgrid,M.rgrid,sqrt(mean(abs(Tz./enddens*M.eps_0/M.qe^2),3)),'edgecolor','none')
 colorbar
 f.PaperOrientation='landscape';
 xlabel('z [m]')
 ylabel('r [m]')
 c = colorbar;
 c.Label.String= '\lambda_z [m]';
 
 [~, name, ~] = fileparts(M.file);
 f.PaperUnits='centimeters';
 papsize=[16 14];
 f.PaperSize=papsize;
 print(f,sprintf('%s_dblength',name),'-dpdf','-fillpage')
 savefig(f,sprintf('%s_dblength',name))
 
 %% Temperature
 f=figure('Name', sprintf('%s Temperature',name));
 ax1=subplot(2,1,1);
 title('Radial temperature')
 surface(M.zgrid,M.rgrid,mean(Tr,3)/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,mean(Tz,3)/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))
 
 
 %% Wells
 M.display2Dpotentialwell(t2dlength-1);
 M.display2Dpotentialwell(t2dlength-1,false);
 
 
 
diff --git a/matlab/coaxial_ion_loss_criteria.m b/matlab/coaxial_ion_loss_criteria.m
index 8d0e55d..57b4700 100644
--- a/matlab/coaxial_ion_loss_criteria.m
+++ b/matlab/coaxial_ion_loss_criteria.m
@@ -1,287 +1,287 @@
 %% 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=6e17;
 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.0079;
 zmin=0.0;
 zmax=-width/2;
 deltar=(R-1)/(R+1)*besseli(1,2*pi*r/width)*(width/2/pi)*(cos(2*pi*zmin/width)-cos(2*pi*zmax/width));
 rb_=0.00764;
 rbplus=0.0842;
 Te=400; %eV
 %1rbplus=0.020;
 
 % Davidson test
 prefix='Davidson';
 phia=-60000;
 phib=0;
 R=1.5;
 B0=0.21;
 width=0.48;
 b=0.16;
 a=0.001;
 rb_=0.0076;
 r=9.7e-3 ;
 zmin=0.0;
 zmax=-width/8;
 deltar=(R-1)/(R+1)*besseli(1,2*pi*r/width)*(width/2/pi)*(cos(2*pi*zmin/width)-cos(2*pi*zmax/width));
 rbplus=rb_ + max(deltar,0.0026);
 rcenter=r+deltar
 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;
+%ITER
+prefix='ITER' 
+R=1.0001;
+phia=-5000;
+phib=0;
+b=0.0792;
+a=0.064;
+r=0.080;
+deltar=(R-1)/(R+1)*r;
+rb_=0.0766;
+rbplus=0.079;
 
 spcharge=@(ne)-qe*ne/2/eps_0;
 
 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((r+deltar)/r);
 
 
 figure
-rforphi=linspace(0.001,0.16,1000);
+rforphi=linspace(a,b,1000);
 phir=zeros(length(rforphi),1);
 for i=1:length(rforphi)
-    phir(i)=Phi(rforphi(i),rb_,rbplus,a,b,phia,phib, 3e16);
+    phir(i)=Phi(rforphi(i),rb_,rbplus,a,b,phia,phib, 6e17);
 end
 plot(rforphi,phir)
 title('Phi')
 figure
-plot((rforphi(1:end-1)+rforphi(2:end))/2,diff(phir)./diff(rforphi))
+plot((rforphi(1:end-1)+rforphi(2:end))/2,-diff(phir)./diff(rforphi))
 title('E_r')
 Emaxcloud=1/(R-1)*deltaphicloud(ne);
 Emaxvacuum=1/(R-1)*deltaphivacuum
 
 rgrid=linspace(1e-3,16e-2,500);
 zgrid=linspace(-0.24,0.24,1000);
 GainedE=zeros(length(rgrid),length(zgrid));
 for i=1:length(zgrid)
     drgrid=(R-1)/(R+1)*besseli(1,2*pi*rgrid/width)*(width/2/pi)*(cos(2*pi*0/width)-cos(2*pi*zgrid(i)/width));
     GainedE(:,i)=-(Phivacuum(rgrid,a,b,phia,phib)-Phivacuum(rgrid+drgrid,a,b,phia,phib));
 end
 figure
 title(sprintf('\\phi_a=%3.2f kV \\phi_b=%3.2f kV R=%3.2f',phia/1e3,phib/1e3, R))
 surface(zgrid,rgrid,GainedE);
 xlabel('z')
 ylabel('r')
 c=colorbar;
 c.Label.String='E [eV]';
 meanE=mean(mean(abs(GainedE(20:375,ceil(length(zgrid)/2):floor(5*length(zgrid)/6))),2))/2
 %%
 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
 h(1)=plot(vpardav,vperdav,'r--','Displayname',sprintf('Dav H0=%#.2g eV P0=%#.2g kg m^2 s^{-1}',H0/qe,P0));
 h(2)=plot(vpargauss,vpergauss,'g--','Displayname',sprintf('Gauss Te=%#.2g eV',Te));
 legend(h)
 legend('location','northoutside')
 xlim([-1 1])
 ylim([-1 1])
 
 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
 h(1)=plot(vpardav,vperdav,'r--','Displayname',sprintf('Dav H0=%#.2g eV P0=%#.2g kg m^2 s^{-1}',H0/qe,P0));
 h(2)=plot(vpargauss,vpergauss,'g--','Displayname',sprintf('Gauss Te=%#.2g eV',Te));
 xlim([-1 1])
 ylim([-1 1])
 
 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])
 
 %%
 figure
 r=double(M.rgrid);
 r(r<1e-3)=1e-3;
 r(r>16e-2)=16e-2;
 phi=zeros(length(r),1);
 for i=1:length(r)
     phi(i)=Phi(r(i),14e-2,15e-2,0.001,16e-2,-60000,0,0);
 end
 plot(r,phi)
 
 phim=phi*ones(1,length(M.zgrid));
 
 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;
 q=-qe;
 
 phibp=1/log(b/a)*(ne/2*q/eps_0*rbm^2*log(b/rbp)*log(rbp/rbm)+(q*ne/2/eps_0*(rbp^2-rbm^2)*(log(rbp/a)-0.5)+phia)*log(b/rbp)+phib*log(rbp/a));
 phibm=1/log(b/a)*(-q/eps_0*ne/2*rbm^2*log(rbm/a)*log(rbp/rbm)+(q/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=-q/eps_0*ne/4*(r^2-rbp^2)+phibp+(phibp-phibm+q/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
 function phi=Phivacuum(r,a,b,phia,phib)
     phi=((phib-phia)*log(r)+phia*log(b)-phib*log(a))/log(b/a);
 end
\ No newline at end of file
diff --git a/matlab/splinevelocity.m b/matlab/splinevelocity.m
index da683ba..8810db2 100644
--- a/matlab/splinevelocity.m
+++ b/matlab/splinevelocity.m
@@ -1,20 +1,20 @@
 classdef splinevelocity < splinequantity
     
     properties
         Partspercell
     end
     
     methods
         function obj=splinevelocity(filename, dataset, knotsr, knotsz, femorder, vnorm, postscale, index)
             obj@splinequantity(filename, dataset, knotsr, knotsz, femorder, vnorm, postscale, index);
-            obj.Partspercell=splinequantity(filename, dataset, knotsr, knotsz, femorder, 1, postscale, 1);
+            obj.Partspercell=splinequantity(filename, dataset, knotsr, knotsz, femorder, 1, 1, 1);
         end
         
         function quantity=coeffs(obj,indices)
             scalefactgrid=coeffs@splinequantity(obj.Partspercell,indices);
             scalefactgrid=1./scalefactgrid;
             scalefactgrid(isinf(scalefactgrid))=0;
             quantity=coeffs@splinequantity(obj,indices).*scalefactgrid;
         end
     end
 end
\ No newline at end of file
diff --git a/matlab/timescales.m b/matlab/timescales.m
index a44527f..d415388 100644
--- a/matlab/timescales.m
+++ b/matlab/timescales.m
@@ -1,93 +1,97 @@
-Pn=1e-8;%mbar
+%%
+Pn=1e-3;%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
+%%
+Ee=[1000,10000,20000,50000]; %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;
+B0=0.28;
 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;
 
+weight=3.83e6;
+ncreated=5e17*nb*sigion*ve*(0.030*2*pi*(0.078^2-0.076^2))
+nbsimucr=ncreated/weight
 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);
+taucycl=2*pi*me/qe/B0;
+taupe=2*pi*sqrt(eps_0*me/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