% 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.plotEnergy

%fieldstart=2300; fieldend=2500;
dens=mean(M.N(:,:,fieldstart:end),3);
pot=mean(M.pot(:,:,fieldstart:end),3);
brillratio=2*dens*M.me./(M.eps_0*M.Bz'.^2);
[R,Z]=meshgrid(M.rgrid,M.zgrid);
Rinv=1./R;
Rinv(:,1)=0;
VTHET=mean(M.fluidUTHET(:,:,fieldstart:end),3);
omegare=(VTHET.*Rinv');
maxdens=max(max(dens));



% f=figure();
% ax3=gca;
% surface(ax3,M.zgrid(1:end-1),M.rgrid(1:end-1),(squeeze(mean(M.Presstens(4,:,:,fieldstart:end),4)))*M.vlight)
% xlabel(ax3,'z [m]')
% ylabel(ax3,'r [m]')
% xlim(ax3,[M.zgrid(1) M.zgrid(end)])
% ylim(ax3,[M.rgrid(1) M.rgrid(50)])
% colormap(ax3,'jet')
% c = colorbar(ax3);
% c.Label.String= 'thermal v_\theta [m/s]';
% %c.Limits=[min(M.fluidUTHET(:)) max(M.fluidUTHET(:))];
% %caxis(ax3,[min(M.fluidUTHET(:)) max(M.fluidUTHET(:))])
% title(ax3,'Azimuthal velocity')
% %set(ax3,'colorscale','log')
% grid(ax3, 'on')
% view(ax3,2)
% f.PaperOrientation='landscape';
% [~, name, ~] = fileparts(M.file);
% f.PaperUnits='centimeters';
% papsize=[16 14];
% f.PaperSize=papsize;
% print(f,sprintf('%sfluid_thermtheta',name),'-dpdf','-fillpage')

%% Psi function
try
    M.displaypsi(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',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))

%% Brillouin
f=figure('Name', sprintf('%s Brillouin',name));
ax1=gca;
h=surface(ax1,M.zgrid,M.rgrid,brillratio);
set(h,'edgecolor','none');
xlim(ax1,[M.zgrid(1) M.zgrid(end)])
ylim(ax1,[M.rgrid(1) M.rgrid(rgridend)])
xlabel(ax1,'z [m]')
ylabel(ax1,'r [m]')
title(ax1,'Position')
c = colorbar(ax1);
c.Label.String= 'Brillouin Ratio';
view(ax1,2)
caxis([0 1.2])
title(ax1,sprintf('Brillouin Ratio t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(fieldstart),M.t2d(end),double(maxdens)))
f.PaperOrientation='landscape';
[~, name, ~] = fileparts(M.file);
f.PaperUnits='centimeters';
papsize=[16 14];
f.PaperSize=papsize;
print(f,sprintf('%sfluid_Brillouin',name),'-dpdf','-fillpage')
savefig(f,sprintf('%sfluid_Brillouin',name))

%% Omegar
f=figure('Name', sprintf('%s Omegar',name));
ax3=gca;
surface(ax3,M.zgrid,M.rgrid,omegare,'edgecolor','None')
xlabel(ax3,'z [m]')
ylabel(ax3,'r [m]')
xlim(ax3,[M.zgrid(1) M.zgrid(end)])
ylim(ax3,[M.rgrid(1) M.rgrid(rgridend)])
colormap(ax3,'jet')
c = colorbar(ax3);
c.Label.String= '|\omega_\theta| [1/s]';
%c.Limits=[min(M.fluidUTHET(:)) max(M.fluidUTHET(:))];
%caxis(ax3,[min(M.fluidUTHET(:)) max(M.fluidUTHET(:))])
title(ax3,sprintf('Azimuthal frequency t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(fieldstart),M.t2d(end),double(maxdens)))
%set(ax3,'colorscale','log')
grid(ax3, 'on')
view(ax3,2)
f.PaperOrientation='landscape';
[~, name, ~] = fileparts(M.file);
f.PaperUnits='centimeters';
papsize=[16 14];
f.PaperSize=papsize;
print(f,sprintf('%sfluid_omegar',name),'-dpdf','-fillpage')
savefig(f,sprintf('%sfluid_omegar',name))

%% HP
if(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
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.Wells(pot,[fieldstart t2dlength],maxdens)
M.Wells(M.pot(:,:,1),1, max(max(M.N(:,:,1))))