diff --git a/matlab/compile_results.m b/matlab/compile_results.m
index 4a0324e22ab567a559b55ed4f82cf2b8a44d7e48..da22009752d189ff4e831874d2b8200bc88cf276 100644
--- a/matlab/compile_results.m
+++ b/matlab/compile_results.m
@@ -253,7 +253,7 @@ else
     KPERP2 = KX.^2+KY.^2;
     %% Add everything in output structure
     % scaling
-    DATA.scale = (1/Nx/Ny)^2;
+    DATA.scale = 1;%(1/Nx/Ny)^2;
     % Fields
     DATA.GGAMMA_RI = GGAMMAI_; DATA.PGAMMA_RI = PGAMMAI_; DATA.HFLUX_X = HFLUXI_;
     DATA.Nipj = Nipj_; DATA.Ni00 = Ni00_; DATA.DENS_I = DENS_I_; DATA.TEMP_I = TEMP_I_;
diff --git a/matlab/compute/compute_polarisation.m b/matlab/compute/compute_polarisation.m
index c1ac9f5bb18c900edf335fb029cbca1409fbfdfb..fc8e1e47feede1670a57b47d3827ac3467139955 100644
--- a/matlab/compute/compute_polarisation.m
+++ b/matlab/compute/compute_polarisation.m
@@ -14,7 +14,7 @@ KP  = sqrt(KX.^2+KY.^2);
 
 GAMMA2_ = 0.*KN_;
    
-for in = 1:Ji+1
+for in = 1:1
     for iz = 1:data.Nz
         GAMMA2_(:,:,iz) = GAMMA2_(:,:,iz) + kernel(in-1,KP*sqrt(2)).^2;
     end
diff --git a/matlab/evolve_tracers.m b/matlab/evolve_tracers.m
index 61256c51ccbf9a3c5f78571d60f487d3ddfe3c8a..c26162876a22e0261473326e88796e039eda7f84 100644
--- a/matlab/evolve_tracers.m
+++ b/matlab/evolve_tracers.m
@@ -1,17 +1,18 @@
 % Options
 SHOW_FILM = 0;
-U_TIME   = 50; % >0 for frozen velocity at a given time, -1 for evolving field
-Evolve_U = 1; % 0 for frozen velocity at a given time, 1 for evolving field
-Tfin   = 50;
-dt_    = 0.05;
-Nstep  = ceil(Tfin/dt_);
+INIT     = 'lin';   % lin (for a line)/ round (for a small round)/ gauss for random
+U_TIME   = 700;     % >0 for frozen velocity at a given time, -1 for evolving field
+Evolve_U = 1;       % 0 for frozen velocity at a given time, 1 for evolving field
+Tfin     = 10;
+dt_      = 0.0005;
+Nstep    = ceil(Tfin/dt_);
 % Init tracers
-Np      = 100; %number of tracers
+Np      = 200; %number of tracers
 % color = tcolors;
 color = jet(Np);
 tcolors = distinguishable_colors(Np); %Their colors
 
-Na = 10000; %length of trace
+Na = 1000000; %length of trace
 
 Traj_x = zeros(Np,Nstep);
 Traj_y = zeros(Np,Nstep);
@@ -21,7 +22,6 @@ Disp_y = zeros(Np,Nstep);
 xmax = max(data.x); xmin = min(data.x);
 ymax = max(data.y); ymin = min(data.y);
 
-INIT = 'round';
 switch INIT
     case 'lin'
         % Evenly distributed initial positions
diff --git a/matlab/plot/spectrum_1D.m b/matlab/plot/spectrum_1D.m
index 6c644301d5f8f83b9972561ea33c2cb9b820ff49..8c9c90b2c5a2c76518dc2a77e12d108ea0b0d5bd 100644
--- a/matlab/plot/spectrum_1D.m
+++ b/matlab/plot/spectrum_1D.m
@@ -5,6 +5,14 @@ options.COMP      = 'avg';
 options.INTERP    = 0;
 options.POLARPLOT = 0;
 options.AXISEQUAL = 1;
+
+switch options.COMPT
+    case 'avg'
+        [~,it0] = min(abs(data.Ts3D-options.TIME(1)));
+        [~,it1] = min(abs(data.Ts3D-options.TIME(end)));
+        options.TIME = data.Ts3D(it0:it1);
+end
+
 toplot = process_field(data,options);
 t = data.Ts3D; frames = toplot.FRAMES;
 
@@ -13,29 +21,40 @@ colors = jet(numel(frames));
 switch options.NAME
     case '\Gamma_x'
     FIGURE.fig = figure; FIGURE.FIGNAME = ['transp_spectrum_',data.PARAMS]; 
-    yname = '$\sum_{k_y}|\Gamma_k|$'; fieldname = 'transport';
+    fieldname = 'transport';
     case '\phi'
     FIGURE.fig = figure; FIGURE.FIGNAME = ['phi_spectrum_',data.PARAMS]; 
-    yname = '$\sum_{k_y}|\phi|$'; fieldname = 'ES pot.';
+    fieldname = 'ES pot.';
     case 'n_i'
     FIGURE.fig = figure; FIGURE.FIGNAME = ['ni_spectrum_',data.PARAMS]; 
-    yname = '$\sum_{k_y}|n_i|$'; fieldname = 'ion dens.';
+    yname = 'n_i'; fieldname = 'ion dens.';
     case 'n_e'
     FIGURE.fig = figure; FIGURE.FIGNAME = ['ne_spectrum_',data.PARAMS]; 
-    yname = '$\sum_{k_y}|n_e|$'; fieldname = 'elec. dens.';
+    fieldname = 'elec. dens.';
 end
 
 PLOT2D = 0;
 switch options.COMPXY
     case 'avg'
-        compx = @(x) mean(x,2);
-        compy = @(x) mean(x,1);
+        compx  = @(x) mean(x,2);
+        compy  = @(x) mean(x,1);
+        ynamecx = ['$\sum_{k_x}',options.NAME,'$'];
+        ynamecy = ['$\sum_{k_y}',options.NAME,'$'];
     case 'sum'
         compx = @(x) sum(x,2);
         compy = @(x) sum(x,1);
+        ynamecx = ['$\sum_{k_x}',options.NAME,'$'];
+        ynamecy = ['$\sum_{k_y}',options.NAME,'$'];
     case 'max'
         compx = @(x) max(x,2);
         compy = @(x) max(x,1);
+        ynamecx = ['$\max_{k_x}',options.NAME,'$'];
+        ynamecy = ['$\max_{k_y}',options.NAME,'$'];
+    case 'zero'
+        compx = @(x) x(:,data.Nx/2);
+        compy = @(x) x(1,:);
+        ynamecx = ['$',options.NAME,'(k_x=',num2str(data.kx(1)),')$'];
+        ynamecy = ['$',options.NAME,'(k_y=',num2str(data.ky(1)),')$'];
     otherwise
         compx =  @(x) x(:,:);
         compy =  @(x) x(:,:);
@@ -85,8 +104,9 @@ if ~PLOT2D
     end
     end
     grid on
-    title(['HeLaZ $k_x$ ',fieldname,' spectrum']); legend('show','Location','eastoutside')
-    xlabel(xname); ylabel(yname)
+    title(['GM $k_x$ ',fieldname]); 
+%     legend('show','Location','eastoutside')
+    xlabel(xname); ylabel(ynamecx)
 
     subplot(1,2,2)
 
@@ -124,9 +144,9 @@ if ~PLOT2D
     end
     end
     grid on
-    %     title('HeLaZ $k_y$ transport spectrum'); 
-    legend('show','Location','eastoutside');
-    xlabel(xname); ylabel(yname)
+    title(['GM $k_y$ ',fieldname]); 
+%     legend('show','Location','eastoutside');
+    xlabel(xname); ylabel(ynamecy)
 else
 %     it0 = toplot.FRAMES(1); it1 = toplot.FRAMES(end);
     Gk    = mean(abs(toplot.FIELD(:,:,:)),3);
diff --git a/wk/analysis_HeLaZ.m b/wk/analysis_HeLaZ.m
index 2c47363b5b599aced6933fd97b7f0da944b379b6..c80985ead6c1d7f9ef4fe6854552f9324df59d22 100644
--- a/wk/analysis_HeLaZ.m
+++ b/wk/analysis_HeLaZ.m
@@ -11,7 +11,7 @@ system(['mkdir -p ',LOCALDIR]);
 CMD = ['rsync ', LOCALDIR,'outputs* ',MISCDIR]; disp(CMD);
 system(CMD);
 % Load outputs from jobnummin up to jobnummax
-JOBNUMMIN = 00; JOBNUMMAX = 00;
+JOBNUMMIN = 00; JOBNUMMAX = 20;
 data = compile_results(MISCDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
 data.localdir = LOCALDIR;
 data.FIGDIR   = LOCALDIR;
@@ -23,7 +23,7 @@ FMT = '.fig';
 
 if 1
 %% Space time diagramm (fig 11 Ivanov 2020)
-options.TAVG_0   = 0.7*data.Ts3D(end); data.scale = (1/data.Nx/data.Ny)^2;
+options.TAVG_0   = 0.7*data.Ts3D(end); data.scale = 1;%(1/data.Nx/data.Ny)^2;
 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)
@@ -45,18 +45,18 @@ if 0
 options.INTERP    = 1;
 options.POLARPLOT = 0;
 % options.NAME      = '\phi';
-options.NAME      = 'N_i^{00}';
-% options.NAME      = 'v_y';
+% options.NAME      = 'N_i^{00}';
+% options.NAME      = 'v_x';
 % options.NAME      = 'n_i^{NZ}';
 % options.NAME      = '\Gamma_x';
-% options.NAME      = 'n_e';
-options.PLAN      = 'kxky';
+options.NAME      = 'n_i';
+options.PLAN      = 'xy';
 % options.NAME      = 'f_i';
 % options.PLAN      = 'sx';
 options.COMP      = 1;
 % options.TIME      = data.Ts5D(end-30:end);
-options.TIME      =  data.Ts3D(50:2:end);
-% options.TIME      = [350:600];
+% options.TIME      =  data.Ts3D(20:2:end);
+options.TIME      = [700:1:750];
 data.EPS          = 0.1;
 data.a = data.EPS * 2000;
 create_film(data,options,'.gif')
@@ -69,7 +69,7 @@ options.INTERP    = 1;
 options.POLARPLOT = 0;
 options.AXISEQUAL = 1;
 % options.NAME      = '\phi';
-options.NAME      = 'n_e';
+% options.NAME      = 'n_e';
 % options.NAME      = 'N_i^{00}';
 % options.NAME      = 'T_i';
 % options.NAME      = '\Gamma_x';
@@ -78,7 +78,7 @@ options.PLAN      = 'xy';
 % options.NAME      'f_i';
 % options.PLAN      = 'sx';
 options.COMP      = 'avg';
-options.TIME      = [600];
+options.TIME      = [200 800 1500];
 data.a = data.EPS * 2e3;
 fig = photomaton(data,options);
 % save_figure(data,fig)
@@ -131,16 +131,16 @@ end
 
 if 0
 %% Time averaged spectrum
-options.TIME   = 1100:10:1400;
+options.TIME   = [500 700];
 options.NORM   =1;
 % options.NAME   = '\phi';
 % options.NAME      = 'N_i^{00}';
-options.NAME      ='\Gamma_x';
+options.NAME   ='\Gamma_x';
 options.PLAN   = 'kxky';
 options.COMPZ  = 'avg';
 options.OK     = 0;
-options.COMPXY = 'avg';
-options.COMPT  = 0;
+options.COMPXY = 'avg'; % avg/sum/max/zero/ 2D plot otherwise
+options.COMPT  = 'avg';
 options.PLOT   = 'semilogy';
 fig = spectrum_1D(data,options);
 % save_figure(data,fig,'.png')
diff --git a/wk/analysis_gene.m b/wk/analysis_gene.m
index 815a21ade423539ba877a649f5edb8875a2324a5..bd2e6202a2d2f3c4b8e943292730372b8c4400db 100644
--- a/wk/analysis_gene.m
+++ b/wk/analysis_gene.m
@@ -10,8 +10,8 @@
 % folder = '/misc/gene_results/HP_fig_2b_mu_5e-2/';
 % folder = '/misc/gene_results/HP_fig_2c_mu_5e-2/';
 % folder = '/misc/gene_results/LD_zpinch_1.6/';
-folder = '/misc/gene_results/ZP_HP_kn_1.6_nuv_3.2/';
-% folder = '/misc/gene_results/ZP_HP_kn_2.5/';
+% folder = '/misc/gene_results/ZP_HP_kn_1.6_nuv_3.2/';
+folder = '/misc/gene_results/ZP_kn_2.5_large_box/';
 gene_data = load_gene_data(folder);
 gene_data = invert_kxky_to_kykx_gene_results(gene_data);
 if 1
@@ -39,15 +39,15 @@ options.INTERP    = 1;
 options.POLARPLOT = 0;
 options.AXISEQUAL = 1;
 % options.NAME      = 'Q_x';
-options.NAME      = '\phi';
-% options.NAME      = 'n_i';
+% options.NAME      = '\phi';
+options.NAME      = 'n_i';
 % options.NAME      = '\Gamma_x';
 % options.NAME      = 'k^2n_e';
 options.PLAN      = 'xy';
 % options.NAME      ='f_e';
 % options.PLAN      = 'sx';
 options.COMP      = 'avg';
-options.TIME      = [600];
+options.TIME      = [1:10];
 gene_data.a = data.EPS * 2000;
 fig = photomaton(gene_data,options);
 save_figure(gene_data,fig,'.png')
@@ -58,12 +58,12 @@ if 0
 % Options
 options.INTERP    = 1;
 options.POLARPLOT = 0;
-options.NAME      = '\phi';
+% 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      = 'n_i';
+options.PLAN      = 'kxky';
 % options.NAME      = 'f_e';
 % options.PLAN      = 'sx';
 options.COMP      = 'avg';
@@ -115,16 +115,16 @@ if 0
 %% Time averaged spectrum
 options.TIME   = 300:600;
 options.NORM   =1;
-options.NAME   = '\phi';
+% options.NAME   = '\phi';
 % options.NAME      = 'n_i';
-% options.NAME      ='\Gamma_x';
+options.NAME      ='\Gamma_x';
 options.PLAN   = 'kxky';
 options.COMPZ  = 'avg';
 options.OK     = 0;
-options.COMPXY = 'avg';
+options.COMPXY = 'zero';
 options.COMPT  = 'avg';
 options.PLOT   = 'semilogy';
-fig = spectrum_1D(data,options);
+fig = spectrum_1D(gene_data,options);
 % save_figure(data,fig)
 end
 
diff --git a/wk/header_3D_results.m b/wk/header_3D_results.m
index 3b3353860062c7a1a6631a065bed43af52cbe743..ef37d89eedfc97ea3018adddb6ad288d24b287b2 100644
--- a/wk/header_3D_results.m
+++ b/wk/header_3D_results.m
@@ -28,7 +28,10 @@ helazdir = '/home/ahoffman/HeLaZ/';
 % outfile = 'shearless_cyclone/64x32x16x5x3_CBC_CO/64x32x16x5x3_CBC_LRGK';
 %% ZPINCH
 % outfile ='Zpinch_rerun/Kn_2.5_200x48x5x3';
-outfile ='Zpinch_rerun/Kn_2.0_200x48x9x5_large_box';
+% outfile ='Zpinch_rerun/Kn_2.5_256x128x5x3_new_NL';
+outfile ='Zpinch_rerun/Kn_2.5_256x64x5x3';
+% outfile ='Zpinch_rerun/Kn_2.0_200x48x9x5_large_box';
+% outfile ='Zpinch_rerun/Kn_2.0_256x64x9x5_Lx_240_Ly_120';
 % outfile ='Zpinch_rerun/Kn_1.6_256x128x7x4';
 % outfile ='Zpinch_rerun/Kn_1.6_200x48x11x6';
 % outfile ='Zpinch_rerun/Kn_1.6_200x48x21x11';
diff --git a/wk/quick_run.m b/wk/quick_run.m
index 687ef2dcc248b4d9e567265f3fe0c05408843e01..63c31a65815cddde71ad535fe0b9150512981cee 100644
--- a/wk/quick_run.m
+++ b/wk/quick_run.m
@@ -13,48 +13,48 @@ EXECNAME = 'helaz3';
 CLUSTER.TIME  = '99:00:00'; % allocation time hh:mm:ss
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% PHYSICAL PARAMETERS
-NU      = 0.0;   % Collision frequency
+NU      = 0.5;   % Collision frequency
 TAU     = 1.0;    % e/i temperature ratio
-K_N     = 1.0;%2.0;   % Density gradient drive
+K_N     = 6.0;%2.0;   % Density gradient drive
 K_T     = 0;%0.25*K_N;   % Temperature '''
 K_E     = 0.0;   % Electrostat '''
 % SIGMA_E = 0.05196152422706632;   % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
 SIGMA_E = 0.0233380;   % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
 KIN_E   = 1;     % 1: kinetic electrons, 2: adiabatic electrons
 %% GRID PARAMETERS
-PMAXE   = 30;     % Hermite basis size of electrons
-JMAXE   = 15;     % Laguerre "
-PMAXI   = 30;     % " ions
-JMAXI   = 15;     % "
-NX      = 1;    % real space x-gridpoints
-NY      = 32;     %     ''     y-gridpoints
+PMAXE   = 4;     % Hermite basis size of electrons
+JMAXE   = 2;     % Laguerre "
+PMAXI   = 4;     % " ions
+JMAXI   = 2;     % "
+NX      = 32;    % real space x-gridpoints
+NY      = 1;     %     ''     y-gridpoints
 LX      = 100;   % Size of the squared frequency domain
-LY      = 30;     % Size of the squared frequency domain
-NZ      = 1;     % number of perpendicular planes (parallel grid)
+LY      = 2*pi/0.5;     % Size of the squared frequency domain
+NZ      = 16;     % number of perpendicular planes (parallel grid)
 NPOL    = 1;
 SG      = 0;     % Staggered z grids option
 %% GEOMETRY
 % GEOMETRY= 'Z-pinch'; % Z-pinch overwrites q0, shear and eps
 GEOMETRY= 's-alpha';
-Q0      = 1.0;    % safety factor
-SHEAR   = 0.0;    % magnetic shear (Not implemented yet)
-EPS     = 0.0;    % inverse aspect ratio
+Q0      = 1.4;    % safety factor
+SHEAR   = 0.796;    % magnetic shear (Not implemented yet)
+EPS     = 0.18;    % inverse aspect ratio
 %% TIME PARMETERS
-TMAX    = 500;  % Maximal time unit
-DT      = 5e-3;   % Time step
+TMAX    = 50;  % Maximal time unit
+DT      = 1e-3;   % Time step
 SPS0D   = 1;      % Sampling per time unit for 2D arrays
 SPS2D   = 0;      % Sampling per time unit for 2D arrays
 SPS3D   = 1;      % Sampling per time unit for 2D arrays
-SPS5D   = 1/2;    % Sampling per time unit for 5D arrays
+SPS5D   = 1/5;    % Sampling per time unit for 5D arrays
 SPSCP   = 0;    % Sampling per time unit for checkpoints
 JOB2LOAD= -1;
 %% OPTIONS
-SIMID   = 'dbg';  % Name of the simulation
+SIMID   = 'shear_testcase_Pan22_linear';  % Name of the simulation
 LINEARITY = 'linear';   % activate non-linearity (is cancelled if KXEQ0 = 1)
 % Collision operator
 % (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
-CO      = 'DG';
-GKCO    = 1; % gyrokinetic operator
+CO      = 'SG';
+GKCO    = 0; % gyrokinetic operator
 ABCO    = 1; % interspecies collisions
 INIT_ZF = 0; ZF_AMP = 0.0;
 CLOS    = 0;   % Closure model (0: =0 truncation, 1: gyrofluid closure (p+2j<=Pmax))s
@@ -77,7 +77,7 @@ INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
 MU_X    = MU;     %
 MU_Y    = MU;     %
 N_HD    = 4;
-MU_Z    = 0.0;     %
+MU_Z    = 0.05;     %
 MU_P    = 0.0;     %
 MU_J    = 0.0;     %
 LAMBDAD = 0.0;
@@ -92,8 +92,8 @@ system(['rm fort*.90']);
 % Run linear simulation
 if RUN
 %     system(['cd ../results/',SIMID,'/',PARAMS,'/; time mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk'])
-    system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',HELAZDIR,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk'])
-%     system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 1 4 0; cd ../../../wk'])
+%     system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',HELAZDIR,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk'])
+    system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ',HELAZDIR,'bin/',EXECNAME,' 1 1 2 0; cd ../../../wk'])
 %     system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk'])
 end