diff --git a/matlab/extract_fig_data.m b/matlab/extract_fig_data.m
index e64b370b5414d03f9f3bc4a9d21da581b81f0456..e7a98b9d7301ab9939cfabe60f3d35ac0155008a 100644
--- a/matlab/extract_fig_data.m
+++ b/matlab/extract_fig_data.m
@@ -7,7 +7,7 @@ for i = 1:numel(dataObjs)
     X_ = [X_ dataObjs(i).XData];
     Y_ = [Y_ dataObjs(i).YData];
 end
-n0 = 250;
+n0 = 1;
 n1 = numel(X_);
 figure;
 mvm = @(x) movmean(x,1);
@@ -16,7 +16,7 @@ shift = X_(n0);
 % plot(X_(n0:end),Y_(n0:end));
 plot(mvm(X_(n0:n1)-shift),mvm(Y_(n0:n1))); hold on;
 
-t0 = ceil(numel(X_)*0.5); t1 = numel(X_);
+t0 = ceil(numel(X_)*0.2); t1 = numel(X_);
 avg= mean(Y_(t0:t1)); dev = std(Y_(t0:t1));
   disp(['AVG =',sprintf('%2.2f',avg),'+-',sprintf('%2.2f',dev)]);
 % 
diff --git a/matlab/load/quick_gene_load.m b/matlab/load/quick_gene_load.m
index 3ac095c94f09c8abae00640542a841853413bdd0..8fb67265b2114225e6482cc0fa57574bf97c1c93 100644
--- a/matlab/load/quick_gene_load.m
+++ b/matlab/load/quick_gene_load.m
@@ -41,13 +41,17 @@
 % fname ='GENE_LIN_Kn_1.6_KT_0.4_nu_0_32x16.txt';
 % fname ='GENE_LIN_Kn_2.5_KT_0.625_nu_0_32x16.txt';
 path = '/home/ahoffman/gene/linear_CBC_results/';
+fname = 'CBC_100_20x1x32x30x14_Lv_3_Lw_12_circ.txt';
 % fname = 'CBC_100_20x1x32x32x12_Lv_3_Lw_12.txt';
 % fname = 'CBC_KT_4_20x1x32x32x12_Lv_3_Lw_12.txt';
 % fname = 'CBC_KT_4_20x1x32x64x24_Lv_6_Lw_24.txt';
 % fname = 'CBC_KT_5.3_20x1x32x32x12_Lv_3_Lw_12.txt';
 % fname = 'CBC_KT_5.3_32x1x48x40x16_Lv_3_Lw_12.txt';
-fname = 'CBC_ky_0.3_20x1x32x32x12_Lv_3_Lw_12.txt';
+% fname = 'CBC_ky_0.3_20x1x32x32x12_Lv_3_Lw_12.txt';
 % fname = 'CBC_ky_0.3_20x1x32x32x12_Lv_3_Lw_12_nuv_1e-3.txt';
+% fname = 'CBC_KT_11_20x1x16x24x10_Lv_3_Lw_12.txt';
+% fname = 'CBC_KT_11_20x1x32x30x14_Lv_3_Lw_12.txt';
+% fname = 'CBC_ky_0.3_20x1x16x24x10_Lv_3_Lw_12_nuv_1e-3.txt';
 data_ = load([path,fname]);
 
 figure
diff --git a/matlab/plot/plot_ballooning.m b/matlab/plot/plot_ballooning.m
index f4deff2ac21791390eb6adfde3bbc1e3212ac67a..a564aa6590c484480de0a9da90036939e4f23674 100644
--- a/matlab/plot/plot_ballooning.m
+++ b/matlab/plot/plot_ballooning.m
@@ -8,14 +8,13 @@ function [FIG] = plot_ballooning(data,options)
     phi_real=real(data.PHI(:,:,:,it1));
     phi_imag=imag(data.PHI(:,:,:,it1));
     % Apply baollooning tranform
+    nexc = round(data.ky(2)*data.SHEAR*2*pi/data.kx(2));
     for iky=ikyarray
         dims = size(phi_real);
         Nkx  = dims(2);
         is   = max(1,iky-1);
-        Npi  = (Nkx-1)-2*(is-1);
-        if(Npi <= 0)
-            break
-        elseif(Npi == 1)
+        Npi  = (Nkx-1)-2*nexc*(is-1);
+        if(Npi <= 1)
             ordered_ikx = 1;
         else
             tmp_ = (Nkx-is+1):-is:(Nkx/2+2);
diff --git a/matlab/plot/plot_radial_transport_and_spacetime.m b/matlab/plot/plot_radial_transport_and_spacetime.m
index 7ad93c0b6bbb408d6f7420c8312e990407d0b3af..de18d2d53a5d0431135166c8b020ea3fa38c76ba 100644
--- a/matlab/plot/plot_radial_transport_and_spacetime.m
+++ b/matlab/plot/plot_radial_transport_and_spacetime.m
@@ -17,7 +17,17 @@ function [FIGURE] = plot_radial_transport_and_spacetime(DATA, OPTIONS)
     ikzf = min([ikzf,DATA.Nky]);
     Ns3D = numel(DATA.Ts3D);
     [KX, KY] = meshgrid(DATA.kx, DATA.ky);
-    
+    %% error estimation
+    DT_       = (tend-tstart)/OPTIONS.NCUT;
+    Qx_ee     = zeros(1,OPTIONS.NCUT);
+    for i = 1:OPTIONS.NCUT
+        [~,its_] = min(abs(DATA.Ts0D - (tstart+(i-1)*DT_)));
+        [~,ite_] = min(abs(DATA.Ts0D - (tstart+ i   *DT_)));
+        Qx_ee(i) = mean(DATA.HFLUX_X(its_:ite_))*SCALE;
+    end
+    Qx_avg    = mean(Qx_ee);
+    Qx_err    =  std(Qx_ee);
+    disp(['Q_avg=',sprintf('%2.2f',Qx_avg),'+-',sprintf('%2.2f',Qx_err)]);
     %% computations
 
     % Compute Gamma from ifft matlab
@@ -66,7 +76,7 @@ mvm = @(x) movmean(x,OPTIONS.NMVA);
         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)]);
+            'DisplayName',['$Q_{avg}=',sprintf('%2.2f',Qx_avg),'\pm',sprintf('%2.2f',Qx_err),'$']); legend('show');
         ylabel('$Q_x$')  
         ylim([0,5*abs(Qx_infty_avg)]); 
         xlim([DATA.Ts0D(1),DATA.Ts0D(end)]);
diff --git a/matlab/setup.m b/matlab/setup.m
index a11f45461138342c4220f61a4c568aa063465f70..163ffe1b0a17be8b06806c55b7584d8459781d29 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -7,6 +7,7 @@ GRID.pmaxi = PMAXI;  % Ion Hermite moments
 GRID.jmaxi = JMAXI;  % Ion Laguerre moments
 GRID.Nx    = NX; % x grid resolution
 GRID.Lx    = LX; % x length
+GRID.Nexc  = NEXC; % to extend Lx when s>0
 GRID.Ny    = NY; % y ''
 GRID.Ly    = LY; % y ''
 GRID.Nz    = NZ;    % z resolution
@@ -126,6 +127,11 @@ if (NZ > 1)  %3D case
        geo_   = [geo_,'_s_',num2str(SHEAR)];
     end
 end
+switch GEOMETRY
+    case 'circular'
+       geo_   = [geo_,'_circ_'];
+end
+        
 % put everything together in the param character chain
 u_ = '_'; % underscore variable
 PARAMS = [res_,HLdeg_,geo_,drives_,coll_,lin_,adiabe_];
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index 90dee735fcacc54f79262fe2fa2ba42e544ac7e3..2e769aefa1e8e04b89bd7402b6992730d41d25bf 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -17,6 +17,7 @@ fprintf(fid,['  pmaxi  = ', num2str(GRID.pmaxi),'\n']);
 fprintf(fid,['  jmaxi  = ', num2str(GRID.jmaxi),'\n']);
 fprintf(fid,['  Nx     = ', num2str(GRID.Nx),'\n']);
 fprintf(fid,['  Lx     = ', num2str(GRID.Lx),'\n']);
+fprintf(fid,['  Nexc   = ', num2str(GRID.Nexc),'\n']);
 fprintf(fid,['  Ny     = ', num2str(GRID.Ny),'\n']);
 fprintf(fid,['  Ly     = ', num2str(GRID.Ly),'\n']);
 fprintf(fid,['  Nz     = ', num2str(GRID.Nz),'\n']);
diff --git a/wk/CBC_ky_PJ_scan.m b/wk/CBC_ky_PJ_scan.m
new file mode 100644
index 0000000000000000000000000000000000000000..d9bc9e0516236adba01b11d213b395df1063d353
--- /dev/null
+++ b/wk/CBC_ky_PJ_scan.m
@@ -0,0 +1,122 @@
+addpath(genpath('../matlab')) % ... add
+default_plots_options
+HELAZDIR = '/home/ahoffman/HeLaZ/';
+EXECNAME = 'helaz3';
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+KY_a = 0.1:0.1:0.8;
+g_max= KY_a*0;
+g_avg= g_max*0;
+g_std= g_max*0;
+k_max= g_max*0;
+
+CO    = 'DG'; GKCO = 0;
+NU    = 0.05;
+DT    = 2e-3;
+TMAX  = 25;
+K_T    = 6.96;
+SIMID = 'linear_CBC_circ_conv';  % Name of the simulation
+RUN   = 0;
+% figure
+% P = 12;
+for P = [12]
+J = P/2;
+
+i=1;
+for ky_ = KY_a
+    
+%Set Up parameters
+for j = 1
+    CLUSTER.TIME  = '99:00:00'; % allocation time hh:mm:ss
+    TAU     = 1.0;    % e/i temperature ratio
+    K_N = 2.22; K_Ne = K_N;
+    K_Te     = K_T;            % Temperature '''
+    SIGMA_E = 0.0233380;   % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
+    KIN_E   = 0;     % 1: kinetic electrons, 2: adiabatic electrons
+    BETA    = 0e-1;     % electron plasma beta 
+    PMAXE   = P; JMAXE   = J;
+    PMAXI   = P; JMAXI   = J;
+    NX      = 12;    % real space x-gridpoints
+    NY      = 2;     %     ''     y-gridpoints
+    LX      = 2*pi/0.1;   % Size of the squared frequency domain
+    LY      = 2*pi/ky_;
+    NZ      = 16;    % number of perpendicular planes (parallel grid)
+    NPOL    = 1; SG = 0;
+%     GEOMETRY= 's-alpha';
+    GEOMETRY= 'circular';
+    Q0      = 1.4;    % safety factor
+    SHEAR   = 0.8;    % magnetic shear (Not implemented yet)
+    EPS     = 0.18;    % inverse aspect ratio
+    SPS0D = 1; SPS2D = 0; SPS3D = 1;SPS5D= 1/5; SPSCP = 0;
+    JOB2LOAD= -1;
+    LINEARITY = 'linear';   % activate non-linearity (is cancelled if KXEQ0 = 1)
+    ABCO    = 1; % interspecies collisions
+    INIT_ZF = 0; ZF_AMP = 0.0;
+    CLOS    = 0;   % Closure model (0: =0 truncation, 1: v^Nmax closure (p+2j<=Pmax))s
+    NL_CLOS = 0;   % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS)
+    KERN    = 0;   % Kernel model (0 : GK)
+    INIT_OPT= 'phi';   % Start simulation with a noisy mom00/phi/allmom
+    W_DOUBLE = 1;
+    W_GAMMA  = 1; W_HF     = 1;
+    W_PHI    = 1; W_NA00   = 1;
+    W_DENS   = 1; W_TEMP   = 1;
+    W_NAPJ   = 1; W_SAPJ   = 0;
+    HD_CO   = 0.0;    % Hyper diffusivity cutoff ratio
+    MU      = 0.0; % Hyperdiffusivity coefficient
+    INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
+    MU_X    = MU;     %
+    MU_Y    = MU;  N_HD    = 4;
+    MU_Z    = 2.0; MU_P    = 0.0;     %
+    MU_J    = 0.0; LAMBDAD = 0.0;
+    NOISE0  = 0.0e-5; % Init noise amplitude
+    BCKGD0  = 1.0;    % Init background
+GRADB   = 1.0;CURVB   = 1.0;
+end
+%%-------------------------------------------------------------------------
+% RUN
+setup
+if RUN
+    system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 2 1 3 0; cd ../../../wk'])
+end
+
+% Load results
+filename = [SIMID,'/',PARAMS,'/'];
+LOCALDIR  = [HELAZDIR,'results/',filename,'/'];
+data = compile_results(LOCALDIR,0,0); %Compile the results from first output found to JOBNUMMAX if existing
+
+%linear growth rate (adapted for 2D zpinch and fluxtube)
+trange = [0.5 1]*data.Ts3D(end);
+nplots = 0;
+lg = compute_fluxtube_growth_rate(data,trange,nplots);
+[gmax,     kmax] = max(lg.g_ky(:,end));
+[gmaxok, kmaxok] = max(lg.g_ky(:,end)./lg.ky);
+msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,lg.ky(kmax)); disp(msg);
+msg = sprintf('gmax/k = %2.2f, kmax/k = %2.2f',gmaxok,lg.ky(kmaxok)); disp(msg);
+    
+    g_max(i) = gmax;
+    k_max(i) = kmax;
+    
+    g_avg(i) = lg.avg_g;
+    g_std(i) = lg.std_g;
+    
+    i = i + 1;
+end
+%%
+
+% plot(KT_a,max(g_max,0));
+y_ = g_avg; 
+e_ = g_std;
+
+y_ = y_.*(y_-e_>0);
+e_ = e_ .* (y_>0);
+% errorbar(KY_a,y_,e_,...
+%     'LineWidth',1.2,...
+%     'DisplayName',['(',num2str(P),',',num2str(J),')']); 
+% hold on;
+% title(['Linear CBC $K_T$ threshold $k_y=$',num2str(ky_),' (CLOS = 1)']);
+% legend('show'); xlabel('$k_y$'); ylabel('$\gamma$');
+% drawnow
+
+fig = plot_ballooning(data,options);
+
+end
+
diff --git a/wk/CBC_nu_CO_scan.m b/wk/CBC_nu_CO_scan.m
index 949847b83d6628c21e9e32a3bfd42ca1a093e36a..d14d69ba31a7db9aad1fe96ad5f9474477d7e13d 100644
--- a/wk/CBC_nu_CO_scan.m
+++ b/wk/CBC_nu_CO_scan.m
@@ -1,20 +1,21 @@
 % NU_a = [0.05 0.15 0.25 0.35 0.45];
-NU_a = [0:0.05:0.5];
+NU_a = [0:0.1:0.5];
 g_max= NU_a*0;
 g_avg= NU_a*0;
 g_std= NU_a*0;
 k_max= NU_a*0;
-CO      = 'LD';
+CO      = 'DG';
 
-K_T   = 5.3;
-DT    = 2e-3;
-TMAX  = 30;
+K_T   = 7;
+DT    = 5e-3;
+TMAX  = 20;
 ky_   = 0.3;
-SIMID = 'linear_CBC_nu_scan_kT_5.3_ky_0.3_LDGK';  % Name of the simulation
-RUN   = 0;
+SIMID = 'linear_CBC_nu_scan_kT_7_ky_0.3_DGGK';  % Name of the simulation
+% SIMID = 'linear_CBC_nu_scan_kT_11_ky_0.3_DGGK';  % Name of the simulation
+RUN   = 1;
 figure
 
-for P = [2 4 6]
+for P = [6 8 10]
 
 i=1;
 for NU = NU_a
diff --git a/wk/Dimits_2000_fig3.m b/wk/Dimits_2000_fig3.m
index 8074cdf5d3ae4f229d20caa9ca5fbf9065904f18..13501cd86b7b315a4a4fef5585ac43d6347076c8 100644
--- a/wk/Dimits_2000_fig3.m
+++ b/wk/Dimits_2000_fig3.m
@@ -1,63 +1,64 @@
 %% Heat flux Qi [R/rhos^2/cs]
 kN = 2.22;
 %-------------- GM ---------------
-%192x96x16x3x2 kymin=0.05
+%(P,J)=(2,1)
 kT_Qi_GM_32 = ...
     [...
-     9.0 1.0e+2 4.3e+1;...
-     7.0 3.6e+1 0;...
-     6.0 2.6e+1 0;...
-     5.0 1.1e+1 0;...
-     4.5 7.8e+0 0;...
+     13. 1.5e+2 1.6e+1;...%192x96x16x3x2 kymin=0.02
+     11. 5.1e+2 3.5e+2;...%192x96x16x3x2 kymin=0.02
+     9.0 9.6e+1 3.0e+1;...%192x96x16x3x2 kymin=0.05
+     7.0 5.0e+1 6.6e+0;...%192x96x16x3x2 kymin=0.05
+     6.0 3.0e+1 4.8e+0;...%192x96x16x3x2 kymin=0.05
+     5.0 1.1e+1 9.4e-1;...%192x96x16x3x2 kymin=0.05
+     4.5 9.2e+0 1.6e+0;...%192x96x16x3x2 kymin=0.05
     ];
-%128x64x16x5x3 kymin=0.05
+%(P,J)=(4,2)
 kT_Qi_GM_53 = ...
     [...
-     13. 1.3e+2 3.5e+1;...
-     11. 9.7e+1 2.2e+1;...
-     9.0 9.0e+1 2.8e+1;...
-     7.0 4.6e+1 0;...
-     5.3 1.9e+1 0;...
-     5.0 1.5e+1 0;...
-     4.5 9.7e+0 0;...9_128
+     13. 2.0e+2 1.2e+1;...%96x64x16x3x2  kymin=0.02 (large box)
+%      13. 1.1e+2 2.0e+1;...%128x64x16x5x3 kymin=0.02 (large box)
+%      13. 1.3e+2 3.5e+1;...%128x64x16x5x3 kymin=0.05
+     11. 1.2e+2 1.6e+1;...%96x64x16x3x2  kymin=0.02 (large box)
+%      11. 1.6e+2 1.8e+1;...%128x64x16x5x3 kymin=0.02
+%      11. 9.7e+1 2.2e+1;...%128x64x16x5x3 kymin=0.05
+     9.0 8.3e+1 2.2e+1;...%128x64x16x5x3 kymin=0.05
+%      9.0 7.6e+1 2.3e+1;...%256x128x16x3x2 kymin=0.05 (high res)
+     7.0 4.6e+1 2.3e+0;...%128x64x16x5x3 kymin=0.05
+     6.0 3.7e+1 6.9e+0;...%128x64x16x5x3 kymin=0.05
+     5.3 1.9e+1 2.0e+0;...%128x64x16x5x3 kymin=0.05
+     5.0 1.3e+1 3.3e+0;...%128x64x16x5x3 kymin=0.05
+     4.5 9.3e+0 1.0e+0;...%128x64x16x5x3 kymin=0.05
     ];
-%128x64x16x5x3 kymin=0.02 (large box)
-kT_Qi_GM_53_LB = ...
+%(P,J)=(12,2) or higher
+kT_Qi_GM_122 = ...
     [...
-     13. 1.1e+2 2.0e+1;...
-    ];
-%256x128x16x3x2 kymin=0.05 (high res)
-kT_Qi_GM_32_HR = ...
-    [...
-     9. 7.6e+1 2.3e+1;...
+     7.0 4.1e+1 6.6e+0;...%192x96x24x13x7 kymin=0.05
+     4.5 9.6e-1 1.5e-1;...%128x64x16x13x2 kymin=0.05
     ];
 %-------------- GENE ---------------
-%128x64x16x24x12 kymin=0.05
 kT_Qi_GENE = ...
     [...
-     13. 2.0e+2 6.6e+1;...
-     11. 3.3e+2 1.6e+2;...
-     9.0 1.1e+2 0;...
-     7.0 5.0e+1 0;...
-     5.3 2.4e+1 0;...
-     4.5 1.9e-1 0;...
-    ];
-%128x64x16x24x12 kymin=0.02 (large box)
-kT_Qi_GENE_53_LB = ...
-    [...
-     13. 2.7e+2 2.2e+1;...
-     11. 1.9e+2 1.7e+1;...
+     13. 2.7e+2 2.2e+1;...%128x64x16x24x12 kymin=0.02 (large box)
+%      13. 2.0e+2 6.6e+1;...%128x64x16x24x12 kymin=0.05
+     11. 1.9e+2 1.7e+1;...%128x64x16x24x12 kymin=0.02 (large box)
+%      11. 3.3e+2 1.6e+2;...%128x64x16x24x12 kymin=0.05
+     9.0 1.1e+2 4.2e+1;...%128x64x16x24x12 kymin=0.05
+     7.0 4.1e+1 2.1e+1;...%128x64x16x24x12 kymin=0.05
+     5.3 1.1e+1 1.8e+1;...%128x64x16x24x12 kymin=0.05
+     4.5 1.9e-1 3.0e-2;...%128x64x16x24x12 kymin=0.05
     ];
 %% Heat conductivity Xi [Ln/rhoi^2/cs] computed as Xi = Qi/kT/kN
 %init
-kT_Xi_GM_32 = kT_Qi_GM_32;
-kT_Xi_GM_53 = kT_Qi_GM_53;
-kT_Xi_GENE  = kT_Qi_GENE;
+kT_Xi_GM_32  = kT_Qi_GM_32;
+kT_Xi_GM_53  = kT_Qi_GM_53;
+kT_Xi_GM_122 = kT_Qi_GM_122;
+kT_Xi_GENE   = kT_Qi_GENE;
 %scale
 for i = 2:3
-kT_Xi_GM_32(:,i) = kT_Qi_GM_32(:,i)./kT_Qi_GM_32(:,1)./kN;
-kT_Xi_GM_53(:,i) = kT_Qi_GM_53(:,i)./kT_Qi_GM_53(:,1)./kN;
-kT_Xi_GENE (:,i) = kT_Qi_GENE (:,i)./kT_Qi_GENE (:,1)./kN;
+kT_Xi_GM_32 (:,i) = kT_Qi_GM_32 (:,i)./kT_Qi_GM_32 (:,1)./kN;
+kT_Xi_GM_53 (:,i) = kT_Qi_GM_53 (:,i)./kT_Qi_GM_53 (:,1)./kN;
+kT_Xi_GM_122(:,i) = kT_Qi_GM_122(:,i)./kT_Qi_GM_122(:,1)./kN;
+kT_Xi_GENE  (:,i) = kT_Qi_GENE  (:,i)./kT_Qi_GENE  (:,1)./kN;
 end
 %% Dimits fig 3 data
 KT_DIM      = [4.0 4.5 5.0 6.0 7.0 9.0 12. 14. 16. 18.];
@@ -80,14 +81,21 @@ GFL__98_DIM = [5.0 1.5;...
                7.0 5.0;...
                9.0 7.5];
 %% Plot
+msz = 8.0; lwt = 2.0;
 figure;
 if 1
 xye = kT_Xi_GM_32;
-errorbar(xye(:,1), xye(:,2),xye(:,3),'o-','DisplayName','192x96x16x3x2','LineWidth',2.0); hold on
+errorbar(xye(:,1), xye(:,2),xye(:,3),'>-','DisplayName','192x96x16x3x2',...
+    'MarkerSize',msz,'LineWidth',lwt); hold on
 xye = kT_Xi_GM_53;
-errorbar(xye(:,1), xye(:,2),xye(:,3),'o-','DisplayName','128x64x16x5x3','LineWidth',2.0); hold on
+errorbar(xye(:,1), xye(:,2),xye(:,3),'<-','DisplayName','128x64x16x5x3',...
+    'MarkerSize',msz,'LineWidth',lwt); hold on
+xye = kT_Xi_GM_122;
+errorbar(xye(:,1), xye(:,2),xye(:,3),'^-','DisplayName','128x64x16x13x3',...
+    'MarkerSize',msz,'LineWidth',lwt); hold on
 xye = kT_Xi_GENE;
-errorbar(xye(:,1), xye(:,2),xye(:,3),'v-.','DisplayName','GENE 128x64x16x24x12','LineWidth',2.0);
+errorbar(xye(:,1), xye(:,2),xye(:,3),'+-.k','DisplayName','GENE 128x64x16x24x12',...
+    'MarkerSize',msz,'LineWidth',lwt); hold on
 end
 if 1
    plot(LLNL_GK_DIM(:,1),LLNL_GK_DIM(:,2),'dk--','DisplayName','Dimits GK, LLNL'); hold on
diff --git a/wk/analysis_HeLaZ.m b/wk/analysis_HeLaZ.m
index da93ddf3fbbc9ce60f40961ff0730d311901dbec..1feabadb7f0c207bef247cf2b4f7f72b9382be01 100644
--- a/wk/analysis_HeLaZ.m
+++ b/wk/analysis_HeLaZ.m
@@ -22,8 +22,10 @@ FMT = '.fig';
 
 if 1
 %% Space time diagramm (fig 11 Ivanov 2020)
-options.TAVG_0   = 0.5*data.Ts3D(end); data.scale = 1;%/(data.Nx*data.Ny)^2;
-options.TAVG_1   = data.Ts3D(end); % Averaging times duration
+% data.scale = 1;%/(data.Nx*data.Ny)^2;
+options.TAVG_0   = 25;%0.4*data.Ts3D(end); 
+options.TAVG_1   = 40;%0.9*data.Ts3D(end); % Averaging times duration
+options.NCUT     = 4;              % Number of cuts for averaging and error estimation
 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)
@@ -49,13 +51,13 @@ options.NAME      = '\phi';
 % options.NAME      = 'n_i^{NZ}';
 % options.NAME      = '\Gamma_x';
 % options.NAME      = 'n_i';
-options.PLAN      = 'xy';
+options.PLAN      = 'kxky';
 % options.NAME      = 'f_i';
 % options.PLAN      = 'sx';
 options.COMP      = 'avg';
 % options.TIME      = data.Ts5D(end-30:end);
-options.TIME      =  data.Ts3D(1:2:end);
-% options.TIME      = [550:2:750];
+options.TIME      =  data.Ts3D;
+% options.TIME      = [850:0.1:1000];
 data.EPS          = 0.1;
 data.a = data.EPS * 2000;
 create_film(data,options,'.gif')
@@ -67,17 +69,17 @@ if 0
 options.INTERP    = 1;
 options.POLARPLOT = 0;
 options.AXISEQUAL = 0;
-% options.NAME      = '\phi';
+options.NAME      = '\phi';
 % options.NAME      = 'n_e';
-options.NAME      = 'N_i^{00}';
+% options.NAME      = 'N_i^{00}';
 % options.NAME      = 'T_i';
 % options.NAME      = '\Gamma_x';
 % options.NAME      = 'k^2n_e';
-options.PLAN      = 'xy';
+options.PLAN      = 'kxky';
 % options.NAME      'f_i';
 % options.PLAN      = 'sx';
 options.COMP      = 'avg';
-options.TIME      = [550 650];
+options.TIME      = [400 440];
 data.a = data.EPS * 2e3;
 fig = photomaton(data,options);
 % save_figure(data,fig)
@@ -101,8 +103,8 @@ if 0
 % options.XPERP     = linspace( 0,6,32);
 options.SPAR      = gene_data.vp';
 options.XPERP     = gene_data.mu';
-options.iz        = 1;
-options.T         = [500 1000];
+options.iz        = 'avg';
+options.T         = [300 600];
 options.PLT_FCT   = 'contour';
 options.ONED      = 0;
 options.non_adiab = 0;
@@ -130,10 +132,10 @@ end
 
 if 0
 %% Time averaged spectrum
-options.TIME   = [100 500];
+options.TIME   = [300 600];
 options.NORM   =1;
-% options.NAME   = '\phi';
-options.NAME      = 'N_i^{00}';
+options.NAME   = '\phi';
+% options.NAME      = 'N_i^{00}';
 % options.NAME   ='\Gamma_x';
 options.PLAN   = 'kxky';
 options.COMPZ  = 'avg';
diff --git a/wk/analysis_gene.m b/wk/analysis_gene.m
index 974681bdf66d94ff629fd986799f4335c61db17b..ac4dc8f0965ede26f6b172ca28a333352fd6089a 100644
--- a/wk/analysis_gene.m
+++ b/wk/analysis_gene.m
@@ -11,14 +11,15 @@
 % 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/Z-pinch/ZP_HP_kn_1.6_HRES/';
 % folder = '/misc/gene_results/ZP_kn_2.5_large_box/';
 % folder = '/misc/gene_results/CBC/128x64x16x24x12/';
-% folder = '/misc/gene_results/CBC/196x96x20x32x16_00/';
+% folder = '/misc/gene_results/CBC/196x96x20x32x16_01/';
 % folder = '/misc/gene_results/CBC/128x64x16x6x4/';
-% folder = '/misc/gene_results/CBC/KT_5.3_128x64x16x24x12_00/';
-% folder = '/misc/gene_results/CBC/KT_4.5_128x64x16x24x12_01/';
+% folder = '/misc/gene_results/CBC/KT_5.3_128x64x16x24x12_01/';
+folder = '/misc/gene_results/CBC/KT_4.5_128x64x16x24x12_01/';
 % folder = '/misc/gene_results/CBC/KT_11_128x64x16x24x12/';
-folder = '/misc/gene_results/CBC/KT_11_large_box_128x64x16x24x12/';
+% folder = '/misc/gene_results/CBC/KT_11_large_box_128x64x16x24x12/';
 gene_data = load_gene_data(folder);
 gene_data = invert_kxky_to_kykx_gene_results(gene_data);
 if 1
@@ -50,11 +51,11 @@ options.NAME      = '\phi';
 % options.NAME      = 'n_i';
 % options.NAME      = '\Gamma_x';
 % options.NAME      = 'k^2n_e';
-options.PLAN      = 'xy';
+options.PLAN      = 'kxky';
 % options.NAME      ='f_e';
 % options.PLAN      = 'sx';
 options.COMP      = 'avg';
-options.TIME      = [20 30 40 50 60];
+options.TIME      = [325 400];
 gene_data.a = data.EPS * 2000;
 fig = photomaton(gene_data,options);
 save_figure(gene_data,fig,'.png')
@@ -107,11 +108,11 @@ end
 
 if 0
 %% Show f_i(vpar,mu)
-options.times   = 20:80;
+options.times   = 250:500;
 options.specie  = 'i';
 options.PLT_FCT = 'contour';
 options.folder  = folder;
-options.iz      = 1;
+options.iz      = 'avg';
 options.FIELD   = '<f_>';
 options.ONED    = 0;
 % options.FIELD   = 'Q_es';
@@ -122,9 +123,9 @@ 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;
diff --git a/wk/header_3D_results.m b/wk/header_3D_results.m
index 8fb75f84d2c4d661758316950c554a01af586c7c..5ccf79448facc6d60186d0f25e437777194955b3 100644
--- a/wk/header_3D_results.m
+++ b/wk/header_3D_results.m
@@ -41,17 +41,23 @@ helazdir = '/home/ahoffman/HeLaZ/';
 % outfile = 'CBC/64x32x16x5x3';
 % outfile = 'CBC/64x128x16x5x3';
 % outfile = 'CBC/128x64x16x5x3';
+% outfile = 'CBC/128x96x16x3x2_Nexc_6';
 % outfile = 'CBC/192x96x16x3x2';
-% outfile = 'CBC/128x64x16x5x3';
-outfile = 'CBC/kT_11_128x64x16x5x3';
+% outfile = 'CBC/192x96x24x13x7';
+% outfile = 'CBC/kT_11_128x64x16x5x3';
 % outfile = 'CBC/kT_9_256x128x16x3x2';
 % outfile = 'CBC/kT_4.5_128x64x16x13x2';
+% outfile = 'CBC/kT_4.5_192x96x24x13x7';
+% outfile = 'CBC/kT_5.3_192x96x24x13x7';
 % outfile = 'CBC/kT_13_large_box_128x64x16x5x3';
-JOBNUMMIN = 00; JOBNUMMAX = 06;
+outfile = 'CBC/kT_13_96x96x16x3x2_Nexc_6';
+% outfile = 'CBC/kT_11_96x64x16x5x3_ky_0.02';
 
 % outfile = 'CBC/kT_scan_128x64x16x5x3';
 % outfile = 'CBC/kT_scan_192x96x16x3x2';
 
 %% Linear CBC
 % outfile = 'linear_CBC/20x2x32_21x11_Lx_62.8319_Ly_31.4159_q0_1.4_e_0.18_s_0.8_kN_2.22_kT_5.3_nu_1e-02_DGDK_adiabe';
+
+JOBNUMMIN = 00; JOBNUMMAX = 20;
 run analysis_HeLaZ
diff --git a/wk/quick_run.m b/wk/quick_run.m
index b1d4218e22a332d9e0717c695e9e9e80aff10010..d78aafece68574167c46bf960ffea03e33cd087f 100644
--- a/wk/quick_run.m
+++ b/wk/quick_run.m
@@ -3,8 +3,10 @@
 % from matlab framework. It is meant to run only small problems in linear
 % for benchmark and debugging purpose since it makes matlab "busy"
 %
-SIMID   = 'linear_CBC';  % Name of the simulation
-RUN     = 0; % To run or just to load
+% SIMID   = 'test_circular_geom';  % Name of the simulation
+% SIMID   = 'linear_CBC';  % Name of the simulation
+SIMID   = 'dbg';  % Name of the simulation
+RUN     = 1; % To run or just to load
 addpath(genpath('../matlab')) % ... add
 default_plots_options
 HELAZDIR = '/home/ahoffman/HeLaZ/';
@@ -19,23 +21,23 @@ NU      = 0.05;   % Collision frequency
 TAU     = 1.0;    % e/i temperature ratio
 K_N     = 2.22;%2.0;   % ion Density gradient drive
 K_Ne    = K_N;        % ele Density gradient drive
-K_T     = 11;%0.25*K_N;   % Temperature '''
+K_T     = 6.96;%0.25*K_N;   % Temperature '''
 K_Te     = K_T;            % Temperature '''
 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   = 0;     % 1: kinetic electrons, 2: adiabatic electrons
 BETA    = 0.0;     % electron plasma beta 
 %% GRID PARAMETERS
-P = 8;
-J = P/2;
+P = 4;
+J = 2;%P/2;
 PMAXE   = P;     % Hermite basis size of electrons
 JMAXE   = J;     % Laguerre "
 PMAXI   = P;     % " ions
 JMAXI   = J;     % "
-NX      = 12;    % real space x-gridpoints
-NY      = 10;     %     ''     y-gridpoints
+NX      = 32;    % real space x-gridpoints
+NY      = 16;     %     ''     y-gridpoints
 LX      = 2*pi/0.1;   % Size of the squared frequency domain
-LY      = 2*pi/0.1;     % Size of the squared frequency domain
+LY      = 2*pi/0.05;     % Size of the squared frequency domain
 NZ      = 16;    % number of perpendicular planes (parallel grid)
 NPOL    = 1;
 SG      = 0;     % Staggered z grids option
@@ -44,8 +46,9 @@ SG      = 0;     % Staggered z grids option
 GEOMETRY= 's-alpha';
 % GEOMETRY= 'circular';
 Q0      = 1.4;    % safety factor
-SHEAR   = 0.8;    % magnetic shear (Not implemented yet)
-EPS     = 0.18;    % inverse aspect ratio
+SHEAR   = 0.8;    % magnetic shear
+NEXC    = 6;      % To extend Lx if needed (Lx = Nexc/(kymin*shear))
+EPS     = 0.18;   % inverse aspect ratio
 %% TIME PARMETERS
 TMAX    = 25;  % Maximal time unit
 DT      = 1e-2;   % Time step
@@ -126,7 +129,7 @@ end
 if 0
 %% Ballooning plot
 options.time_2_plot = [120];
-options.kymodes     = [0.1 0.2 0.3];
+options.kymodes     = [0.1];
 options.normalized  = 1;
 options.field       = 'phi';
 fig = plot_ballooning(data,options);