%file='teststablegauss_13_traject2.h5';
%file='unigauss.h5';
file='Dav_5e18fineloss_wr+.h5';
%file='stablegauss_2e19fine.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=floor(0.8*t2dlength);

tmin=2;
tmax=length(M.ekin);
f=figure();
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]')
grid on

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)*100),'h-')
catch
    semilogy(M.t0d(tmin:tmax),abs(M.eerr(tmin:tmax)/M.etot(2)*100),'h-')
end
xlabel('t [s]')
ylabel('E_{err} %')
xlim([M.t0d(tmin),M.t0d(tmax)])
grid on
[~, name, ~] = fileparts(M.file);
f.PaperUnits='centimeters';
papsize=[14 16];
f.PaperSize=papsize;
print(f,sprintf('%sEnergy',name),'-dpdf','-fillpage')

maxdens=max(max(mean(M.N(:,:,fieldstart:end),3)));



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

f=figure();
ax1=gca;
dens=mean(M.N(:,:,fieldstart:end),3);
h=surface(ax1,M.zgrid,M.rgrid,dens);
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,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')

f=figure();
ax1=gca;
brillratio=2*mean(M.N(:,:,fieldstart:end),3)*M.me./(M.eps_0*M.Bz'.^2);
surface(ax1,M.zgrid,M.rgrid,brillratio);
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)
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')

[R,Z]=meshgrid([0.5*(M.rgrid(2:end)+M.rgrid(1:end-1)); 0],M.zgrid);
Rinv=1./R;
Rinv(:,1)=0;
f=figure();
ax3=gca;
omegare=(mean(M.fluidUTHET(:,:,fieldstart:end),3).*Rinv');
surface(ax3,M.zgrid,M.rgrid,omegare)
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')


f=figure();
tstart=floor(0.8*size(M.tpart,1));
legtext=sprintf("t=%2.1f - %2.1f [ns]",M.tpart(tstart)*1e9,M.tpart(end)*1e9);
subplot(1,2,1)
H=M.H(:,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(:,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);
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(:,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')


BrillouinRatio=2*M.omepe^2/M.omece^2