% 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-10);%floor(0.95*t2dlength);

tmin=2;
tmax=length(M.ekin);
[~, name, ~] = fileparts(M.file);
f=figure('Name', sprintf('%s Energy',name));
subplot(2,1,1)
plot(M.t0d(tmin:tmax),M.ekin(tmin:tmax),'o-',...
    M.t0d(tmin:tmax),M.epot(tmin:tmax),'d-',...
    M.t0d(tmin:tmax),M.etot(tmin:tmax),'h-',...
    M.t0d(tmin:tmax),M.eerr(tmin:tmax),'x--')
legend('ekin', 'epot', 'etot','eerr')
xlabel('Time [s]')
ylabel('Energies [J]')
%ylim([0 1.2])
grid on
xlimits=xlim();

%[~, rgridend]=min(abs(M.rgrid-0.045));
rgridend=sum(M.nnr(1:2));

subplot(2,1,2)
try
    semilogy(M.t0d(tmin:tmax),abs(M.eerr(tmin:tmax)./M.etot0(tmin:tmax)),'h-')
catch
    semilogy(M.t0d(tmin:tmax),abs(M.eerr(tmin:tmax)/M.etot(2)),'h-')
end
xlabel('t [s]')
ylabel('E_{err}/E_{tot}')
xlim(xlimits)
grid on
try
    yyaxis right
    plot(M.t0d(tmin:tmax),abs(M.nbparts(tmin:tmax)./M.nbparts(1)*100),'d--')
    ylabel('Nparts %')
    ylim([0 110])
catch
end

f.PaperUnits='centimeters';
papsize=[14 16];
f.PaperSize=papsize;
print(f,sprintf('%sEnergy',name),'-dpdf','-fillpage')
savefig(f,sprintf('%sEnergy',name))

%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;
VR=mean(M.fluidUR(:,:,fieldstart:end),3);
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)

%% Density
f=figure('Name', sprintf('%sfluid_dens',name));
ax1=gca;
h=surface(ax1,M.zgrid,M.rgrid,dens);
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;
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;
h=contourf(ax1,M.zgrid,M.rgrid,sign(M.qsim)*M.qe*pot+0.5*M.me*(VR.^2+VTHET.^2));
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('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 [J]';
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))

%% 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
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)
H=M.H(1:M.nbparts(1),1);
h1=histogram(H,20,'BinLimits',[min(H(:)) max(H(:))],'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9));
hold on
H=mean(M.H(1:M.nbparts(end),tstart:end),2);
h1=histogram(H,'BinWidth',h1.BinWidth,'DisplayName',legtext);
ylabel('counts')
xlabel('H [J]')
legend

subplot(1,2,2)
P=M.P(1:M.nbparts(1),1);
h2=histogram(P,20,'BinLimits',[min(P(:)) max(P(:))],'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9));
hold on
P=mean(M.P(1:M.nbparts(end),tstart:end),2);
histogram(P,'BinWidth',h2.BinWidth,'DisplayName',legtext);
ylabel('counts')
xlabel('P [kg\cdotm^2\cdots^{-1}]')
xlim(h2.BinLimits)
%clear P
%clear H
legend
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))


%% 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 temperatura')
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])

function displaypsi(M)
%% Psi function
[~, name, ~] = fileparts(M.file);
f=figure('Name', sprintf('%s Psi',name));
zpos=floor(length(M.zgrid)/2);
tinit=2;
tend=length(M.t2d);
deltat=200;
H0=3.2e-14;
P0=8.66e-25;
%H0=mean(M.H(:,end));
%P0=mean(M.P(:,end));
lw=1.5;
Mirrorratio=(M.Rcurv-1)/(M.Rcurv+1);
locpot=M.pot(:,zpos,tend);
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);
I=[I(1)-2; I(1)-1; I; I(end)+1; I(end)+2];
rq=linspace(M.rgrid(I(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/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