diff --git a/matlab/Testflux.m b/matlab/Testflux.m
new file mode 100644
index 0000000..945861b
--- /dev/null
+++ b/matlab/Testflux.m
@@ -0,0 +1,101 @@
+%%
+rAthetpos=30;
+
+fieldstep=length(M.t2d);%-100:5:length(M.t2d);
+Vpar=M.fluidUR(:,:,fieldstep).*(M.Br./M.B)' + (M.Bz./M.B)'.*M.fluidUZ(:,:,fieldstep);
+N=M.N(:,:,fieldstep);
+Gammapar=mean(N.*Vpar,3);
+
+
+levels=M.rAthet(rAthetpos,floor(M.nz/2)+1)*[1 1];
+
+rAthetposend=find(levels(1)>=M.rAthet(:,1),1,'last');
+rleft=sqrt(levels(1)/(M.rAthet(rAthetposend,1)/M.rgrid(rAthetposend)^2));
+
+thefig=figure;
+sp(1)=subplot(1,3,1);
+contourf(M.zgrid,M.rgrid,mean(N,3),'edgecolor','none');
+hold on;
+contour(M.zgrid,M.rgrid,M.rAthet,levels,'r-','linewidth',1.5,'Displayname','Magnetic field lines');
+c=colorbar;
+c.Label.String='n_e [m^{-3}]';
+sp(2)=subplot(1,3,2);
+contourf(M.zgrid,M.rgrid,Gammapar,'edgecolor','none');
+hold on
+contour(M.zgrid,M.rgrid,M.rAthet,levels,'r-','linewidth',1.5,'Displayname','Magnetic field lines');
+c=colorbar;
+c.Label.String='\Gamma_{par} [m^{-2}s^{-1}]';
+sp(3)=subplot(1,3,3);
+contourf(M.zgrid,M.rgrid,mean(Vpar,3),'edgecolor','none');
+hold on
+contour(M.zgrid,M.rgrid,M.rAthet,levels,'r-','linewidth',1.5,'Displayname','Magnetic field line');
+c=colorbar;
+c.Label.String='V_{par} [ms^{-1}]';
+linkaxes(sp);
+xlabel(sp,'z [m]')
+ylabel(sp,'r [m]')
+ylim(sp,[M.rgrid(1) M.rgrid(sum(M.nnr(1:2)))])
+M.savegraph(thefig,sprintf('%s/%s_finalwellfluxline',M.folder,M.name),[12,10])
+
+
+%%
+dN=1e13*M.weight;
+dN=0;
+rp=9.85e-3;
+rm=7e-3;
+V=2*pi*(rp^2-rm^2);
+n2=mean(mean(mean(N(rAthetposend+(0:1),[1 end],:),3)));
+vpar2=(dN/V/n2)^2;
+actvpar2=mean(mean(abs(mean(Vpar(rAthetposend+(0:1),[1 end],:),3))))^2;
+b=M.rgrid(end);
+a=M.rgrid(1);
+vd1=((M.potout-M.potinn)*M.phinorm/M.rgrid(rAthetpos))/M.Bz(floor(M.nz/2)+1,rAthetpos)/log(b/a);
+%vd1=-M.Er(rAthetpos+5,M.nz/2+1,end)/M.Bz(M.nz/2+1,rAthetpos+5)
+vinit=sqrt(M.kb*10000/M.me)
+%vinit=sqrt(80*M.qe/M.me)
+vd2=((M.potout-M.potinn)*M.phinorm/rleft)/M.Bz(1,rAthetposend)/log(b/a);
+deltaphi=-0.5*M.me/M.qe*(vpar2-(vd1+vinit)^2*(1-M.Rcurv))
+ratioparper=vpar2/((vd1+vinit)^2*(1-M.Rcurv))
+
+%%
+
+
+thefig=figure;
+%fieldstep=fieldstep(end);
+plotaxes(1)=subplot(1,2,1,'Parent',thefig);
+plotaxes(2)=subplot(1,2,2,'Parent',thefig);
+dens=mean(M.N(:,:,fieldstep),3);
+model=M.potentialwellmodel(fieldstep);
+z=model.z;
+r=model.r;
+pot=model.pot;
+rathet=model.rathet;
+[Zmesh,Rmesh]=meshgrid(M.zgrid,M.rAthet(:,floor(M.nz/2)+1));
+densplot=zeros(size(Zmesh,1),size(Zmesh,2),length(fieldstep));
+for i=1:length(fieldstep)
+densplot(:,:,i)=griddata(Zmesh,M.rAthet,dens,Zmesh,Rmesh,'natural');
+end
+plot(plotaxes(1),M.zgrid,mean(densplot(rAthetpos,1:end,:),3));
+xlim(plotaxes(1),[M.zgrid(1) M.zgrid(end)])
+xlabel(plotaxes(1),'z [m]')
+ylabel(plotaxes(1),'n [m^{-3}]')
+title(plotaxes(1),'Density')
+plotpot=zeros(size(Zmesh,1),size(Zmesh,2),length(fieldstep));
+for i=1:length(fieldstep)
+plotpot(:,:,i)=griddata(z,rathet,pot(:,i),Zmesh,Rmesh,'natural');
+end
+plotpot=mean(plotpot,3);
+plot(plotaxes(2),M.zgrid(1:end),plotpot(rAthetpos,1:end));
+hold on
+plot(plotaxes(2),M.zgrid(1:end),deltaphi*ones(size(M.zgrid)),'k--','displayname','Analytical')
+xlim(plotaxes(2),[M.zgrid(1) M.zgrid(end)])
+xlabel(plotaxes(2),'z [m]')
+ylabel(plotaxes(2),'depth [eV]')
+title(plotaxes(2),'Well')
+sgtitle(thefig,sprintf('r(0)=%1.2f [mm] t= %1.2f [ns]',M.rgrid(rAthetpos)*1e3,M.t2d(fieldstep(end))*1e9))
+
+M.savegraph(thefig,sprintf('%s/%s_finalwell',M.folder,M.name),[15,10])
+
+
+
+
diff --git a/matlab/analyse_espicfields.m b/matlab/analyse_espicfields.m
index 0043573..230a3f1 100644
--- a/matlab/analyse_espicfields.m
+++ b/matlab/analyse_espicfields.m
@@ -1,337 +1,343 @@
 % 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-500);%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.displayenergy
 
+
+
+%% Radial profile 
+try
+    M.displayrprofile(deltat)
+catch
+end
+
+
 %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
-try
-    M.displayrprofile(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)
 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',real(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',real(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))
 
 %% 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))
 
 %% VTHET
 f=figure('Name', sprintf('%s Vthet',name));
 ax3=gca;
 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)])
 ylim(ax3,[M.rgrid(1) M.rgrid(rgridend)])
 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(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_vthet',name),'-dpdf','-fillpage')
 savefig(f,sprintf('%sfluid_vthet',name))
 
 %% VDRIFT
 f=figure('Name', sprintf('%s V drift',name));
 ax3=gca;
 vdrift=-mean(M.Er(:,:,fieldstart:end),3).*M.Bz'./(M.B.^2)';
 v=abs(vdrift-VTHET);
 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)])
 ylim(ax3,[M.rgrid(1) M.rgrid(rgridend)])
 colormap(ax3,'jet')
 c = colorbar(ax3);
 set(ax3,'zscale','log')
 set(ax3,'colorscale','log')
 c.Label.String= '|V_{drift}| [m/s]';
 %c.Limits=[min(M.fluidUTHET(:)) max(M.fluidUTHET(:))];
 %caxis(ax3,[min(M.fluidUTHET(:)) max(M.fluidUTHET(:))])
 title(ax3,sprintf('Azimuthal drift velocity 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_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
 
 %% 0D diagnostics
-BrillouinRatio=2*M.omepe^2/M.omece^2
+BrillouinRatioth=2*M.omepe^2/M.omece^2
+BrillouinMax=max(max(brillratio))
+
 NpartsTrapped=mean(M.nbparts(end-10:end))*M.msim/M.me
 dblengthiter=fieldstart:length(M.t2d);
 
+
 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)))
 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)))
 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.displaypotentialwell(t2dlength);
 
 
 
diff --git a/matlab/distribution_function.m b/matlab/distribution_function.m
index cbb6a4d..5921055 100644
--- a/matlab/distribution_function.m
+++ b/matlab/distribution_function.m
@@ -1,254 +1,357 @@
 %% Show the particles velocity histogram at the position set using ginput
 % The histogram is compared between timesteppart and end
 
 timesteppart=length(M.tpart);
 
 fig=figure;
 ax1=gca;
 
 timestepN=find(M.tpart(end)==M.t2d,1);
 Ndistrib=M.N(:,:,timestepN);
 h=contourf(ax1,M.zgrid,M.rgrid,Ndistrib);
 %set(h, 'EdgeColor', 'none');
 hold on
 [r,z]=find(Ndistrib~=0);
 xlim(ax1,[M.zgrid(min(z)) M.zgrid(max(z))])
 ylim(ax1,[M.rgrid(min(r)) M.rgrid(max(r))])
 xlabel(ax1,'Z [m]')
 ylabel(ax1,'R [m]')
 c = colorbar(ax1);
 c.Label.String= 'n [m^{-3}]';
 view(ax1,2)
 %set(ax1,'colorscale','log')
 
 
-[x,y]=ginput(1);
-Zindex=find(x>M.zgrid,1,'last');
-Rindex=find(y>M.rgrid,1,'last');
+% [x,y]=ginput(1);
+% Zindex=find(x>M.zgrid,1,'last');
+% Rindex=find(y>M.rgrid,1,'last');
+Rindex=29
+Zindex=193;
 
 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;
 
 Rp=M.R(:,1,false);
 Zp=M.Z(:,1,false);
 Indices=Rp>=M.rgrid(Rindex) & Rp<M.rgrid(Rindex+1) & Zp>=M.zgrid(Zindex) & Zp<M.zgrid(Zindex+1);
 
 Vr=M.VR(Indices,1,false)/M.vlight;
 Vz=M.VZ(Indices,1,false)/M.vlight;
 Vthet=M.VTHET(Indices,1,false)/M.vlight;
 Vperp=sqrt(Vthet.^2+Vr.^2);
 Rp=M.R(Indices,1,false);
 Zp=M.Z(Indices,1,false);
 Thetp=M.THET(Indices,1,false);
 
 Rend=M.R(:,timesteppart,false);
 Zend=M.Z(:,timesteppart,false);
 Indicesend=Rend>=M.rgrid(Rindex) & Rend<M.rgrid(Rindex+1) & Zend>=M.zgrid(Zindex) & Zend<M.zgrid(Zindex+1);
 
 Vrend=M.VR(Indicesend,timesteppart,false)/M.vlight;
 Vzend=M.VZ(Indicesend,timesteppart,false)/M.vlight;
 Vthetend=M.VTHET(Indicesend,timesteppart,false)/M.vlight;
 Vperpend=sqrt(Vthetend.^2+Vrend.^2);
 Rend=M.R(Indicesend,timesteppart,false);
 Zend=M.Z(Indicesend,timesteppart,false);
 Thetend=M.THET(Indicesend,timesteppart,false);
 
 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=histfit(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=histfit(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('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))
 
 
 %% Shows the particles phase space at position (rindex, zindex) 
 % and times tinit and timesteppart
 timeinit=1;
 timestepNinit=find(M.tpart(timeinit)==M.t2d);
 timesteppart=length(M.tpart);
 timestepNend=find(M.tpart(timesteppart)==M.t2d);
 
 Rp=M.R(:,timeinit,false);
 Zp=M.Z(:,timeinit,false);
 Indices=Rp>=M.rgrid(Rindex) & Rp<M.rgrid(Rindex+1) & Zp>=M.zgrid(Zindex) & Zp<M.zgrid(Zindex+1);
 
 costhet=M.Bz(Zindex,Rindex)/M.B(Zindex,Rindex);
 sinthet=M.Br(Zindex,Rindex)/M.B(Zindex,Rindex);
 
 Vr=M.VR(Indices,timeinit,false)/M.vlight;
 Vz=M.VZ(Indices,timeinit,false)/M.vlight;
 Vthet=M.VTHET(Indices,timeinit,false)/M.vlight;
 Vperp=M.Vperp(Indices,timeinit,false)/M.vlight;
 Vpar=M.Vpar(Indices,timeinit,false)/M.vlight;
 Rp=M.R(Indices,timeinit,false);
 Zp=M.Z(Indices,timeinit,false);
 Thetp=M.THET(Indices,timeinit,false);
 
 Rend=M.R(:,timesteppart,false);
 Zend=M.Z(:,timesteppart,false);
 Indicesend=Rend>=M.rgrid(Rindex) & Rend<M.rgrid(Rindex+1) & Zend>=M.zgrid(Zindex) & Zend<M.zgrid(Zindex+1);
 
 Vrend=M.VR(Indicesend,timesteppart,false)/M.vlight;
 Vzend=M.VZ(Indicesend,timesteppart,false)/M.vlight;
 Vthetend=M.VTHET(Indicesend,timesteppart,false)/M.vlight;
 Vperpend=M.Vperp(Indicesend,timesteppart,false)/M.vlight;
 Vparend=M.Vpar(Indicesend,timesteppart,false)/M.vlight;
 Rend=Rend(Indicesend);
 Zend=Zend(Indicesend);
 Thetend=M.THET(Indicesend,timesteppart,false);
 
 legendinit=sprintf('t=%1.3g [s]',M.tpart(timeinit));
 legendend=sprintf('t=%1.3g [s]',M.tpart(timesteppart));
 
 Tpar=std(0.5*M.me*Vpar.^2*M.vlight^2)/M.qe
 Tparend=std(0.5*M.me*Vparend.^2*M.vlight^2)/M.qe
 Tperp=std(0.5*M.me*Vperp.^2*M.vlight^2)/M.qe
 Tperpend=std(0.5*M.me*Vperpend.^2*M.vlight^2)/M.qe
 
 f=figure;
-subplot(2,2,1)
+subplot(1,2,1)
 p(1)=scatter(Vpar,Vperp,'.','displayname',legendinit);
 hold on
 p(2)=scatter(Vparend,Vperpend,'.','displayname',legendend);
-xlabel('\beta_{par}')
-ylabel('\beta_\perp')
+xlabel('v_{par}/c')
+ylabel('v_\perp/c')
 legend('location','southoutside')
 
-R=M.Rcurv;
+zleftlim=1;
+
 vpar=linspace(1,M.vlight,1000);
-dN=1e12*M.weight;
+dN=1e13*M.weight;
+dN=0;
 rp=9.85e-3;
 rm=7e-3;
 V=2*pi*(rp^2-rm^2);
 N=M.N(:,:,timestepNend);
-n2=mean(mean(N(25,[1 end],:),3));
+levels=M.rAthet(Rindex,Zindex)*[1 1];
+rAthetposend=find(levels(1)>=M.rAthet(:,zleftlim),1,'last');
+rleft=sqrt(levels(1)/(M.rAthet(rAthetposend,1)/M.rgrid(rAthetposend)^2));
+rposleft=find(M.rgrid>=rleft,1,'first');
+n2=mean(mean(N(rposleft,[1 end],:),3));
 vpar2=(dN/V/n2)^2;
 if n2==0
     vpar2=0;
 end
 b=M.rgrid(end);
 a=M.rgrid(1);
 vd1=((M.potout-M.potinn)*M.phinorm/M.rgrid(Rindex))/M.Bz(Zindex,Rindex)/log(b/a);
-vd1=-M.Er(Rindex,Zindex,timestepNend)/M.Bz(Zindex,Rindex);
+vd1=-M.Er(Rindex,Zindex,1)/M.Bz(Zindex,Rindex);
+%vd1=-M.Er(Rindex,Zindex,timestepNend)/M.Bz(Zindex,Rindex);
 vinit=sqrt(M.kb*10000/M.me);
 %vinit=sqrt(M.qe*80/M.me)
 %vinit=sqrt(120*M.qe/M.me)
 %vd2=((M.potout-M.potinn)*M.phinorm/M.rgrid(25))/M.Bz(1,25)/log(b/a);
-deltaphi=-0.5*M.me/M.qe*(vpar2-(vd1+vinit)^2*(1-M.Rcurv))
+[Zmesh,Rmesh]=meshgrid(M.zgrid,M.rAthet(:,floor(M.nz/2)+1));
+Rcurv=griddata(Zmesh,M.rAthet,M.B',M.zgrid(zleftlim),M.rAthet(Rindex,Zindex),'natural')/M.B(Zindex,Rindex);
+deltaphicomp=-0.5*M.me/M.qe*(vpar2-(vd1+vinit)^2*(1-M.Rcurv))
+ratioparper=vpar2/((vd1+vinit)^2*(1-M.Rcurv))
 
 Ndistrib=M.N(:,:,timestepNend);
 model=M.potentialwellmodel(timestepNend);
 z=model.z;
 r=model.r;
 pot=model.pot;
 rathet=model.rathet;
-zleftlim=1;
-Zmesh=[M.zgrid(zleftlim) M.zgrid(Zindex)];
-Psimesh=M.rAthet(Rindex,Zindex)*[1 1];
-deltaphi=-diff(griddata(z,rathet,pot,Zmesh,Psimesh,'natural'));
+Zeval=[M.zgrid(zleftlim) M.zgrid(Zindex)];
+Psieval=M.rAthet(Rindex,Zindex)*[1 1];
+phis=griddata(Zmesh,M.rAthet,M.pot(:,:,end),Zeval,Psieval,'natural');
+deltaphi=-diff(phis)
 
 %deltaphi=-215;
-[Zmesh,Rmesh]=meshgrid(M.zgrid,M.rAthet(:,1));
 R=griddata(Zmesh,M.rAthet,M.B',M.zgrid(zleftlim),M.rAthet(Rindex,Zindex),'natural')/M.B(Zindex,Rindex);
-vper=sqrt(-2*M.qe/M.me*deltaphi/(R-1)+vpar.^2/(R-1))/M.vlight;
+vper=sqrt(2*M.qe/M.me*deltaphi/(R-1)+vpar.^2/(R-1))/M.vlight;
 vpar=vpar/M.vlight;
 vpar=vpar(real(vper)~=0);
 vper=vper(real(vper)~=0);
 xlimits=xlim;
 ylimits=ylim;
-p(3)=plot(vpar,vper,'b-','displayname','Energy boundary');
+p(3)=plot(vpar,vper,'b-','displayname','Loss parabolla simul');
 hold on
 plot(-vpar,vper,'b-')
+
+vpar=linspace(1,M.vlight,1000);
+vper=sqrt(-2*M.qe/M.me*deltaphicomp/(R-1)+vpar.^2/(R-1))/M.vlight;
+vpar=vpar/M.vlight;
+vpar=vpar(real(vper)~=0);
+vper=vper(real(vper)~=0);
+p(4)=plot(vpar,vper,'k--','displayname','Loss parabolla prediction');
+hold on
+plot(-vpar,vper,'k--')
 xlim(xlimits)
 ylim(ylimits)
 legend(p)
 axis equal
 
-subplot(2,2,2)
-scatter(Zp,Vpar,'.','displayname','Init')
-hold on
-scatter(Zend,Vparend,'.','displayname','End')
-xlabel('z [m]')
-ylabel('\beta_{par}')
-
-subplot(2,2,3)
-scatter(Vperp,Rp,'.','displayname','Init')
-hold on
-scatter(Vperpend,Rend,'.','displayname','End')
-xlabel('\beta_\perp')
-ylabel('R [m]')
-
-ax1=subplot(2,2,4);
+% subplot(2,2,2)
+% scatter(Zp,Vpar,'.','displayname','Init')
+% hold on
+% scatter(Zend,Vparend,'.','displayname','End')
+% xlabel('z [m]')
+% ylabel('\beta_{par}')
+% 
+% subplot(2,2,3)
+% scatter(Vperp,Rp,'.','displayname','Init')
+% hold on
+% scatter(Vperpend,Rend,'.','displayname','End')
+% xlabel('\beta_\perp')
+% ylabel('R [m]')
+
+ax1=subplot(1,2,2);
 h=contourf(ax1,M.zgrid,M.rgrid,Ndistrib);
 hold on
 [r,z]=find(Ndistrib~=0);
 xlim(ax1,[M.zgrid(min(z)) M.zgrid(max(z))])
 ylim(ax1,[M.rgrid(min(r)) M.rgrid(max(r))])
 xlabel(ax1,'Z [m]')
 ylabel(ax1,'R [m]')
 c = colorbar(ax1);
 c.Label.String= 'n [m^{-3}]';
 Zx=(M.zgrid(Zindex));
 Rx=(M.rgrid(Rindex));
 
 plot(ax1,Zx,Rx,'rx','Markersize',12);
 
-sgtitle(sprintf('r=%1.3g [m] z=%1.3g [m] \\Delta\\phi=%1.3g[kV]',M.rgrid(Rindex),M.zgrid(Zindex),(M.potout-M.potinn)*M.phinorm))
+sgtitle(sprintf('r=%1.3g [m] z=%1.3g [m] \\Delta\\phi=%1.3g[kV] R=%1.3g',M.rgrid(Rindex),M.zgrid(Zindex),(M.potout-M.potinn)*M.phinorm,M.Rcurv))
+
+M.savegraph(f,sprintf('%s/%s_phasespaceR%dZ%dpos',M.folder,M.name,Rindex,Zindex),[16,12])
+
+
+%%
+f=figure;
+p(1)=scatter(Vpar,Vperp,'.','displayname','Initial');
+hold on
+p(2)=scatter(Vparend,Vperpend,'.','displayname','Final');
+xlabel('v_{par}/c')
+ylabel('v_\perp/c')
+
+
+zleftlim=1;
+
+vpar=linspace(1,M.vlight,1000);
+dN=1e13*M.weight;
+dN=0;
+rp=9.85e-3;
+rm=7e-3;
+V=2*pi*(rp^2-rm^2);
+N=M.N(:,:,timestepNend);
+levels=M.rAthet(Rindex,Zindex)*[1 1];
+rAthetposend=find(levels(1)>=M.rAthet(:,zleftlim),1,'last');
+rleft=sqrt(levels(1)/(M.rAthet(rAthetposend,1)/M.rgrid(rAthetposend)^2));
+rposleft=find(M.rgrid>=rleft,1,'first');
+n2=mean(mean(N(rposleft,[1 end],:),3));
+vpar2=(dN/V/n2)^2;
+if n2==0
+    vpar2=0;
+end
+b=M.rgrid(end);
+a=M.rgrid(1);
+vd1=((M.potout-M.potinn)*M.phinorm/M.rgrid(Rindex))/M.Bz(Zindex,Rindex)/log(b/a);
+vd1=-M.Er(Rindex,Zindex,1)/M.Bz(Zindex,Rindex);
+%vd1=-M.Er(Rindex,Zindex,timestepNend)/M.Bz(Zindex,Rindex);
+vinit=sqrt(M.kb*10000/M.me);
+%vinit=sqrt(M.qe*80/M.me)
+%vinit=sqrt(120*M.qe/M.me)
+%vd2=((M.potout-M.potinn)*M.phinorm/M.rgrid(25))/M.Bz(1,25)/log(b/a);
+[Zmesh,Rmesh]=meshgrid(M.zgrid,M.rAthet(:,floor(M.nz/2)+1));
+Rcurv=griddata(Zmesh,M.rAthet,M.B',M.zgrid(zleftlim),M.rAthet(Rindex,Zindex),'natural')/M.B(Zindex,Rindex);
+deltaphicomp=-0.5*M.me/M.qe*(vpar2-(vd1+vinit)^2*(1-M.Rcurv))
+ratioparper=vpar2/((vd1+vinit)^2*(1-M.Rcurv))
+
+Ndistrib=M.N(:,:,timestepNend);
+model=M.potentialwellmodel(timestepNend);
+z=model.z;
+r=model.r;
+pot=model.pot;
+rathet=model.rathet;
+Zeval=[M.zgrid(zleftlim) M.zgrid(Zindex)];
+Psieval=M.rAthet(Rindex,Zindex)*[1 1];
+phis=griddata(Zmesh,M.rAthet,M.pot(:,:,end),Zeval,Psieval,'natural');
+deltaphi=-diff(phis)
+
+%deltaphi=-215;
+R=griddata(Zmesh,M.rAthet,M.B',M.zgrid(zleftlim),M.rAthet(Rindex,Zindex),'natural')/M.B(Zindex,Rindex);
+vper=sqrt(2*M.qe/M.me*deltaphi/(R-1)+vpar.^2/(R-1))/M.vlight;
+vpar=vpar/M.vlight;
+vpar=vpar(real(vper)~=0);
+vper=vper(real(vper)~=0);
+xlimits=xlim;
+ylimits=ylim;
+p(3)=plot(vpar,vper,'b-','displayname','Simulation');
+hold on
+plot(-vpar,vper,'b-')
+
+vpar=linspace(1,M.vlight,1000);
+vper=sqrt(-2*M.qe/M.me*deltaphicomp/(R-1)+vpar.^2/(R-1))/M.vlight;
+vpar=vpar/M.vlight;
+vpar=vpar(real(vper)~=0);
+vper=vper(real(vper)~=0);
+p(4)=plot(vpar,vper,'k--','displayname','Prediction');
+hold on
+plot(-vpar,vper,'k--')
+xlim(xlimits)
+ylim(ylimits)
+legend(p)
+axis equal
+legend('location','northeast')
+
+title(sprintf('r=%1.3g [m] z=%1.3g [m] \\Delta\\phi=%1.3g[kV] R=%1.3g',M.rgrid(Rindex),M.zgrid(Zindex),(M.potout-M.potinn)*M.phinorm,M.Rcurv))
 
-M.savegraph(f,sprintf('%s/%s_phasespaceR%dZ%d',M.folder,M.name,Rindex,Zindex),[18,14])
\ No newline at end of file
+M.savegraph(f,sprintf('%s/%s_phasespaceR%dZ%d',M.folder,M.name,Rindex,Zindex),[12,8])
\ No newline at end of file