diff --git a/matlab/compute/compute_fa_2D.m b/matlab/compute/compute_fa_2D.m
index 50b3593..c66ebba 100644
--- a/matlab/compute/compute_fa_2D.m
+++ b/matlab/compute/compute_fa_2D.m
@@ -1,107 +1,107 @@
 function [SS,XX,FF] = compute_fa_2D(data, options)
 %% Compute the dispersion function from the moment hierarchi decomp.
 % Normalized Hermite
 Hp = @(p,s) polyval(HermitePoly(p),s)./sqrt(2.^p.*factorial(p));
 % Hp = @(p,s) hermiteH(p,s)./sqrt(2.^p.*factorial(p));
 % Laguerre
 Lj = @(j,x) polyval(LaguerrePoly(j),x);
 % Maxwellian factor
 FaM = @(s,x) exp(-s.^2-x);
 
 %meshgrid for efficient evaluation
 [SS, XX] = meshgrid(options.SPAR, options.XPERP);
 
 [~, ikx0] = min(abs(data.kx));
 [~, iky0] = min(abs(data.ky));
 kx_ = data.kx; ky_ = data.ky;
 
 switch options.SPECIE
     case 'e'
         Napj_     = data.Nepj;
         parray    = double(data.Pe);
         jarray    = double(data.Je);
     case 'i'
         Napj_     = data.Nipj;
         parray    = double(data.Pi);
         jarray    = double(data.Ji);
 end
 
 switch options.Z
     case 'avg'
         Napj_     = mean(Napj_,5);
         phi_      = mean(data.PHI,3);
     otherwise
         iz        = options.Z; 
         Napj_     = Napj_(:,:,:,:,iz,:);
         phi_      = data.PHI(:,:,iz);
 end
 Napj_ = squeeze(Napj_);
 
 frames = options.T;
 for it = 1:numel(options.T)
     [~,frames(it)] = min(abs(options.T(it)-data.Ts5D)); 
 end
 frames = unique(frames);
 
 Napj_     = mean(Napj_(:,:,:,:,frames),5);
 
 Napj_ = squeeze(Napj_);
 
 
 Np = numel(parray); Nj = numel(jarray);
 
 if options.non_adiab
     for ij_ = 1:Nj
         for ikx = 1:data.Nkx
             for iky = 1:data.Nky    
                 kp_ = sqrt(kx_(ikx)^2 + ky_(iky)^2);
-                Napj_(1,ij_,ikx,iky) = Napj_(1,ij_,ikx,iky) + kernel(ij_,kp_)*phi_(ikx,iky);
+                Napj_(1,ij_,iky,ikx) = Napj_(1,ij_,iky,ikx) + kernel(ij_,kp_)*phi_(iky,ikx);
             end
         end
     end
 end
 
 if options.RMS
-    FF = zeros(data.Nkx,data.Nky,numel(options.XPERP),numel(options.SPAR));
+    FF = zeros(data.Nky,data.Nkx,numel(options.XPERP),numel(options.SPAR));
     FAM = FaM(SS,XX);
     for ip_ = 1:Np
         p_ = parray(ip_);
         HH = Hp(p_,SS);
         for ij_ = 1:Nj
             j_ = jarray(ij_);
             LL = Lj(j_,XX);
             HLF = HH.*LL.*FAM;
             for ikx = 1:data.Nkx
                 for iky = 1:data.Nky
-                    FF(ikx,iky,:,:) = squeeze(FF(ikx,iky,:,:)) + Napj_(ip_,ij_,ikx,iky)*HLF;
+                    FF(iky,ikx,:,:) = squeeze(FF(iky,ikx,:,:)) + Napj_(ip_,ij_,iky,ikx)*HLF;
                 end
             end
        end
     end
 else
     FF = zeros(numel(options.XPERP),numel(options.SPAR));
     FAM = FaM(SS,XX);
     for ip_ = 1:Np
         p_ = parray(ip_);
         HH = Hp(p_,SS);
         for ij_ = 1:Nj
             j_ = jarray(ij_);
             LL = Lj(j_,XX);
             FF = FF + squeeze(Napj_(ip_,ij_,ikx0,iky0))*HH.*LL.*FAM;
         end
     end
 end
 FF = (FF.*conj(FF)); %|f_a|^2
 % FF = abs(FF); %|f_a|
 if options.RMS
 %     FF = squeeze(mean(mean(sqrt(FF),1),2)); %sqrt(<|f_a|^2>kx,ky)
     FF = sqrt(squeeze(mean(mean(FF,1),2))); %<|f_a|>kx,ky
 else
     FF = sqrt(squeeze(FF)); %sqrt(<|f_a|>x,y)
 end
 
 FF = FF./max(max(FF));
 FF = FF';
 % FF = sqrt(FF);
 % FF = FF';
 end
\ No newline at end of file
diff --git a/wk/analysis_3D.m b/wk/analysis_3D.m
index 0f2cdef..34dbcad 100644
--- a/wk/analysis_3D.m
+++ b/wk/analysis_3D.m
@@ -1,195 +1,198 @@
 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 = 00; JOBNUMMAX = 10;
 data = compile_results(MISCDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
 
 
 %% PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 default_plots_options
 disp('Plots')
 FMT = '.fig';
 
 if 1
 %% Space time diagramm (fig 11 Ivanov 2020)
 options.TAVG_0   = 0.8*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      = 'xy';
 % options.NAME      = 'f_e';
 % options.PLAN      = 'sx';
 options.COMP      = 'avg';
 % options.TIME      = dat.Ts5D;
 options.TIME      = 00:1:250;
 data.EPS          = 0.1;
 data.a = data.EPS * 2000;
 create_film(data,options,'.gif')
 end
 
 if 0
 %% 2D snapshots
 % Options
 options.INTERP    = 0;
 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      = 'xy';
 % options.NAME      = 'f_e';
 % options.PLAN      = 'sx';
 options.COMP      = 1;
 options.TIME      = [150];
 data.a = data.EPS * 1000;
 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;
 options.TIME      = [100 200];
 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      = vp';
-% options.XPERP     = mu';
-options.Z         = 1;
+% 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.Z         = 'avg';
 options.T         = 200;
-options.CTR       = 1;
+options.PLT_FCT   = 'pcolor';
 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.NORMALIZED = 0;
 fig = show_moments_spectrum(data,options);
 save_figure(data,fig)
 end
 
 if 0
 %% Time averaged spectrum
 options.TIME   = 1000:5000;
 options.NORM   = 0;
 options.NAME   = '\phi';
 options.NAME      = 'n_i';
 options.NAME      ='\Gamma_x';
 options.PLAN   = 'kxky';
 options.COMPZ  = 'avg';
 options.OK     = 0;
 options.COMPXY = 0;
 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 = 1;
 options.K2PLOT = 1;
 options.TIME   = 5:1:15;
 options.NMA    = 1;
 options.NMODES = 5;
 options.iz     = 1;
 fig = mode_growth_meter(data,options);
 save_figure(data,fig)
 end
 
 if 0
 %% ZF caracteristics (space time diagrams)
 TAVG_0 = 200; TAVG_1 = 3000; % 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/gene_analysis_3D.m b/wk/gene_analysis_3D.m
index ed4d356..7479f55 100644
--- a/wk/gene_analysis_3D.m
+++ b/wk/gene_analysis_3D.m
@@ -1,181 +1,181 @@
 % 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.2/';
+folder = '/misc/gene_results/shearless_cyclone/s_alpha_output_1.0/';
 % 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 = rotate_c_plane_nxnky_to_nkxny(gene_data);
 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;
 fig = plot_radial_transport_and_spacetime(gene_data,options);
 % save_figure(data,fig)
 end
 
 if 0
 %% 2D snapshots
 % Options
 options.INTERP    = 1;
 options.POLARPLOT = 0;
 options.AXISEQUAL = 1;
 options.NAME      = '\phi';
 % options.NAME      = 'n_i';
 % options.NAME      = 'T_i';
 % options.NAME      = '\Gamma_x';
 % options.NAME      = 'k^2n_e';
 options.PLAN      = 'xz';
 % options.NAME      ='f_e';
 % options.PLAN      = 'sx';
 options.COMP      = 'avg';
 options.TIME      = [100 300 900];
 gene_data.a = data.EPS * 2000;
 fig = photomaton(gene_data,options);
 save_figure(gene_data,fig)
 end
 
 if 0
 %% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Options
 options.INTERP    = 0;
 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      = 'xy';
 % options.NAME      = 'f_e';
 % options.PLAN      = 'sx';
 options.COMP      = 'avg';
 options.TIME      = 000:200;
 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
 % folder = '/misc/gene_results/HP_fig_2b_mu_5e-2/';
 folder = '/misc/gene_results/shearless_cyclone/s_alpha_output_1.0/';
 % folder = '/misc/gene_results/cyclone/s_alpha_output_1.0/';
 gene_data = load_gene_data(folder);
 gene_data = rotate_c_plane_nxnky_to_nkxny(gene_data);
 if 1
 %% Space time diagramm (fig 11 Ivanov 2020)
 TAVG_0 = 0.8*gene_data.Ts3D(end); TAVG_1 = gene_data.Ts3D(end); % Averaging times duration
-% chose your field to plot in spacetime diag (uzf,szf,Gx)
-field = 'phi';
-compz = 'avg';
-nmvm  = 1;
-fig = plot_radial_transport_and_spacetime(gene_data,TAVG_0,TAVG_1,field,nmvm,compz);
+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(gene_data,options);
 % save_figure(data,fig)
 end
 
 if 0
 %% 2D snapshots
 % Options
 options.INTERP    = 1;
 options.POLARPLOT = 0;
 options.AXISEQUAL = 1;
 options.NAME      = '\phi';
 % options.NAME      = 'n_i';
 % options.NAME      = 'T_i';
 % options.NAME      = '\Gamma_x';
 % options.NAME      = 'k^2n_e';
 options.PLAN      = 'xz';
 % options.NAME      ='f_e';
 % options.PLAN      = 'sx';
 options.COMP      = 'avg';
 options.TIME      = [100 300 900];
 gene_data.a = data.EPS * 2000;
 fig = photomaton(gene_data,options);
 save_figure(gene_data,fig)
 end
 
 if 0
 %% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Options
 options.INTERP    = 0;
 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      = 'xy';
 % options.NAME      = 'f_e';
 % options.PLAN      = 'sx';
 options.COMP      = 'avg';
 options.TIME      = 000:170;
 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:600;
 options.specie  = 'i';
 options.PLT_FCT = 'pcolor';
 options.folder  = folder;
 options.Z       = 'avg';
 options.FIELD   = '<f_>';
 options.ONED    = 0;
 % options.FIELD   = 'Q_es';
 plot_fa_gene(options);
 end