%% 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);
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;
pot=mean(M.pot(:,:,fieldstart:fieldend),3);
brillratio=2*dens*M.me./(M.eps_0*M.Bz'.^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);
[~,h1]=contour(ax1,M.zgrid*1000,M.rgrid*1000,Blines',real(levels),'r-.','linewidth',1.5,'Displayname','Magnetic field lines');

[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);
f.PaperUnits='centimeters';
papsize=[20 16];
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;
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;
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;
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)';
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))

end


%% 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);