diff --git a/matlab/compute/mode_growth_meter.m b/matlab/compute/mode_growth_meter.m
index b8c229d..8c2bce8 100644
--- a/matlab/compute/mode_growth_meter.m
+++ b/matlab/compute/mode_growth_meter.m
@@ -1,100 +1,106 @@
 function [FIGURE] = mode_growth_meter(DATA,OPTIONS)
 
 NORMALIZED = OPTIONS.NORMALIZED;
 Nma   = OPTIONS.NMA; %Number moving average
 
-iz = OPTIONS.iz;
-[~,ikzf] = max(squeeze(mean(abs(squeeze(DATA.PHI(1,:,iz,:))),2)));
+switch OPTIONS.iz
+    case 'avg'
+        field = squeeze(mean(DATA.PHI,3));
+    otherwise
+        field = squeeze(DATA.PHI(:,:,OPTIONS.iz,:));
+end
 
 FRAMES = zeros(size(OPTIONS.TIME));
 for i = 1:numel(OPTIONS.TIME)
     [~,FRAMES(i)] =min(abs(OPTIONS.TIME(i)-DATA.Ts3D));
 end
 FRAMES = unique(FRAMES);
 t  = DATA.Ts3D(FRAMES);
 
+[~,ikzf] = max(mean(squeeze(abs(field(1,:,FRAMES))),2));
+
 FIGURE.fig = figure; set(gcf, 'Position',  [100 100 1200 700])
 FIGURE.FIGNAME = 'mode_growth_meter';
 for i = 1:2
     MODES_SELECTOR = i; %(2:Zonal, 1: NZonal, 3: ky=kx)
 
     if MODES_SELECTOR == 2
         if NORMALIZED
-            plt = @(x,ik) movmean(abs(squeeze(x(1,ik,iz,FRAMES)))./max(abs(squeeze(x(1,ik,iz,FRAMES)))),Nma);
+            plt = @(x,ik) movmean(abs(squeeze(x(1,ik,FRAMES)))./max(abs(squeeze(x(1,ik,FRAMES)))),Nma);
         else
-            plt = @(x,ik) movmean(abs(squeeze(x(1,ik,iz,FRAMES))),Nma);
+            plt = @(x,ik) movmean(abs(squeeze(x(1,ik,FRAMES))),Nma);
         end
         kstr = 'k_x';
         k = DATA.kx;
         MODESTR = 'Zonal modes';
     elseif MODES_SELECTOR == 1
         if NORMALIZED
-            plt = @(x,ik) movmean(abs(squeeze(x(ik,1,iz,FRAMES)))./max(abs(squeeze(x(ik,1,iz,FRAMES)))),Nma);
+            plt = @(x,ik) movmean(abs(squeeze(x(ik,1,FRAMES)))./max(abs(squeeze(x(ik,1,FRAMES)))),Nma);
         else
-            plt = @(x,ik) movmean(abs(squeeze(x(ik,1,iz,FRAMES))),Nma);
+            plt = @(x,ik) movmean(abs(squeeze(x(ik,1,FRAMES))),Nma);
         end
         kstr = 'k_y';
         k = DATA.ky;
         MODESTR = 'NZ modes';
     elseif MODES_SELECTOR == 3
         if NORMALIZED
-            plt = @(x,ik) movmean(abs(squeeze(x(ik,ik,iz,FRAMES)))./max(abs(squeeze(x(ik,ik,iz,FRAMES)))),Nma);
+            plt = @(x,ik) movmean(abs(squeeze(x(ik,ik,FRAMES)))./max(abs(squeeze(x(ik,ik,FRAMES)))),Nma);
         else
-            plt = @(x,ik) movmean(abs(squeeze(x(ik,ik,iz,FRAMES))),Nma);
+            plt = @(x,ik) movmean(abs(squeeze(x(ik,ik,FRAMES))),Nma);
         end
         kstr = 'k_y=k_x';
         k = DATA.ky;
         MODESTR = 'Diag modes';
     end 
 
     MODES = 1:numel(k);
     % MODES = zeros(size(OPTIONS.K2PLOT));
     % for i = 1:numel(OPTIONS.K2PLOT)
     %     [~,MODES(i)] =min(abs(OPTIONS.K2PLOT(i)-k));
     % end
 
 
     % plt = @(x,ik) abs(squeeze(x(1,ik,iz,FRAMES)))./max(abs(squeeze(x(1,ik,iz,FRAMES))));
 
     gamma = MODES;
     amp   = MODES;
     i_=1;
     for ik = MODES
-        gr = polyfit(t,log(plt(DATA.PHI,ik)),1);
+        gr = polyfit(t,log(plt(field,ik)),1);
         gamma(i_) = gr(1);
         amp(i_)   = gr(2);
         i_=i_+1;
     end
 
     %plot
     subplot(2,3,1+3*(i-1))
         [YY,XX] = meshgrid(t,fftshift(k));
-        pclr = pcolor(XX,YY,abs(plt(fftshift(DATA.PHI,MODES_SELECTOR),1:numel(k))));set(pclr, 'edgecolor','none');  hold on;
+        pclr = pcolor(XX,YY,abs(plt(fftshift(field,MODES_SELECTOR),1:numel(k))));set(pclr, 'edgecolor','none');  hold on;
         set(gca,'YDir','normal')
     %     xlim([t(1) t(end)]); %ylim([1e-5 1])
         xlabel(['$',kstr,'\rho_s$']); ylabel('$t c_s /\rho_s$');
         title(MODESTR)  
 
     subplot(2,3,2+3*(i-1))
         mod2plot = [2:OPTIONS.NMODES+1];
         for i_ = mod2plot
-            semilogy(t,plt(DATA.PHI,MODES(i_))); hold on;
+            semilogy(t,plt(field,MODES(i_))); hold on;
     %         semilogy(t,exp(gamma(i_).*t+amp(i_)),'--k')
         end
         if MODES_SELECTOR == 2
-            semilogy(t,plt(DATA.PHI,ikzf),'--k');
+            semilogy(t,plt(field,ikzf),'--k');
         end
         xlim([t(1) t(end)]); %ylim([1e-5 1])
         xlabel('$t c_s /\rho_s$'); ylabel(['$|\phi_{',kstr,'}|$']);
         title('Measure time window')
 
     subplot(2,3,3+3*(i-1))
         plot(k(MODES),gamma); hold on;
         plot(k(MODES(mod2plot)),gamma(mod2plot),'x')
         if MODES_SELECTOR == 2
             plot(k(ikzf),gamma(ikzf),'ok');
         end
         xlabel(['$',kstr,'\rho_s$']); ylabel('$\gamma$');
         title('Growth rates')
 end
 end
\ No newline at end of file
diff --git a/matlab/plot/plot_radial_transport_and_spacetime.m b/matlab/plot/plot_radial_transport_and_spacetime.m
index 93e0c5c..d052fc1 100644
--- a/matlab/plot/plot_radial_transport_and_spacetime.m
+++ b/matlab/plot/plot_radial_transport_and_spacetime.m
@@ -1,120 +1,122 @@
 function [FIGURE] = plot_radial_transport_and_spacetime(DATA, OPTIONS)
     %Compute steady radial transport
     tend = OPTIONS.TAVG_1; tstart = OPTIONS.TAVG_0;
     [~,its0D] = min(abs(DATA.Ts0D-tstart));
     [~,ite0D]   = min(abs(DATA.Ts0D-tend));
+    [~,its3D] = min(abs(DATA.Ts3D-tstart));
+    [~,ite3D]   = min(abs(DATA.Ts3D-tend));
     SCALE = DATA.scale;%(1/DATA.Nx/DATA.Ny)^2;
     Gx_infty_avg = mean(DATA.PGAMMA_RI(its0D:ite0D))*SCALE;
     Gx_infty_std = std (DATA.PGAMMA_RI(its0D:ite0D))*SCALE;
     Qx_infty_avg = mean(DATA.HFLUX_X(its0D:ite0D))*SCALE;
     Qx_infty_std = std (DATA.HFLUX_X(its0D:ite0D))*SCALE;
-    [~,ikzf] = max(squeeze(mean(abs(squeeze(DATA.PHI(:,1,1,:))),2)));
+    f_avg_z      = squeeze(mean(DATA.PHI(:,:,:,:),3));
+    [~,ikzf] = max(squeeze(mean(abs(f_avg_z(1,:,its3D:ite3D)),3)));
     Ns3D = numel(DATA.Ts3D);
     [KX, KY] = meshgrid(DATA.kx, DATA.ky);
     
     %% computations
 
     % Compute Gamma from ifft matlab
     Gx = zeros(DATA.Ny,DATA.Nx,numel(DATA.Ts3D));
-    for it = 1:numel(DATA.Ts3D)
-        for iz = 1:DATA.Nz
-            Gx(:,:,it)  = Gx(:,:,it) + ifourier_GENE(-1i*KY.*(DATA.PHI(:,:,iz,it)))...
-                          .*ifourier_GENE(DATA.DENS_I(:,:,iz,it));
-        end
-        Gx(:,:,it)  = Gx(:,:,it)/DATA.Nz;
-    end
+%     for it = 1:numel(DATA.Ts3D)
+%         for iz = 1:DATA.Nz
+%             Gx(:,:,it)  = Gx(:,:,it) + ifourier_GENE(-1i*KY.*(DATA.PHI(:,:,iz,it)))...
+%                           .*ifourier_GENE(DATA.DENS_I(:,:,iz,it));
+%         end
+%         Gx(:,:,it)  = Gx(:,:,it)/DATA.Nz;
+%     end
     Gx_t_mtlb = squeeze(mean(mean(Gx,1),2)); 
     % Compute Heat flux from ifft matlab
     Qx = zeros(DATA.Ny,DATA.Nx,numel(DATA.Ts3D));
-    for it = 1:numel(DATA.Ts3D)
-        for iz = 1:DATA.Nz
-            Qx(:,:,it)  = Qx(:,:,it) + ifourier_GENE(-1i*KY.*(DATA.PHI(:,:,iz,it)))...
-                          .*ifourier_GENE(DATA.TEMP_I(:,:,iz,it));
-        end
-        Qx(:,:,it)  = Qx(:,:,it)/DATA.Nz;
-    end
+%     for it = 1:numel(DATA.Ts3D)
+%         for iz = 1:DATA.Nz
+%             Qx(:,:,it)  = Qx(:,:,it) + ifourier_GENE(-1i*KY.*(DATA.PHI(:,:,iz,it)))...
+%                           .*ifourier_GENE(DATA.TEMP_I(:,:,iz,it));
+%         end
+%         Qx(:,:,it)  = Qx(:,:,it)/DATA.Nz;
+%     end
     Qx_t_mtlb = squeeze(mean(mean(Qx,1),2)); 
     % zonal vs nonzonal energies for phi(t)
 
     E_Zmode_SK       = zeros(1,Ns3D);
     E_NZmode_SK      = zeros(1,Ns3D);
     for it = 1:numel(DATA.Ts3D)
-        E_Zmode_SK(it)   = squeeze(DATA.ky(ikzf).^2.*abs(DATA.PHI(ikzf,1,1,it))^2);
-        E_NZmode_SK(it)  = squeeze(sum(sum(((1+KX.^2+KY.^2).*abs(DATA.PHI(:,:,1,it)).^2.*(KY~=0)))));
+        E_Zmode_SK(it)   = squeeze(DATA.ky(ikzf).^2.*abs(squeeze(f_avg_z(ikzf,1,it))).^2);
+        E_NZmode_SK(it)  = squeeze(sum(sum(((1+KX.^2+KY.^2).*abs(squeeze(f_avg_z(:,:,it))).^2.*(KY~=0)))));
     end
-    [~,its3D] = min(abs(DATA.Ts3D-tstart));
-    [~,ite3D]   = min(abs(DATA.Ts3D-tend));
+
 
 %% Figure    
 mvm = @(x) movmean(x,OPTIONS.NMVA);
     FIGURE.fig = figure; FIGURE.FIGNAME = ['ZF_transport_drphi','_',DATA.PARAMS]; set(gcf, 'Position',  [100, 100, 1000, 600])
     subplot(311)
     yyaxis left
         plot(mvm(DATA.Ts0D),mvm(DATA.PGAMMA_RI*SCALE),'DisplayName','$\langle n_i \partial_y\phi \rangle_y$'); hold on;
         plot(mvm(DATA.Ts3D),mvm(Gx_t_mtlb),'DisplayName','matlab comp.'); hold on;
         plot(DATA.Ts0D(its0D:ite0D),ones(ite0D-its0D+1,1)*Gx_infty_avg, '-k',...
             'DisplayName',['$\Gamma^{\infty} = $',num2str(Gx_infty_avg),'$\pm$',num2str(Gx_infty_std)]);
         ylabel('$\Gamma_x$')
         ylim([0,5*abs(Gx_infty_avg)]); xlim([DATA.Ts0D(1),DATA.Ts0D(end)]);
     yyaxis right
         plot(mvm(DATA.Ts0D),mvm(DATA.HFLUX_X*SCALE),'DisplayName','$\langle n_i \partial_y\phi \rangle_y$'); hold on;
         plot(mvm(DATA.Ts3D),mvm(Qx_t_mtlb),'DisplayName','matlab comp.'); hold on;
         plot(DATA.Ts0D(its0D:ite0D),ones(ite0D-its0D+1,1)*Qx_infty_avg, '-k',...
             'DisplayName',['$Q_x^{\infty} = $',num2str(Qx_infty_avg),'$\pm$',num2str(Qx_infty_std)]);
         ylabel('$Q_x$')  
         ylim([0,5*abs(Qx_infty_avg)]); xlim([DATA.Ts0D(1),DATA.Ts0D(end)]);
     grid on; set(gca,'xticklabel',[]); 
     title({DATA.param_title,...
         ['$\Gamma^{\infty} = $',num2str(Gx_infty_avg),'$, Q^{\infty} = $',num2str(Qx_infty_avg)]});
     %% radial shear radial profile
         % computation
     Ns3D = numel(DATA.Ts3D);
     [KX, KY] = meshgrid(DATA.kx, DATA.ky);
     plt = @(x) mean(x(:,:,:),1);
     kycut = max(DATA.ky);
     kxcut = max(DATA.kx);
     LP = (abs(KY)<kycut).*(abs(KX)<kxcut); %Low pass filter
     
     OPTIONS.NAME = OPTIONS.ST_FIELD;
     OPTIONS.PLAN = 'xy';
     OPTIONS.COMP = 'avg';
     OPTIONS.TIME = DATA.Ts3D;
     OPTIONS.POLARPLOT = 0;
     toplot = process_field(DATA,OPTIONS);
     f2plot = toplot.FIELD;
     
-    clim = max(max(max(abs(plt(f2plot(:,:,:))))));
+    clim = max(max(max(abs(plt(f2plot(:,:,its3D:ite3D))))));
     subplot(313)
         [TY,TX] = meshgrid(DATA.x,DATA.Ts3D(toplot.FRAMES));
         pclr = pcolor(TX,TY,squeeze(plt(f2plot))'); 
         set(pclr, 'edgecolor','none'); 
         legend(['$\langle ',OPTIONS.ST_FIELD,' \rangle_{y,z}$']) %colorbar;
         caxis(clim*[-1 1]);
         cmap = bluewhitered(256);
         colormap(cmap)
         xlabel('$t c_s/R$'), ylabel('$x/\rho_s$'); 
         if OPTIONS.INTERP
             shading interp
         end
     if strcmp(OPTIONS.ST_FIELD,'Gx')
         subplot(311)
         plot(DATA.Ts3D,squeeze(mean(plt(f2plot),1)));
     end
 %% Zonal vs NZonal energies    
     subplot(312)
     it0 = 1; itend = Ns3D;
     trange = toplot.FRAMES;
     plt1 = @(x) x;%-x(1);
     plt2 = @(x) x./max((x(:)));
     toplot = sum(squeeze(plt(f2plot)).^2,1); %ST from before
 %     plty = @(x) x(500:end)./max(squeeze(x(500:end)));
         yyaxis left
 %         plot(plt1(DATA.Ts3D(trange)),plt2(E_Zmode_SK(trange)),'DisplayName','$k_{zf}^2|\phi_{kzf}|^2$');
         plot(plt1(DATA.Ts3D(trange)),plt2(toplot(:)),'DisplayName','Sum $A^2$');
         ylim([-0.1, 1.5]); ylabel('$E_{Z}$')
         yyaxis right
         plot(plt1(DATA.Ts3D(trange)),plt2(E_NZmode_SK(trange)),'DisplayName','$(1+k^2)|\phi^k|^2$');
         xlim([DATA.Ts3D(it0), DATA.Ts3D(itend)]);
         ylim([-0.1, 1.5]); ylabel('$E_{NZ}$')
         xlabel('$t c_s/R$'); grid on; set(gca,'xticklabel',[]);% xlim([0 500]);
 end
\ No newline at end of file
diff --git a/wk/analysis_HeLaZ.m b/wk/analysis_HeLaZ.m
index 6fba72b..4a06a7a 100644
--- a/wk/analysis_HeLaZ.m
+++ b/wk/analysis_HeLaZ.m
@@ -1,205 +1,205 @@
 addpath(genpath([helazdir,'matlab'])) % ... add
 addpath(genpath([helazdir,'matlab/plot'])) % ... add
 addpath(genpath([helazdir,'matlab/compute'])) % ... add
 addpath(genpath([helazdir,'matlab/load'])) % ... add
 
 %% Load the results
 LOCALDIR  = [helazdir,'results/',outfile,'/'];
 MISCDIR   = ['/misc/HeLaZ_outputs/results/',outfile,'/'];
 system(['mkdir -p ',MISCDIR]);
 CMD = ['rsync ', LOCALDIR,'outputs* ',MISCDIR]; disp(CMD);
 system(CMD);
 % Load outputs from jobnummin up to jobnummax
-JOBNUMMIN = 10; JOBNUMMAX = 20;
+JOBNUMMIN = 15; JOBNUMMAX = 16;
 data = compile_results(MISCDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
 data.localdir = LOCALDIR;
 data.FIGDIR   = LOCALDIR;
 
 %% PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 default_plots_options
 disp('Plots')
 FMT = '.fig';
 
 if 1
 %% Space time diagramm (fig 11 Ivanov 2020)
 options.TAVG_0   = 0.98*data.Ts3D(end); 
 options.TAVG_1   = data.Ts3D(end); % Averaging times duration
 options.NMVA     = 1;              % Moving average for time traces
 % options.ST_FIELD = '\Gamma_x';          % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
 options.ST_FIELD = '\phi';          % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
 options.INTERP   = 1;
 fig = plot_radial_transport_and_spacetime(data,options);
 save_figure(data,fig)
 end
 
 if 0
 %% statistical transport averaging
 options.T = [16000 17000];
 fig = statistical_transport_averaging(data,options);
 end
 
 if 0
 %% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Options
 options.INTERP    = 1;
 options.POLARPLOT = 0;
 options.NAME      = '\phi';
 % options.NAME      = 'N_i^{00}';
 % options.NAME      = 'v_y';
 % options.NAME      = 'n_i^{NZ}';
 % options.NAME      = '\Gamma_x';
 % options.NAME      = 'n_i';
-options.PLAN      = 'xz';
+options.PLAN      = 'xy';
 % options.NAME      = 'f_e';
 % options.PLAN      = 'sx';
-options.COMP      = 'avg';
+options.COMP      = 9;
 % options.TIME      = dat.Ts5D;
-options.TIME      = 1250:1:1500;
+options.TIME      = [1800:1:2000];
 data.EPS          = 0.1;
 data.a = data.EPS * 2000;
 create_film(data,options,'.gif')
 end
 
 if 0
 %% 2D snapshots
 % Options
-options.INTERP    = 0;
+options.INTERP    = 1;
 options.POLARPLOT = 0;
 options.AXISEQUAL = 1;
 options.NAME      = '\phi';
 % options.NAME      = 'n_i';
 % options.NAME      = 'N_i^{00}';
 % options.NAME      = 'T_i';
 % options.NAME      = '\Gamma_x';
 % options.NAME      = 'k^2n_e';
-options.PLAN      = 'kxky';
+options.PLAN      = 'xy';
 % options.NAME      = 'f_i';
 % options.PLAN      = 'sx';
 options.COMP      = 'avg';
-options.TIME      = [1200 1300 1400 1500];
+options.TIME      = [1600 1800 2000];
 data.a = data.EPS * 2e3;
 fig = photomaton(data,options);
 save_figure(data,fig)
 end
 
 if 0
 %% 3D plot on the geometry
 options.INTERP    = 1;
 options.NAME      = 'n_i';
 options.PLANES    = [1:2:12];
 options.TIME      = [60];
 options.PLT_MTOPO = 1;
 data.rho_o_R      = 2e-3; % Sound larmor radius over Machine size ratio
 fig = show_geometry(data,options);
 save_figure(data,fig)
 end
 
 if 0
 %% Kinetic distribution function sqrt(<f_a^2>xy) (GENE vsp)
 % options.SPAR      = linspace(-3,3,64)+(6/127/2);
 % options.XPERP     = linspace( 0,6,64);
 options.SPAR      = gene_data.vp';
 options.XPERP     = gene_data.mu';
 options.iz        = 9;
 options.T         = 30;
 options.PLT_FCT   = 'contour';
 options.ONED      = 0;
 options.non_adiab = 1;
 options.SPECIE    = 'i';
 options.RMS       = 1; % Root mean square i.e. sqrt(sum_k|f_k|^2) as in Gene
 fig = plot_fa(data,options);
 save_figure(data,fig)
 end
 
 if 0
 %% Hermite-Laguerre spectrum
 % options.TIME = 'avg';
 options.P2J        = 1;
 options.ST         = 0;
 options.PLOT_TYPE  = 'space-time';
 options.NORMALIZED = 1;
 options.JOBNUM     = 0;
 options.TIME       = [1300 1400];
 options.specie     = 'i';
 options.compz      = 'avg';
 fig = show_moments_spectrum(data,options);
 % fig = show_napjz(data,options);
 save_figure(data,fig)
 end
 
 if 0
 %% Time averaged spectrum
 options.TIME   = 1000:1200;
 options.NORM   =1;
 options.NAME   = '\phi';
 % options.NAME      = 'n_i';
 % options.NAME      ='\Gamma_x';
 options.PLAN   = 'kxky';
 options.COMPZ  = 'avg';
 options.OK     = 0;
 options.COMPXY = 'avg';
 options.COMPT  = 'avg';
 options.PLOT   = 'semilogy';
 fig = spectrum_1D(data,options);
 % save_figure(data,fig)
 end
 
 if 0
 %% 1D real plot
 options.TIME   = [50 100 200];
 options.NORM   = 0;
 options.NAME   = '\phi';
 % options.NAME      = 'n_i';
 % options.NAME      ='\Gamma_x';
 % options.NAME      ='s_y';
 options.COMPX  = 'avg';
 options.COMPY  = 'avg';
 options.COMPZ  = 1;
 options.COMPT  = 1;
 options.MOVMT  = 1;
 fig = real_plot_1D(data,options);
 % save_figure(data,fig)
 end
 
 if 0
 %% Mode evolution
 options.NORMALIZED = 0;
 options.K2PLOT = 1;
-options.TIME   = 1200:1300;
+options.TIME   = 1800:2000;
 options.NMA    = 1;
 options.NMODES = 15;
-options.iz     = 9;
+options.iz     = 'avg';
 fig = mode_growth_meter(data,options);
 save_figure(data,fig)
 end
 
 if 0
 %% ZF caracteristics (space time diagrams)
 TAVG_0 = 1200; TAVG_1 = 1500; % Averaging times duration
 % chose your field to plot in spacetime diag (uzf,szf,Gx)
 fig = ZF_spacetime(data,TAVG_0,TAVG_1);
 save_figure(data,fig)
 end
 
 if 0
 %% Metric infos
 fig = plot_metric(data);
 end
 
 if 0
 %% linear growth rate for 3D fluxtube
 trange = [0 100];
 nplots = 1;
 lg = compute_fluxtube_growth_rate(data,trange,nplots);
 end
 
 if 0
 %% linear growth rate for 3D Zpinch
 trange = [5 15];
 options.keq0 = 1; % chose to plot planes at k=0 or max
 options.kxky = 1;
 options.kzkx = 0;
 options.kzky = 1;
 [lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options);
 save_figure(data,fig)
 end
diff --git a/wk/analysis_gene.m b/wk/analysis_gene.m
index 281885d..86e064c 100644
--- a/wk/analysis_gene.m
+++ b/wk/analysis_gene.m
@@ -1,113 +1,113 @@
 % folder = '/misc/gene_results/shearless_cyclone/miller_output_1.0/';
 % folder = '/misc/gene_results/shearless_cyclone/miller_output_0.8/';
 folder = '/misc/gene_results/shearless_cyclone/s_alpha_output_1.0/';
 % folder = '/misc/gene_results/shearless_cyclone/linear_s_alpha_CBC_100/';
 % folder = '/misc/gene_results/shearless_cyclone/s_alpha_output_0.5/';
 % folder = '/misc/gene_results/shearless_cyclone/LD_s_alpha_output_1.0/';
 % folder = '/misc/gene_results/shearless_cyclone/LD_s_alpha_output_0.8/';
 % folder = '/misc/gene_results/HP_fig_2b_mu_5e-2/';
 % folder = '/misc/gene_results/HP_fig_2c_mu_5e-2/';
 gene_data = load_gene_data(folder);
 gene_data = invert_kxky_to_kykx_gene_results(gene_data);
 if 1
 %% Space time diagramm (fig 11 Ivanov 2020)
 options.TAVG_0   = 0.8*gene_data.Ts3D(end);
 options.TAVG_1   = gene_data.Ts3D(end); % Averaging times duration
 options.NMVA     = 1;              % Moving average for time traces
 options.ST_FIELD = '\phi';          % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x, Q_x)
 options.INTERP   = 1;
-gene_data.FIGDIR = folder;;
+gene_data.FIGDIR = folder;
 fig = plot_radial_transport_and_spacetime(gene_data,options);
 save_figure(gene_data,fig)
 end
 
 if 0
 %% 2D snapshots
 % Options
 options.INTERP    = 0;
 options.POLARPLOT = 0;
 options.AXISEQUAL = 1;
 % options.NAME      = 'Q_x';
 options.NAME      = '\phi';
 % options.NAME      = 'T_i';
 % options.NAME      = '\Gamma_x';
 % options.NAME      = 'k^2n_e';
-options.PLAN      = 'kxky';
+options.PLAN      = 'xy';
 % options.NAME      ='f_e';
 % options.PLAN      = 'sx';
-options.COMP      = 9;
-options.TIME      = [0 50 100 200 300];
+options.COMP      = 'avg';
+options.TIME      = [500];
 gene_data.a = data.EPS * 2000;
 fig = photomaton(gene_data,options);
 save_figure(gene_data,fig)
 end
 
 if 0
 %% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Options
 options.INTERP    = 1;
 options.POLARPLOT = 0;
 options.NAME      = '\phi';
 % options.NAME      = 'v_y';
 % options.NAME      = 'n_i^{NZ}';
 % options.NAME      = '\Gamma_x';
 % options.NAME      = 'n_i';
-options.PLAN      = 'xz';
+options.PLAN      = 'xy';
 % options.NAME      = 'f_e';
 % options.PLAN      = 'sx';
 options.COMP      = 'avg';
-options.TIME      = 000:300;
+options.TIME      = gene_data.Ts3D;
 gene_data.a = data.EPS * 2000;
 create_film(gene_data,options,'.gif')
 end
 
 if 0
 %% Geometry
 names = {'$g^{xx}$','$g^{xy}$','$g^{xz}$','$g^{yy}$','$g^{yz}$','$g^{zz}$',...
          '$B_0$','$\partial_x B_0$','$\partial_y B_0$','$\partial_z B_0$',...
          '$J$','$R$','$\phi$','$Z$','$\partial_R x$','$\partial_Z x$'};
 figure;
 subplot(311)
     for i = 1:6
     plot(gene_data.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on;
     end
     xlim([min(gene_data.z),max(gene_data.z)]); legend('show'); title('GENE geometry');
 
 subplot(312)
     for i = 7:10
     plot(gene_data.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on;
     end
     xlim([min(gene_data.z),max(gene_data.z)]); legend('show');
 
 subplot(313)
     for i = 11:16
     plot(gene_data.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on;
     end
     xlim([min(gene_data.z),max(gene_data.z)]); legend('show');
 
 end
 
 if 0
 %% Show f_i(vpar,mu)
 options.times   = 200:300;
 options.specie  = 'i';
 options.PLT_FCT = 'contour';
 options.folder  = folder;
 options.iz      = 9;
 options.FIELD   = '<f_>';
 options.ONED    = 0;
 % options.FIELD   = 'Q_es';
 plot_fa_gene(options);
 end
 
 if 0
 %% Mode evolution
 options.NORMALIZED = 0;
 options.K2PLOT = 1;
-options.TIME   = 100:200;
+options.TIME   = 100:700;
 options.NMA    = 1;
 options.NMODES = 15;
-options.iz     = 9;
+options.iz     = 'avg';
 fig = mode_growth_meter(gene_data,options);
 save_figure(gene_data,fig)
 end
diff --git a/wk/header_3D_results.m b/wk/header_3D_results.m
index 05a0724..9b790c1 100644
--- a/wk/header_3D_results.m
+++ b/wk/header_3D_results.m
@@ -1,17 +1,18 @@
 % Directory of the code "mypathtoHeLaZ/HeLaZ/"
 helazdir = '/home/ahoffman/HeLaZ/';
 % Directory of the simulation (from results)
 % if 1% Local results
 outfile ='';
 outfile ='';
 outfile ='';
-outfile ='shearless_cyclone/128x128x16xdmax_6_L_120_CBC_1.0';
+outfile ='shearless_cyclone/128x128x16x5x3_start';
+% outfile ='shearless_cyclone/128x128x16xdmax_6_L_120_CBC_1.0';
 % outfile ='shearless_cyclone/128x128x16xdmax_L_120_CBC_1.0';
 % outfile ='shearless_cyclone/180x180x20x4x2_L_120_CBC_0.8_to_1.0';
 % outfile ='quick_run/CLOS_1_64x64_5x3_L_120_kN_2.0_kT_0.5_nu_1e-01_SGGK';
 % outfile ='pedestal/64x64x16x2x1_L_300_LnT_20_nu_0.1';
 % outfile ='quick_run/32x32x16_5x3_L_300_q0_2.5_e_0.18_kN_20_kT_20_nu_1e-01_DGGK';
 % outfile ='shearless_cyclone/128x128x16x8x4_L_120_CTC_1.0/';
 % outfile ='shearless_cyclone/180x180x20x4x2_L_120_CBC_0.8_to_1.0/';
 % outfile ='pedestal/128x128x16x4x2_L_120_LnT_40_nuDG_0.1';
 run analysis_HeLaZ