diff --git a/matlab/analyse_espicfields.m b/matlab/analyse_espicfields.m
new file mode 100644
index 0000000..97013ba
--- /dev/null
+++ b/matlab/analyse_espicfields.m
@@ -0,0 +1,173 @@
+%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