From 8787f7eab857e45d22b6d36d49e19689c92b32e6 Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Thu, 9 Sep 2021 15:20:32 +0200
Subject: [PATCH] post processing scripts update

---
 matlab/LinearFit_s.m                          |   1 +
 matlab/compile_results.m                      |   4 +
 matlab/create_gif.m                           |   1 +
 matlab/load_params.m                          |   2 +-
 matlab/plot_param_evol.m                      |  33 +++-
 matlab/plots/plot_kperp_spectrum.m            |   7 +-
 .../plots/plot_radial_transport_and_shear.m   |   6 +-
 matlab/plots/plot_space_time_diagrams.m       |  16 +-
 matlab/setup.m                                |   7 +-
 matlab/write_fort90.m                         |   3 +-
 wk/HD_study.m                                 |  40 ++--
 wk/ZF_fourier_analysis.m                      |  94 +++------
 wk/analysis_3D.m                              |  87 ++++++---
 wk/continue_multiple_runs_marconi.m           |  20 +-
 wk/linear_study.m                             | 179 ++++++------------
 wk/load_multiple_outputs_marconi.m            |   4 +
 wk/local_run.m                                |   2 +-
 wk/marconi_run.m                              |  16 +-
 wk/new_flux_results.m                         |  31 +++
 wk/open_figure_script.m                       |   8 +-
 wk/plot_cosol_mat.m                           |   5 +-
 21 files changed, 300 insertions(+), 266 deletions(-)
 create mode 100644 wk/new_flux_results.m

diff --git a/matlab/LinearFit_s.m b/matlab/LinearFit_s.m
index dc424911..dd5f1129 100644
--- a/matlab/LinearFit_s.m
+++ b/matlab/LinearFit_s.m
@@ -30,5 +30,6 @@ function [gamma,fit] = LinearFit_s(time,Na00abs)
 
     % Return gamma(t) for amplitude ratio method
     fit.gammaoft = gammaoft;
+    fit.t        = time;
 
 end % ... end function
diff --git a/matlab/compile_results.m b/matlab/compile_results.m
index 44177cf2..8e367c4e 100644
--- a/matlab/compile_results.m
+++ b/matlab/compile_results.m
@@ -2,8 +2,10 @@ CONTINUE = 1;
 JOBNUM   = JOBNUMMIN; JOBFOUND = 0;
 TJOB_SE  = []; % Start and end times of jobs
 NU_EVOL  = []; % evolution of parameter nu between jobs
+CO_EVOL  = []; % evolution of CO
 MU_EVOL  = []; % evolution of parameter mu between jobs
 ETAN_EVOL= []; %
+ETAB_EVOL= []; %
 L_EVOL   = []; % 
 DT_EVOL  = []; %
 % FIELDS
@@ -112,8 +114,10 @@ while(CONTINUE)
         load_params
         TJOB_SE   = [TJOB_SE Ts0D(1) Ts0D(end)]; 
         NU_EVOL   = [NU_EVOL NU NU];
+        CO_EVOL   = [CO_EVOL CO CO];
         MU_EVOL   = [MU_EVOL MU MU];
         ETAN_EVOL = [ETAN_EVOL ETAN ETAN];
+        ETAB_EVOL = [ETAB_EVOL ETAB ETAB];
         L_EVOL    = [L_EVOL L L];
         DT_EVOL   = [DT_EVOL DT_SIM DT_SIM];
     
diff --git a/matlab/create_gif.m b/matlab/create_gif.m
index cd020f8d..299d65cd 100644
--- a/matlab/create_gif.m
+++ b/matlab/create_gif.m
@@ -32,6 +32,7 @@ fig  = figure('Color','white','Position', [100, 100, 400, 400]);
         end
         set(pclr, 'edgecolor','none'); axis square;
         caxis([-1,1]);
+        colormap(bluewhitered)
         xlabel(XNAME); ylabel(YNAME); %colorbar;
         title([FIELDNAME,', $t \approx$', sprintf('%.3d',ceil(T(n)))...
             ,', scaling = ',sprintf('%.1e',scale)]);
diff --git a/matlab/load_params.m b/matlab/load_params.m
index 336c3057..12819418 100644
--- a/matlab/load_params.m
+++ b/matlab/load_params.m
@@ -30,7 +30,7 @@ end
 if    (CO == -3); CONAME = 'PADK';
 elseif(CO == -2); CONAME = 'SGDK';
 elseif(CO == -1); CONAME = 'DGDK';
-elseif(CO ==  0); CONAME = 'LB';
+elseif(CO ==  0); CONAME = 'LBGK';
 elseif(CO ==  1); CONAME = 'DGGK';
 elseif(CO ==  2); CONAME = 'SGGK';
 elseif(CO ==  3); CONAME = 'PAGK';
diff --git a/matlab/plot_param_evol.m b/matlab/plot_param_evol.m
index fc3c8f8e..718bad2b 100644
--- a/matlab/plot_param_evol.m
+++ b/matlab/plot_param_evol.m
@@ -1,18 +1,39 @@
-figure; hold on; set(gcf, 'Position',  [100, 100, 400, 1500])
+fig = figure; FIGNAME = 'linear_study';
+hold on; set(gcf, 'Position',  [100, 100, 400, 1500])
 subplot(4,1,1);
-    plot(TJOB_SE,NU_EVOL,'DisplayName','$\nu$');
+    title('Parameter evolution'); hold on;
+    yyaxis left
+    plot(TJOB_SE,NU_EVOL,'DisplayName','$\nu$'); 
+    yyaxis right
+    plot(TJOB_SE,CO_EVOL,'DisplayName','CO'); 
     xlim([TJOB_SE(1) TJOB_SE(end)]);legend('show');grid on;
-    
+    xticks([]);
+    plot_tjob_lines(TJOB_SE,ylim)    
 subplot(4,1,2);
-    plot(TJOB_SE,MU_EVOL,'DisplayName','$\mu$');
+    plot(TJOB_SE,MU_EVOL,'DisplayName','$\mu$'); hold on;
     xlim([TJOB_SE(1) TJOB_SE(end)]);legend('show');grid on;
+    xticks([]);
+    plot_tjob_lines(TJOB_SE,ylim)    
     
 subplot(4,1,3);
-    plot(TJOB_SE,1./ETAN_EVOL,'DisplayName','$\eta$');
+    yyaxis left
+    plot(TJOB_SE,ETAN_EVOL,'DisplayName','$\nabla n$'); hold on;
+    yyaxis right
+    plot(TJOB_SE,ETAB_EVOL,'DisplayName','$\nabla B$'); hold on;
     xlim([TJOB_SE(1) TJOB_SE(end)]);legend('show');grid on;
+    xticks([]);
+    plot_tjob_lines(TJOB_SE,ylim)    
 
 subplot(4,1,4);
-    plot(TJOB_SE,L_EVOL,'DisplayName','$L$');
+    plot(TJOB_SE,L_EVOL,'DisplayName','$L$'); hold on;
     xlim([TJOB_SE(1) TJOB_SE(end)]);legend('show');grid on;
+    plot_tjob_lines(TJOB_SE,ylim)    
 
  xlabel('$t c_s/R$');
+ saveas(fig,[BASIC.RESDIR,'param_evol.png']);
+ 
+ function [] = plot_tjob_lines(TJOB,limits)
+     for i = 2:numel(TJOB)-1
+        plot(TJOB(i)*[1 1],limits,'--k') 
+     end
+ end
diff --git a/matlab/plots/plot_kperp_spectrum.m b/matlab/plots/plot_kperp_spectrum.m
index 238faa03..43eadb46 100644
--- a/matlab/plots/plot_kperp_spectrum.m
+++ b/matlab/plots/plot_kperp_spectrum.m
@@ -18,7 +18,12 @@ a1 = field_kp(2)*kp_ip(2).^(13/3);
 a2 = field_kp(2)*kp_ip(2).^(3)./(1+kp_ip(2).^2).^(-2);
 fig = figure; FIGNAME = ['cascade','_',FNAME,'_',PARAMS];set(gcf, 'Position',  [100, 100, 800, 300])
 % scatter(kperp,phi_k_2,'.k','MarkerEdgeAlpha',0.4,'DisplayName','$|\phi_k|^2$'); hold on; grid on;
-plot(kp_ip,field_kp,'^','DisplayName',['$\langle|',FIELDLTX,'|^2\rangle_{k_\perp}$']); hold on;
+if NORMALIZED
+   plt = @(x) x./max(x);
+else
+   plt = @(x) x;
+end
+plot(kp_ip,plt(field_kp),'^','DisplayName',['$\langle|',FIELDLTX,'|^2\rangle_{k_\perp}$']); hold on;
 if TRENDS
 plot(kp_ip,a1*kp_ip.^(-13/3),'-','DisplayName','$k^{-13/3}$');
 plot(kp_ip,a2/100*kp_ip.^(-3)./(1+kp_ip.^2).^2,'-','DisplayName','$k^{-3}/(1+k^2)^2$');
diff --git a/matlab/plots/plot_radial_transport_and_shear.m b/matlab/plots/plot_radial_transport_and_shear.m
index 734551cc..dc94c2cf 100644
--- a/matlab/plots/plot_radial_transport_and_shear.m
+++ b/matlab/plots/plot_radial_transport_and_shear.m
@@ -1,12 +1,12 @@
 %Compute steady radial transport
-tend = Ts0D(end); tstart = tend - TAVG;
+tend = TAVG_1; tstart = TAVG_0;
 [~,its0D] = min(abs(Ts0D-tstart));
 [~,ite0D]   = min(abs(Ts0D-tend));
 SCALE = (2*pi/Nx/Ny)^2;
 gamma_infty_avg = mean(PGAMMA_RI(its0D:ite0D))*SCALE;
 gamma_infty_std = std (PGAMMA_RI(its0D:ite0D))*SCALE;
 % Compute steady shearing rate
-tend = Ts3D(end); tstart = tend - TAVG;
+tend = TAVG_1; tstart = TAVG_0;
 [~,its2D] = min(abs(Ts3D-tstart));
 [~,ite2D]   = min(abs(Ts3D-tend));
 shear_infty_avg = mean(mean(shear_maxx_avgy(:,its2D:ite2D),1));
@@ -22,7 +22,7 @@ fig = figure; FIGNAME = ['ZF_transport_drphi','_',PARAMS];set(gcf, 'Position',
         ylim([0,5*abs(gamma_infty_avg)]); xlim([Ts0D(1),Ts0D(end)]);
         title(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta=$',num2str(ETAB/ETAN),...
         ', $L=',num2str(L),'$, $N=',num2str(Nx),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',...
-        ' $\mu_{hd}=$',num2str(MU)]);
+        ' $\mu_{hd}=$',num2str(MU),', $\Gamma^{\infty} \approx $',num2str(gamma_infty_avg)]);
         %         
     subplot(312)
         clr      = line_colors(1,:);
diff --git a/matlab/plots/plot_space_time_diagrams.m b/matlab/plots/plot_space_time_diagrams.m
index 3e263a0a..582f5220 100644
--- a/matlab/plots/plot_space_time_diagrams.m
+++ b/matlab/plots/plot_space_time_diagrams.m
@@ -4,21 +4,23 @@ trange = itstart:itend;
 [TY,TX] = meshgrid(x,Ts3D(trange));
 fig = figure; FIGNAME = ['space_time','_',PARAMS];set(gcf, 'Position',  [100, 100, 1200, 600])
     subplot(211)
-%         pclr = pcolor(TX,TY,squeeze(mean(dens_i(:,:,trange).*dzphi(:,:,trange),2))'); set(pclr, 'edgecolor','none'); colorbar;
-        pclr = pcolor(TX,TY,squeeze(mean(ni00(:,:,trange).*dzphi(:,:,trange),2))'); set(pclr, 'edgecolor','none'); colorbar;
+%         pclr = pcolor(TX,TY,squeeze(mean(dens_i(:,:,trange).*dyphi(:,:,trange),2))'); set(pclr, 'edgecolor','none'); colorbar;
+        pclr = pcolor(TX,TY,squeeze(mean(ni00(:,:,1,trange).*dyphi(:,:,1,trange),2))'); set(pclr, 'edgecolor','none',...
+            'DisplayName','$\langle n_i\partial_z\phi\rangle_z$'); colorbar;
         shading interp
         colormap hot;
-        caxis([0.0,0.05*max(max(mean(ni00(:,:,its2D:ite2D).*dzphi(:,:,its2D:ite2D),2)))]);
-        caxis([0.0,cmax]); c = colorbar; c.Label.String ='\langle n_i\partial_z\phi\rangle_z';
+%         caxis([0.0,0.05*max(max(mean(ni00(:,:,its2D:ite2D).*dyphi(:,:,1,its2D:ite2D),2)))]);
+        caxis([0.0,cmax]); c = colorbar; c.Label.String ='\langle\Gamma_{x}\rangle_{z}';
          xticks([]); ylabel('$x/\rho_s$')
 %         legend('Radial part. transport $\langle n_i\partial_z\phi\rangle_z$')
         title(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta=$',num2str(ETAB/ETAN),...
         ', $L=',num2str(L),'$, $N=',num2str(Nx),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',...
         ' $\mu_{hd}=$',num2str(MU)]);
     subplot(212)
-        pclr = pcolor(TX,TY,squeeze(mean(drphi(:,:,1,trange),2))'); set(pclr, 'edgecolor','none'); colorbar;
-        fieldmax = max(max(mean(abs(drphi(:,:,1,its2D:ite2D)),2)));
-        caxis([-fieldmax,fieldmax]); c = colorbar; c.Label.String ='\langle \partial_r\phi\rangle_z';
+        pclr = pcolor(TX,TY,squeeze(mean(dxphi(:,:,1,trange),2))'); set(pclr, 'edgecolor','none',...
+            'DisplayName','$\langle \partial_r\phi\rangle_z$'); colorbar;
+        fieldmax = max(max(mean(abs(dxphi(:,:,1,its2D:ite2D)),2)));
+        caxis([-fieldmax,fieldmax]); c = colorbar; c.Label.String ='\langle v_{E\times B,z}\rangle_z';
         xlabel('$t c_s/R$'), ylabel('$x/\rho_s$')
 %         legend('Zonal flow $\langle \partial_r\phi\rangle_z$')        
 save_figure
\ No newline at end of file
diff --git a/matlab/setup.m b/matlab/setup.m
index b04e92c6..f753c54d 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -43,9 +43,10 @@ MODEL.lambdaD = LAMBDAD;
 TIME_INTEGRATION.numerical_scheme  = '''RK4''';
 if (INIT_PHI); INITIAL.init_noisy_phi = '.true.'; else; INITIAL.init_noisy_phi = '.false.';end;
 INITIAL.INIT_ZF = INIT_ZF;
-if (WIPE_TURB); INITIAL.wipe_turb = '.true.'; else; INITIAL.wipe_turb = '.false.';end;
+INITIAL.wipe_turb = WIPE_TURB;
+INITIAL.wipe_zf   = WIPE_ZF;
 if (INIT_BLOB); INITIAL.init_blob = '.true.'; else; INITIAL.init_blob = '.false.';end;
-INITIAL.init_background  = (INIT_ZF>0)*ZF_AMP;
+INITIAL.init_background  = (INIT_ZF>0)*ZF_AMP + BCKGD0;
 INITIAL.init_noiselvl = NOISE0;
 INITIAL.iseed             = 42;
 INITIAL.mat_file   = '''null''';
@@ -70,7 +71,7 @@ switch abs(CO)
     case 4; CONAME = 'FC';
     otherwise; CONAME ='UK';
 end
-if (CO <= 0); CONAME = [CONAME,'DK'];
+if (CO < 0); CONAME = [CONAME,'DK'];
 else;         CONAME = [CONAME,'GK'];
 end
 if    (CLOS == 0); CLOSNAME = 'Trunc.';
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index dcf61884..d23fdec6 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -71,7 +71,8 @@ fprintf(fid,'/\n');
 fprintf(fid,'&INITIAL_CON\n');
 fprintf(fid,['  INIT_NOISY_PHI    = ', INITIAL.init_noisy_phi,'\n']);
 fprintf(fid,['  INIT_ZF       = ', num2str(INITIAL.INIT_ZF),'\n']);
-fprintf(fid,['  WIPE_TURB     = ', INITIAL.wipe_turb,'\n']);
+fprintf(fid,['  WIPE_ZF       = ', num2str(INITIAL.wipe_zf),'\n']);
+fprintf(fid,['  WIPE_TURB     = ', num2str(INITIAL.wipe_turb),'\n']);
 fprintf(fid,['  INIT_BLOB     = ', INITIAL.init_blob,'\n']);
 fprintf(fid,['  init_background  = ', num2str(INITIAL.init_background),'\n']);
 fprintf(fid,['  init_noiselvl = ', num2str(INITIAL.init_noiselvl),'\n']);
diff --git a/wk/HD_study.m b/wk/HD_study.m
index ca79206e..9ccde4cd 100644
--- a/wk/HD_study.m
+++ b/wk/HD_study.m
@@ -1,5 +1,9 @@
+red_ = [0.6350 0.0780 0.1840];
+gre_ = [0.4660 0.6740 0.1880];
+
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% DOUGHERTY
+% nu mu
 % bursts simulations
 b_=[...
     5e-1, 5e+0;...
@@ -10,7 +14,9 @@ b_=[...
     5e-2, 3e-3;...
     5e-2, 2e-3;...
     1e-2, 1e-2;... % v2.7_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e-02_DGGK_CLOS_0_mu_1e-02/
-    1e-2, 3e-3;...
+    1e-2, 3e-3;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_5e-02_DGGK_mu_3e-03
+    1e-2, 5e-4;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-02_DGGK_mu_5e-04
+    1e-2, 1e-4;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-02_DGGK_mu_3e-03
     1e-3, 2e-2;... % v2.7_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.6_nu_1e-03_DGGK_CLOS_0_mu_2e-02/
     ];
 % converged turb plateau simulations
@@ -24,7 +30,6 @@ cp_=[...
     1e-1, 1e-3;...
     5e-2, 1e-3;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-01_DGGK_mu_3e-02
     5e-2, 5e-4;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_5e-02_DGGK_mu_5e-04
-    1e-2, 5e-4;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-02_DGGK_mu_5e-04
     ];
 % moving/no turb plateau
 dp_=[...
@@ -39,7 +44,7 @@ dp_=[...
     ];
 % not sure
 rp_=[...
-    1e-2, 1e-4;...
+    0,0;...
     ];
 figure; set(gcf, 'Position',  [100, 100, 900, 400])
 title('Hyperdiffusion study, Dougherty GK')
@@ -95,21 +100,22 @@ nu_ = [1e-1];
 % plot(nu_,mu_,'x--','DisplayName','N=200, L=120, P,J=6,3');
 end
 scatter( b_(:,1), b_(:,2),'o',...
-    'MarkerFaceColor',[0.6350 0.0780 0.1840],'MarkerEdgeColor',[0 0 0],'SizeData',50,...
+    'MarkerFaceColor',red_,'MarkerEdgeColor',[0 0 0],'SizeData',50,...
     'DisplayName','Bursts'); 
 scatter(cp_(:,1),cp_(:,2),'s',...
-    'MarkerFaceColor',[0.4660 0.6740 0.1880],'MarkerEdgeColor',[0 0 0],'SizeData',50,...
+    'MarkerFaceColor',gre_,'MarkerEdgeColor',[0 0 0],'SizeData',50,...
     'DisplayName','Converged Plateau'); 
 scatter(dp_(:,1),dp_(:,2),'d',...
-    'MarkerFaceColor',[0.4660 0.6740 0.1880],'MarkerEdgeColor',[0 0 0],'SizeData',50,...
+    'MarkerFaceColor',gre_,'MarkerEdgeColor',[0 0 0],'SizeData',50,...
     'DisplayName','Moving Plateau'); 
 scatter(rp_(:,1),rp_(:,2),'h',...
     'MarkerFaceColor',[0.9290 0.6940 0.1250],'MarkerEdgeColor',[0 0 0],'SizeData',50,...
     'DisplayName','not sure'); 
 plot(0,0,'v','MarkerFaceColor',[0 0 0],'MarkerEdgeColor',[0 0 0], 'DisplayName','$\mu=0$');
 legend('show','Location','NorthWest')
-scatter(1,5e-5,80,'v','MarkerFaceColor',[0.4660 0.6740 0.1880],'MarkerEdgeColor',[0 0 0]);
-
+scatter(1,5e-5,80,'v','MarkerFaceColor',gre_,'MarkerEdgeColor',[0 0 0]);
+% HD_study/150x75_L_100_P_2_J_1_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_3e-02/
+scatter(0.1,5e-5,80,'v','MarkerFaceColor',gre_,'MarkerEdgeColor',[0 0 0]);
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% SUGAMA
 % nu mu
@@ -121,7 +127,7 @@ b_=[...
     5e-1, 1.6e-2;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_5e-01_SGGK_mu_3e-02
     5e-1, 0;...      % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_5e-01_SGGK_mu_0e+00/out.txt
     1e-1, 2.0e-2;... % v2.7_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e-01_SGGK_CLOS_0_mu_2e-02/
-    1e-1, 1.6e-2;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-01_DGGK_mu_3e-02/
+    1e-1, 1.6e-2;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-01_GGK_mu_3e-02/
     1e-2, 5.0e-3;... % v2.7_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e-02_SGGK_CLOS_0_mu_1e-02/
     1e-2, 3.0e-3;... % v2.7_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e-02_SGGK_CLOS_0_mu_1e-02/
     ];
@@ -167,13 +173,13 @@ plot(nu_,mu_,'x--','DisplayName','N=200, L=050, P,J=2,1');
 end
 
 scatter( b_(:,1), b_(:,2),'o',...
-    'MarkerFaceColor',[0.6350 0.0780 0.1840],'MarkerEdgeColor',[0 0 0],'SizeData',80,...
+    'MarkerFaceColor',red_,'MarkerEdgeColor',[0 0 0],'SizeData',80,...
     'DisplayName','Bursts'); 
 % scatter(cp_(:,1),cp_(:,2),'s',...
-%     'MarkerFaceColor',[0.4660 0.6740 0.1880],'MarkerEdgeColor',[0 0 0],'SizeData',60,...
+%     'MarkerFaceColor',gre_,'MarkerEdgeColor',[0 0 0],'SizeData',60,...
 %     'DisplayName','Converged Plateau'); 
 scatter(dp_(:,1),dp_(:,2),'d',...
-    'MarkerFaceColor',[0.4660 0.6740 0.1880],'MarkerEdgeColor',[0 0 0],'SizeData',60,...
+    'MarkerFaceColor',gre_,'MarkerEdgeColor',[0 0 0],'SizeData',60,...
     'DisplayName','Plateau'); 
 % scatter(rp_(:,1),rp_(:,2),'h',...
 %     'MarkerFaceColor',[0.9290 0.6940 0.1250],'MarkerEdgeColor',[0 0 0],'SizeData',60,...
@@ -182,7 +188,13 @@ scatter(0,0,'v','MarkerFaceColor',[0 0 0],'MarkerEdgeColor',[0 0 0],...
     'DisplayName','$\mu=0$');
 legend('show','Location','NorthWest')
 
-scatter(0.1,5e-5,80,'v','MarkerFaceColor',[0.4660 0.6740 0.1880],'MarkerEdgeColor',[0 0 0],...
+scatter(0.075,5e-5,80,'v','MarkerFaceColor',gre_,'MarkerEdgeColor',[0 0 0],...
+    'DisplayName','$\mu=0$');
+scatter(0.1,5e-5,80,'v','MarkerFaceColor',gre_,'MarkerEdgeColor',[0 0 0],...
+    'DisplayName','$\mu=0$');
+scatter(0.25,5e-5,80,'v','MarkerFaceColor',red_,'MarkerEdgeColor',[0 0 0],...
+    'DisplayName','$\mu=0$');
+scatter(0.5,5e-5,80,'v','MarkerFaceColor',red_,'MarkerEdgeColor',[0 0 0],...
     'DisplayName','$\mu=0$');
-scatter(0.5,5e-5,80,'v','MarkerFaceColor',[0.6350 0.0780 0.1840],'MarkerEdgeColor',[0 0 0],...
+scatter(1.0,5e-5,80,'v','MarkerFaceColor',red_,'MarkerEdgeColor',[0 0 0],...
     'DisplayName','$\mu=0$');
diff --git a/wk/ZF_fourier_analysis.m b/wk/ZF_fourier_analysis.m
index da4cf878..b64a5af4 100644
--- a/wk/ZF_fourier_analysis.m
+++ b/wk/ZF_fourier_analysis.m
@@ -1,10 +1,11 @@
 %% Zonal flow spectral analysis
 fig = figure; FIGNAME = ['zonal_flow_spectral_analysis_',PARAMS];
-tend = Ts0D(end); tstart = tend-TAVG ;
+tend = TAVG_1; tstart = TAVG_0;
 [~,its0D] = min(abs(Ts0D-tstart));
 [~,ite0D]   = min(abs(Ts0D-tend));
-[~,its2D] = min(abs(Ts2D-tstart));
-[~,ite2D]   = min(abs(Ts2D-tend));
+[~,its3D] = min(abs(Ts3D-tstart));
+[~,ite3D]   = min(abs(Ts3D-tend));
+TAVG = Ts3D(ite3D)-Ts3D(its3D);
 set(gcf, 'Position',  [100, 100, 800, 400])
     % Time series analysis (burst period and time frequencies spectrum)
     subplot(121)
@@ -17,87 +18,54 @@ set(gcf, 'Position',  [100, 100, 800, 400])
     Pot(1) = 0;
     nmax = min(20,round(n/2));
     [amax, itmax] = max(Pot);
-    plot((0:nmax-1) , Pot(1:nmax)/amax,'DisplayName','$\Gamma_r(\omega)$');hold on;
-    plot([itmax-1,itmax-1],[0,1],'--k', 'DisplayName',['$T_{per}\approx',num2str(round(1/Freq(itmax))),'L_\perp/c_s$']);
-    legend('show'); grid on; box on; xlabel('Period number'); yticks([]);
-    title('$\Gamma_r$ temporal spectrum')
+%     plot((0:nmax-1) , Pot(1:nmax)/amax,'DisplayName','$\Gamma_x(\omega)$');hold on;
+%     plot([itmax-1,itmax-1],[0,1],'--k', 'DisplayName',['$T_{per}\approx',num2str(round(1/Freq(itmax))),'L_\perp/c_s$']);
+    semilogx(TAVG./(0:nmax-1) , Pot(1:nmax)/amax,'o-','DisplayName','$\Gamma_x(\omega)$');hold on;
+    semilogx(TAVG./[itmax-1,itmax-1],[0,1],'--k', 'DisplayName',['$T_{max}\approx',num2str(round(1/Freq(itmax-1))),'L_\perp/c_s$']);
+    legend('show'); grid on; box on; xlabel('Period $Tc_s/R$'); yticks([]);
+    title('$\Gamma_x$ temporal spectrum')
     % Space analysis (spatial period of ZF)
     subplot(122)
-    nmax = 20; n = numel(r);
-    [TT,NN] = meshgrid(Ts2D(its2D:ite2D),0:n-1);
+    nmax = 20; n = numel(x);
+    [TT,NN] = meshgrid(Ts3D(its3D:ite3D),0:n-1);
     Pot = NN;
-    for it = 1:ite2D-its2D+1
-        Y = mean(real(drphi(:,:,it)),2);
+    for it = 1:ite3D-its3D+1
+        Y = mean(real(dxphi(:,:,it)),2);
         Yy = fft(Y);
         [n,~] = size(Yy);
         Pot(:,it) = Yy .* conj(Yy) / n;
     end
     [amax, ikZF] = max(mean(Pot,2));
 %     pclr = pcolor(NN(1:nmax,:),TT(1:nmax,:),Pot(1:nmax,:)); set(pclr, 'edgecolor','none'); hold on;
-    plot(0:nmax,mean(Pot(1:nmax+1,:),2)/amax,'DisplayName','$\langle\partial_r\phi\rangle_z (k_r)$'); hold on;
-    plot([ikZF-1,ikZF-1],[0,1],'--k', 'DisplayName',['$L_z=',num2str(2*pi/kx(ikZF)),'\rho_s$']);
+%     plot(0:nmax,mean(Pot(1:nmax+1,:),2)/amax,'DisplayName','$\langle\partial_x\phi\rangle_y (k_x)$'); hold on;
+%     plot([ikZF-1,ikZF-1],[0,1],'--k', 'DisplayName',['$L_x=',num2str(2*pi/kx(ikZF)),'\rho_s$']);
+    semilogx(Lx./(0:nmax),mean(Pot(1:nmax+1,:),2)/amax,'o-','DisplayName','$\langle\partial_x\phi\rangle_y (k_x)$'); hold on;
+    semilogx(Lx./[ikZF-1,ikZF-1],[0,1],'--k', 'DisplayName',['$L_x=',num2str(2*pi/kx(ikZF)),'\rho_s$']);
     grid on; box on;
     title('ZF spatial spectrum')
-    xlabel('radial mode number');  yticks([]); legend('show')
+    xlabel('Period $\lambda/\rho_s$');  yticks([]); legend('show')
 save_figure
 
 %% Pred-Pray phase space (A Zonal Flow review, Diamond 2005, Fig 15, Kobayashi 2015)
-
-E_turb           = zeros(1,Ns2D);    % Time evol. of the turbulence energy (Pred in Kobayashi 2015)
-E_ZF             = zeros(1,Ns2D);    % Time evol. of the ZF energy (Pray in Kobayashi 2015)
-for it = 1:numel(Ts2D)
-    E_turb(it) = sum(sum((1+kx.^2+ky.^2).*abs(PHI(:,:,it)).^2))- sum((1+kx.^2).*abs(PHI(:,1,it)).^2);
-    E_ZF(it)   = kx(ikZF)^2*abs(PHI(ikZF,1,it)).^2;
+% Time evol. of the turbulence energy (Pred in Kobayashi 2015, N = sum phi_k^2 (1+k^2) Non zonal)
+E_turb           = zeros(1,Ns3D);
+% Time evol. of the ZF energy (Pray in Kobayashi 2015, Ev = phi_q^2 q^2)
+E_ZF             = zeros(1,Ns3D);    
+for it = 1:numel(Ts3D)
+    E_turb(it) = sum(sum(((KY~=0).*(1+KX.^2+KY.^2).*abs(PHI(:,:,1,it)).^2)));
+%     E_ZF(it)   = kx(ikZF)^2*abs(PHI(ikZF,1,1,it)).^2;
+    E_ZF(it)   = sum(sum(((KY==0).*(1+KX.^2+KY.^2).*abs(PHI(:,:,1,it)).^2)));
 end
 fig = figure; FIGNAME = ['phi_shear_phase_space_',PARAMS];
 set(gcf, 'Position',  [100, 100, 700, 500])
-scatter(E_ZF*SCALE,E_turb*SCALE,35,Ts2D,'.',...
+scatter(E_ZF*SCALE,E_turb*SCALE,80,Ts3D,'.',...
     'DisplayName',PARAMS); cbar = colorbar;ylabel(cbar,'$t c_s/\rho_s$','Interpreter','LaTeX')
 hold on
 % xlabel('$\langle \phi \rangle_z^r$'); ylabel('$\langle dV_E/dr \rangle_z^r$')
 xlabel('$E_v$'); ylabel('$N$')
 grid on; title('ES pot. vs Shear phase space')
-% plot(phi_avgr_maxz(its2D:ite2D),shear_avgr_maxz(its2D:ite2D),'-')
-% plot(phi_maxr_maxz(its2D:ite2D),shear_maxr_maxz(its2D:ite2D),'-')
-% plot(phi_avgr_avgz(its2D:ite2D),shear_avgr_avgz(its2D:ite2D),'-')
+% plot(phi_avgr_maxz(iTs3D:ite2D),shear_avgr_maxz(iTs3D:ite2D),'-')
+% plot(phi_maxr_maxz(iTs3D:ite2D),shear_maxr_maxz(iTs3D:ite2D),'-')
+% plot(phi_avgr_avgz(iTs3D:ite2D),shear_avgr_avgz(iTs3D:ite2D),'-')
 save_figure
 clear x_ y_
-if 0
-%% density and phi phase space
-fig = figure; FIGNAME = ['phi_ni_phase_space_',PARAMS];
-set(gcf, 'Position',  [100, 100, 700, 500])
-t1 = Ts2D(end); t0 = 0;
-[~,its2D] = min(abs(Ts2D-t0)); [~,ite2D] = min(abs(Ts2D-t1));
-scatter3(max(mean(ni00(:,:,its2D:ite2D),2),[],1),phi_maxr_avgz(its2D:ite2D),shear_maxr_avgz(its2D:ite2D),35,Ts2D(its2D:ite2D),'.',...
-    'DisplayName',PARAMS); cbar = colorbar;ylabel(cbar,'$t c_s/\rho_s$','Interpreter','LaTeX')
-hold on
-xlabel('$\langle n_i^{00} \rangle_z^r$'); ylabel('$\langle \phi \rangle_z^r$'); zlabel('$\langle dV_E/dr \rangle_z^r$')
-grid on; title('ES pot. vs Shear phase space')
-% plot(phi_avgr_maxz(its2D:ite2D),shear_avgr_maxz(its2D:ite2D),'-')
-% plot(phi_maxr_maxz(its2D:ite2D),shear_maxr_maxz(its2D:ite2D),'-')
-% plot(phi_avgr_avgz(its2D:ite2D),shear_avgr_avgz(its2D:ite2D),'-')
-% save_figure
-end
-%% Non zonal quantities
-PHI_NZ = PHI;
-PHI_NZ(ikZF-1:ikZF+1,:,:) = 0;
-
-phi_nz    = zeros(Nx,Ny,Ns2D);
-for it = 1:numel(Ts2D)
-    PH_ = PHI_NZ(:,:,it);
-    phi_nz (:,:,it)  = real(fftshift(ifft2((PH_),Nx,Ny)));
-end
-%%
-t0    = 1000;
-[~, it02D] = min(abs(Ts2D-t0));
-[~, it05D] = min(abs(Ts5D-t0));
-skip_ = 10; 
-DELAY = 0.005*skip_;
-FRAMES_2D = it02D:skip_:numel(Ts2D);
-if 0
-%% Phi non zonal real space
-GIFNAME = ['phi_nz',sprintf('_%.2d',JOBNUM),'_',PARAMS];INTERP = 0;
-FIELD = real(phi_nz); X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D;
-FIELDNAME = '$\phi_{Ny}$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$';
-create_gif
-end
\ No newline at end of file
diff --git a/wk/analysis_3D.m b/wk/analysis_3D.m
index 63441138..c0bd1598 100644
--- a/wk/analysis_3D.m
+++ b/wk/analysis_3D.m
@@ -6,9 +6,11 @@ outfile ='';
 outfile ='';
 outfile ='';
 outfile ='';
-outfile ='HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-02_DGGK_mu_3e-03';
+outfile ='simulation_A/cw_SGGK_mu_1e-2';
+% outfile ='simulation_A/LBDK_damping_150x75_L_100_P_4_J_2_eta_0.6_nu_1e-01_SGGK_mu_0e+00';
     BASIC.RESDIR      = ['../results/',outfile,'/'];
     BASIC.MISCDIR     = ['/misc/HeLaZ_outputs/results/',outfile,'/'];
+    system(['mkdir -p ',BASIC.MISCDIR]);
     CMD = ['cp ', BASIC.RESDIR,'outputs* ',BASIC.MISCDIR]; disp(CMD);
     system(CMD);
 else% Marconi results
@@ -17,16 +19,14 @@ outfile ='';
 outfile ='';
 outfile ='';
 outfile ='';
-outfile ='';
-% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-01_SGGK_mu_0e+00/out.txt';
-outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-01_SGGK_mu_3e-02/out.txt';
 BASIC.RESDIR      = ['../',outfile(46:end-8),'/'];
 BASIC.MISCDIR     = ['/misc/HeLaZ_outputs/',outfile(46:end-8),'/'];
 end
 
 %% Load the results
 % Load outputs from jobnummin up to jobnummax
-JOBNUMMIN = 00; JOBNUMMAX = 20; 
+JOBNUMMIN = 06; JOBNUMMAX = 20; 
+% JOBNUMMIN = 07; JOBNUMMAX = 20; % For CO damping sim A 
 compile_results %Compile the results from first output found to JOBNUMMAX if existing
 
 %% Post-processing
@@ -44,13 +44,13 @@ end
 
 if 1
 %% Space time diagramm (fig 11 Ivanov 2020)
-TAVG = 1000; % Averaging time duration
+TAVG_0 = 1.4e4; TAVG_1 = 1.5e4; % Averaging times duration
 plot_radial_transport_and_shear
 end
 
 if 0
 %% Space time diagramms
-cmax = 0.01 % max of the colorbar for transport
+cmax = 0.0001 % max of the colorbar for transport
 tstart = 0; tend = Ts3D(end); % time window
 plot_space_time_diagrams
 end
@@ -58,21 +58,22 @@ end
 if 0
 %% |phi_k|^2 spectra (Kobayashi 2015 fig 3)
 % tstart = 0.8*Ts3D(end); tend = Ts3D(end); % Time window
-tstart = 5000; tend = tstart+1000;
+tstart = 14000; tend = 15000;
+% tstart = 10000; tend = 12000;
 % Chose the field to plot
 % FIELD = Ni00;   FNAME = 'Ni00'; FIELDLTX = 'N_i^{00}';
-% FIELD = Ne00;   FNAME = 'Ne00'; FIELDLTX = 'N_e^{00}'
+FIELD = Ne00;   FNAME = 'Ne00'; FIELDLTX = 'N_e^{00}'
 % FIELD = PHI; FNAME = 'PHI'; FIELDLTX = '\tilde\phi';
-FIELD_ = fft2(Gamma_x); FIELD = FIELD_(1:76,:,:,:); FNAME = 'Gamma_x'; FIELDLTX = '\tilde\Gamma_x';
-LOGSCALE = 1; TRENDS = 0;
+% FIELD_ = fft2(Gamma_x); FIELD = FIELD_(1:76,:,:,:); FNAME = 'Gamma_x'; FIELDLTX = '\tilde\Gamma_x';
+LOGSCALE = 0; TRENDS = 0; NORMALIZED = 0;
 plot_kperp_spectrum
 end
 
 if 0
 %% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Options
-t0    =3000; iz = 1; ix = 1; iy = 1;
-skip_ =1; DELAY = 1e-2*skip_;
+t0    =10950; iz = 1; ix = 1; iy = 1;
+skip_ =1; DELAY = 2e-3*skip_;
 [~, it03D] = min(abs(Ts3D-t0)); FRAMES_3D = it03D:skip_:numel(Ts3D);
 [~, it05D] = min(abs(Ts5D-t0)); FRAMES_5D = it05D:skip_:numel(Ts5D);
 INTERP = 0; T = Ts3D; FRAMES = FRAMES_3D;
@@ -118,17 +119,17 @@ if 0
 % FIELD = ne00;          FNAME = 'ne00';    FIELDLTX = 'n_e^{00}'
 % FIELD = dens_i;        FNAME = 'ni';      FIELDLTX = 'n_i';
 % FIELD = dens_e;        FNAME = 'ne';      FIELDLTX = 'n_e';
-FIELD = dens_e-Z_n_e;   FNAME = 'ne_NZ';   FIELDLTX = 'n_e^{NZ}';
+% FIELD = dens_e-Z_n_e;   FNAME = 'ne_NZ';   FIELDLTX = 'n_e^{NZ}';
 % FIELD = dens_i-Z_n_i;   FNAME = 'ni_NZ';   FIELDLTX = 'n_i^{NZ}';
 % FIELD = temp_i;        FNAME = 'Ti';      FIELDLTX = 'T_i';
 % FIELD = temp_e;        FNAME = 'Te';      FIELDLTX = 'T_e';
-% FIELD = phi;           FNAME = 'phi';     FIELDLTX = '\phi';
+FIELD = phi;           FNAME = 'phi';     FIELDLTX = '\phi';
 % FIELD = Z_phi-phi;     FNAME = 'phi_NZ';  FIELDLTX = '\phi^{NZ}';
 % FIELD = Gamma_x;       FNAME = 'Gamma_x'; FIELDLTX = '\Gamma_x';
 % FIELD = dens_e-Z_n_e-(Z_phi-phi);       FNAME = 'Non_adiab_part'; FIELDLTX = 'n_e^{NZ}-\phi^{NZ}';
 
 % Chose when to plot it
-tf = 500:500:2500;
+tf = 11000:50:11300;
 
 % Sliced
 ix = 1; iy = 1; iz = 1;
@@ -145,13 +146,13 @@ plt = @(x,it) real(x( :, :,iz,it)); X = X_XY; Y = Y_XY; XNAME = 'x'; YNAME = 'y'
 %
 TNAME = [];
 fig = figure; FIGNAME = [FNAME,TNAME,'_snaps','_',PARAMS]; set(gcf, 'Position',  [100, 100, 1500, 350])
-plt_2 = @(x) x./max(max(x));
+plt_2 = @(x) x;%./max(max(x));
     for i_ = 1:numel(tf)
     [~,it] = min(abs(Ts3D-tf(i_))); TNAME = [TNAME,'_',num2str(Ts3D(it))];
     subplot(1,numel(tf),i_)
         DATA = plt_2(squeeze(plt(FIELD,it)));
         pclr = pcolor((X),(Y),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
-        colormap(bluewhitered); caxis([-1,1]);
+        colormap(bluewhitered); caxis([-30,30]);
         xlabel(latexize(XNAME)); ylabel(latexize(YNAME));set(gca,'ytick',[]); 
         title(sprintf('$t c_s/R=%.0f$',Ts3D(it)));
     end
@@ -163,14 +164,14 @@ if 0
 %% Photomaton : k space
 
 % Chose the field to plot
-% FIELD = Ni00;   FNAME = 'Ni00'; FIELDLTX = 'N_i^{00}';
-FIELD = Ne00;   FNAME = 'Ne00'; FIELDLTX = 'N_e^{00}'
+FIELD = Ni00;   FNAME = 'Ni00'; FIELDLTX = 'N_i^{00}';
+% FIELD = Ne00;   FNAME = 'Ne00'; FIELDLTX = 'N_e^{00}'
 % FIELD = PHI; FNAME = 'PHI'; FIELDLTX = '\tilde\phi';
 % FIELD_ = fft2(Gamma_x); FIELD = FIELD_(1:Nx/2+1,:,:,:); FNAME = 'Gamma_x'; FIELDLTX = '\tilde\Gamma_x';
 % FIELD_ = fft2(dens_e); FIELD = FIELD_(1:Nx/2+1,:,:,:); FNAME = 'FFT_Dens_e'; FIELDLTX = '\tilde n_e';
 
 % Chose when to plot it
-tf = 500:500:2500;
+tf = 14000:50:14200;
 
 % Sliced
 ix = 1; iy = 1; iz = 1;
@@ -196,5 +197,47 @@ plt_2 = @(x) (fftshift(x,2));
 save_figure
 end
 
+if 0
+%%
+TAVG_0 = 10000; TAVG_1 = 15000; % Averaging times duration
+ZF_fourier_analysis
+end
+if 1
 %%
-% ZF_fourier_analysis
\ No newline at end of file
+plot_param_evol
+end
+
+if 0
+%%
+figure
+plot(Ts3D,shear_maxx_avgy);
+
+end
+
+if 0
+%% zonal vs nonzonal energies for phi(t)
+Ephi_Z           = zeros(1,Ns3D);
+Ephi_NZ          = zeros(1,Ns3D);    
+for it = 1:numel(Ts3D)
+    Ephi_NZ(it) = sum(sum(((KY~=0).*abs(PHI(:,:,1,it)).^2)));
+    Ephi_Z(it)  = sum(sum(((KY==0).*abs(PHI(:,:,1,it)).^2)));
+end
+pltx = @(x) x-x(1);
+plty = @(x) x./max(x);
+fig = figure; FIGNAME = ['ZF_turb_energy_vs_time_',PARAMS];
+set(gcf, 'Position',  [100, 100, 1400, 500])
+subplot(131)
+    semilogy(pltx(Ts3D),plty(Ephi_Z),'DisplayName',['$|\phi_k|^2$ ',CONAME]); hold on;
+    title('Zonal Energy'); legend('show')
+    xlabel('$t c_s/R$'); grid on;% xlim([0 500]);
+
+subplot(132)
+    semilogy(pltx(Ts3D),plty(Ephi_NZ));
+    title('Non Zonal Energy'); legend(CONAME)
+    xlabel('$t c_s/R$'); grid on;% xlim([0 500]);
+    
+subplot(133)
+    semilogy(pltx(Ts0D),plty(abs(PGAMMA_RI)*SCALE));
+    title('Radial particle flux'); legend(CONAME)
+    xlabel('$t c_s/R$'); grid on;% xlim([0 500]);
+end
\ No newline at end of file
diff --git a/wk/continue_multiple_runs_marconi.m b/wk/continue_multiple_runs_marconi.m
index c3d7e300..97346616 100644
--- a/wk/continue_multiple_runs_marconi.m
+++ b/wk/continue_multiple_runs_marconi.m
@@ -1,12 +1,10 @@
 %% Paste the list of continue_run calls
-continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_5e-01_SGGK_mu_3e-02/out.txt')
-continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-01_SGGK_mu_3e-02/out.txt')
-continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_5e-01_SGGK_mu_0e+00/out.txt')
-continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-01_SGGK_mu_0e+00/out.txt')
+
+continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_7e-02_SGGK_mu_0e+00/out.txt')
 
 %% Functions to modify preexisting fort.90 input file and launch on marconi
 function [] = continue_run(outfilename)
-    EXECNAME = 'helaz_3.1';
+    EXECNAME = 'helaz_3.2';
     %% CLUSTER PARAMETERS
     CLUSTER.PART  = 'prod';     % dbg or prod
     CLUSTER.TIME  = '24:00:00'; % allocation time hh:mm:ss
@@ -50,12 +48,16 @@ function [] = continue_run(outfilename)
         J2L = str2num(line) + 1;
     end
     % Change job 2 load in fort.90
-    A{39} = ['  job2load      = ',num2str(3)];
+    A{39} = ['  job2load      = ',num2str(J2L)];
     disp(A{39})
     % Change time step
-    A{3} = ['  dt     = 0.01'];
+    A{3} = ['  dt     = 0.005'];
     % Increase endtime
-    A{4} = ['  tmax      = 10000'];
+    A{4} = ['  tmax      = 20000'];
+    % Change collision operator
+    line_= A{43};
+    CO_old = str2num(line_(13:end));
+    A{43} = ['  CO      = ',num2str(1)];
     % Put non linear term back
     A{45} = ['  NL_CLOS = -1'];
     % change HD
@@ -89,4 +91,4 @@ function [] = continue_run(outfilename)
     write_sbash_marconi
     % Launch the job
     system('ssh ahoffman@login.marconi.cineca.it sh HeLaZ/wk/setup_and_run.sh');
-end
\ No newline at end of file
+end
diff --git a/wk/linear_study.m b/wk/linear_study.m
index efbbcce1..3f174a96 100644
--- a/wk/linear_study.m
+++ b/wk/linear_study.m
@@ -1,4 +1,5 @@
-%clear all;
+for CO = [1]
+    RUN = 1; % To run or just to load
 addpath(genpath('../matlab')) % ... add
 default_plots_options
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -6,45 +7,49 @@ default_plots_options
 CLUSTER.TIME  = '99:00:00'; % allocation time hh:mm:ss
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% PHYSICAL PARAMETERS
-% NU      = 1.0;   % Collision frequency
+NU      = 1.0;   % Collision frequency
 TAU     = 1.0;    % e/i temperature ratio
 ETAB    = 1.0;
-% ETAN    = ETAN;    % Density gradient
+ETAN    = 1/0.6;    % Density gradient
 ETAT    = 0.0;    % Temperature gradient
 NU_HYP  = 0.0;   % Hyperdiffusivity coefficient
 LAMBDAD = 0.0;
-NOISE0  = 1.0e-5;
+NOISE0  = 1.0e-5; % Init noise amplitude
+BCKGD0  = 0.0;    % Init background
 %% GRID PARAMETERS
-N       = 50;     % Frequency gridpoints (Nkx = N/2)
-L       = 75;     % Size of the squared frequency domain
+N       = 150;     % Frequency gridpoints (Nkx = N/2)
+L       = 100;     % Size of the squared frequency domain
 KXEQ0   = 1;      % put kx = 0
 MU_P    = 0.0;     % Hermite  hyperdiffusivity -mu_p*(d/dvpar)^4 f
 MU_J    = 0.0;     % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
 %% TIME PARMETERS
-TMAX    = 200;  % Maximal time unit
-DT      = 1e-2;   % Time step
+TMAX    = 400;  % Maximal time unit
+DT      = 8e-3;   % Time step
 SPS0D   = 1;      % Sampling per time unit for 2D arrays
-SPS2D   = 1;      % Sampling per time unit for 2D arrays
-SPS5D   = 1/5;    % Sampling per time unit for 5D arrays
+SPS2D   = 0;      % Sampling per time unit for 2D arrays
+SPS3D   = 1;      % Sampling per time unit for 2D arrays
+SPS5D   = 1/50;    % Sampling per time unit for 5D arrays
 SPSCP   = 0;    % Sampling per time unit for checkpoints
 RESTART = 0;      % To restart from last checkpoint
 JOB2LOAD= 00;
 %% OPTIONS
-SIMID   = 'v3.1_lin_analysis';  % Name of the simulation
+% SIMID   = 'v3.6_kobayashi_lin';  % Name of the simulation
+% SIMID   = 'v3.2_CO_damping';  % Name of the simulation
+% SIMID   = 'CO_Patchwork_damping';  % Name of the simulation
+SIMID   = 'v3.4_entropy_mode_linear';  % Name of the simulation
 NON_LIN = 0 *(1-KXEQ0);   % activate non-linearity (is cancelled if KXEQ0 = 1)
 % Collision operator
 % (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Pitch angle, 4 : Full Couloumb ; +/- for GK/DK)
-% CO      = 2;
+% CO      = 1;
 INIT_ZF = 0; ZF_AMP = 0.0;
 CLOS    = 0;   % Closure model (0: =0 truncation, 1: semi coll, 2: Copy closure J+1 = J, P+2 = P)
 NL_CLOS = 0;   % nonlinear closure model (0: =0 nmax = jmax, 1: nmax = jmax-j, >1 : nmax = NL_CLOS)
 KERN    = 0;   % Kernel model (0 : GK)
 INIT_PHI= 0;   % Start simulation with a noisy phi
-INIT_ZF = 0;   % Start simulation with a noisy phi
 %% OUTPUTS
 W_DOUBLE = 0;
 W_GAMMA  = 0;
-W_PHI    = 0;
+W_PHI    = 1;
 W_NA00   = 1;
 W_NAPJ   = 1;
 W_SAPJ   = 0;
@@ -63,14 +68,15 @@ Nz      = 1;      % number of perpendicular planes (parallel grid)
 q0      = 1.0;    % safety factor
 shear   = 0.0;    % magnetic shear
 eps     = 0.0;    % inverse aspect ratio
+INIT_BLOB = 0; WIPE_TURB = 0; WIPE_ZF = 0;
 %% PARAMETER SCANS
 
-if 0
+if 1
 %% Parameter scan over PJ
-% PA = [2 4 6 10];
-% JA = [1 2 3  5];
-PA = [6];
-JA = [3];
+% PA = [2 4];
+% JA = [1 2];
+PA = [4];
+JA = [2];
 DTA= DT*ones(size(JA));%./sqrt(JA);
 % DTA= DT;
 mup_ = MU_P;
@@ -79,8 +85,9 @@ Nparam = numel(PA);
 param_name = 'PJ';
 gamma_Ni00 = zeros(Nparam,floor(N/2)+1);
 gamma_Nipj = zeros(Nparam,floor(N/2)+1);
-Bohm_transport = zeros(Nparam,1);
-Ni00_ST  = zeros(Nparam,floor(N/2)+1,floor(SPS2D*TMAX));
+gamma_phi  = zeros(Nparam,floor(N/2)+1);
+Ni00_ST  = zeros(Nparam,floor(N/2)+1,TMAX/SPS3D);
+ PHI_ST  = zeros(Nparam,floor(N/2)+1,TMAX/SPS3D);
 for i = 1:Nparam
     % Change scan parameter
     PMAXE = PA(i); PMAXI = PA(i);
@@ -88,32 +95,35 @@ for i = 1:Nparam
     DT = DTA(i);
     setup
     % Run linear simulation
-%     system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ./../../../bin/helaz_3 1 1; cd ../../../wk'])
-    system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz_3 1 6; cd ../../../wk'])
-%     system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz_3 2 3; cd ../../../wk'])
+    if RUN
+        system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz_3.6 1 6; cd ../../../wk'])
+    end
 %     Load and process results
     %%
     filename = ['../results/',SIMID,'/',PARAMS,'/outputs_00.h5'];
     load_results
-    tend   = Ts3D(end); tstart   = 0.4*tend;
-    [~,itstart] = min(abs(Ts3D-tstart));
-    [~,itend]   = min(abs(Ts3D-tend));
     for ikx = 1:N/2+1
-        gamma_Ni00(i,ikx) = (LinearFit_s(Ts3D(itstart:itend)',(squeeze(abs(Ni00(ikx,1,itstart:itend))))));
-        Ni00_ST(i,ikx,1:numel(Ts3D)) = squeeze((Ni00(ikx,1,:)));
+        %find if there is a tail of 0 (highest damped mode)
+
+        tend   = max(Ts3D(abs(Ni00(ikx,1,1,:))~=0));
+        tstart   = 0.8*tend;
+        [~,itstart] = min(abs(Ts3D-tstart));
+        [~,itend]   = min(abs(Ts3D-tend));
+        gamma_Ni00(i,ikx) = (LinearFit_s(Ts3D(itstart:itend)',(squeeze(abs(Ni00(ikx,1,1,itstart:itend))))));
+        gamma_phi (i,ikx) = (LinearFit_s(Ts3D(itstart:itend)',(squeeze(abs(PHI (ikx,1,1,itstart:itend))))));
+        Ni00_ST(i,ikx,:) = squeeze((Ni00(ikx,1,1,1:TMAX/SPS3D)));
+         PHI_ST(i,ikx,:) = squeeze((PHI (ikx,1,1,1:TMAX/SPS3D)));
     end
     tend   = Ts5D(end); tstart   = 0.4*tend;
     [~,itstart] = min(abs(Ts5D-tstart));
     [~,itend]   = min(abs(Ts5D-tend));
     for ikx = 1:N/2+1
-        gamma_Nipj(i,ikx) = LinearFit_s(Ts5D(itstart:itend)',squeeze(max(max(abs(Nipj(:,:,ikx,1,itstart:itend)),[],1),[],2)));
+        gamma_Nipj(i,ikx) = LinearFit_s(Ts5D(itstart:itend)',squeeze(max(max(abs(Nipj(:,:,ikx,1,1,itstart:itend)),[],1),[],2)));
     end
-    gamma_Ni00(i,:) = real(gamma_Ni00(i,:) .* (gamma_Ni00(i,:)>=0.0));
-    gamma_Nipj(i,:) = real(gamma_Nipj(i,:) .* (gamma_Nipj(i,:)>=0.0));
-%     kymax = abs(kx(ikymax));
-%     Bohm_transport(i) = ETAB/ETAN*gmax/kymax^2;
+    gamma_Ni00(i,:) = real(gamma_Ni00(i,:));% .* (gamma_Ni00(i,:)>=0.0));
+    gamma_Nipj(i,:) = real(gamma_Nipj(i,:));% .* (gamma_Nipj(i,:)>=0.0));
     % Clean output
-    system(['rm -r ',BASIC.RESDIR]);
+%     system(['rm -r ',BASIC.RESDIR]);
 end
 
 if 1
@@ -125,102 +135,27 @@ plt = @(x) x;
     for i = 1:Nparam
         clr       = line_colors(mod(i-1,numel(line_colors(:,1)))+1,:);
         linestyle = line_styles(floor((i-1)/numel(line_colors(:,1)))+1);
-        semilogx(plt(SCALE*kx(2:numel(kx))),plt(gamma_Ni00(i,2:end)),...
+        semilogx(plt(SCALE*kx(1:numel(kx))),plt(gamma_Ni00(i,1:end)),...
             'Color',clr,...
             'LineStyle',linestyle{1},'Marker','^',...
-            'DisplayName',['$\eta=',num2str(ETAB/ETAN),'$, $\nu_{',CONAME,'}=',num2str(NU),'$, $P=',num2str(PA(i)),'$, $J=',num2str(JA(i)),'$']);
+...%             'DisplayName',['$\eta=',num2str(ETAB/ETAN),'$, $\nu_{',CONAME,'}=',num2str(NU),'$, $P=',num2str(PA(i)),'$, $J=',num2str(JA(i)),'$']);
+            'DisplayName',[CONAME,', $P,J=',num2str(PA(i)),',',num2str(JA(i)),'$']);
         hold on;
     end
     grid on; xlabel('$k_z\rho_s^{R}$'); ylabel('$\gamma(N_i^{00})L_\perp/c_s$'); xlim([0.0,max(kx)]);
     title(['$\eta=',num2str(ETAB/ETAN),'$, $\nu_{',CONAME,'}=',num2str(NU),'$'])
-    legend('show'); xlim([0.01,10])
+%     title(['$\nabla N = 0$', ', $\nu=',num2str(NU),'$'])
+    legend('show'); %xlim([0.01,10])
 saveas(fig,[SIMDIR,'gamma_Ni_vs_',param_name,'_',PARAMS,'.fig']);
 saveas(fig,[SIMDIR,'gamma_Ni_vs_',param_name,'_',PARAMS,'.png']);
 end
 end
-
-if 1
-%% Parameter scan over CO
-P=2; J=1;
-N       = 20;     % Frequency gridpoints (Nkx = N/2)
-L       = 75;     % Size of the squared frequency domain
-TMAX  = 200;
-DT    = 0.01;
-CO_A = [4];
-CONAME_A = {};
-Nparam = numel(CO_A);
-param_name = 'CO';
-ETAN    = 2.0;
-NU      = 1e-1;   % Collision frequency
-PMAXE = P; PMAXI = P;
-JMAXE = J; JMAXI = J;
-Bohm_transport = zeros(Nparam,1);
-gamma_Ni00 = zeros(Nparam,floor(N/2)+1);
-gamma_Nipj = zeros(Nparam,floor(N/2)+1);
-
-for i = 1:Nparam
-    % Change scan parameter
-    CO = CO_A(i);
-    setup
-    CONAME_A{i} = CONAME;
-    % Run linear simulation
-    system(...
-        ['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz_3 1 6; cd ../../../wk']...
-    )
-    %% Load an process results
-    filename = ['../results/',SIMID,'/',PARAMS,'/outputs_00.h5'];
-    load_results
-    tend   = Ts3D(end); tstart   = 0.4*tend;
-    [~,itstart] = min(abs(Ts3D-tstart));
-    [~,itend]   = min(abs(Ts3D-tend));
-    for ikx = 1:N/2+1
-        gamma_Ni00(i,ikx) = (LinearFit_s(Ts3D(itstart:itend)',(squeeze(abs(Ni00(ikx,1,itstart:itend))))));
-        Ni00_ST(i,ikx,1:numel(Ts3D)) = squeeze((Ni00(ikx,1,:)));
-    end
-    tend   = Ts5D(end); tstart   = 0.4*tend;
-    [~,itstart] = min(abs(Ts5D-tstart));
-    [~,itend]   = min(abs(Ts5D-tend));
-    for ikx = 1:N/2+1
-        gamma_Nipj(i,ikx) = LinearFit_s(Ts5D(itstart:itend)',squeeze(max(max(abs(Nipj(:,:,ikx,1,itstart:itend)),[],1),[],2)));
-    end
-    gamma_Ni00(i,:) = real(gamma_Ni00(i,:) .* (gamma_Ni00(i,:)>=0.0));
-    gamma_Nipj(i,:) = real(gamma_Nipj(i,:) .* (gamma_Nipj(i,:)>=0.0));
-%     kymax = abs(kx(ikymax));
-%     Bohm_transport(i) = ETAB/ETAN*gmax/kymax^2;
-    % Clean output
-    system(['rm -r ',BASIC.RESDIR]);
-end
-
-if 1
-%% Plot
-SCALE = 1;%sqrt(2);
-fig = figure; FIGNAME = 'linear_study';
-plt = @(x) x;
-% subplot(211)
-    for i = 1:Nparam
-        clr       = line_colors(mod(i-1,numel(line_colors(:,1)))+1,:);
-        linestyle = line_styles(floor((i-1)/numel(line_colors(:,1)))+1);
-        semilogx(plt(SCALE*kx(2:numel(kx))),plt(gamma_Ni00(i,2:end)),...
-            'Color',clr,...
-            'LineStyle',linestyle{1},'Marker','^',...
-            'DisplayName',[CONAME_A{i}]);
-        hold on;
-    end
-    grid on; xlabel('$k_z\rho_s^{R}$'); ylabel('$\gamma(N_i^{00})L_\perp/c_s$'); xlim([0.0,max(kx)]);
-    title(['$\eta=',num2str(ETAB/ETAN),'$, $\nu=',num2str(NU),'$',', $P=',num2str(PMAXI),'$, $J=',num2str(JMAXI),'$'])
-    legend('show'); xlim([0.01,10])
-saveas(fig,[SIMDIR,'gamma_Ni_vs_',param_name,'_',PARAMS,'.fig']);
-saveas(fig,[SIMDIR,'gamma_Ni_vs_',param_name,'_',PARAMS,'.png']);
-end
-
 if 0
-%% Plot
-fig = figure; FIGNAME = 'mixing_length';
-plot(eta_B, Bohm_transport)
-grid on; xlabel('$L_n/L_B$'); ylabel('$\eta\gamma_{max}/k_{max}^2$');
-title(['$P_e=',num2str(PMAXE),'$',', $J_e=',num2str(JMAXE),'$',...
-       ', $P_i=',num2str(PMAXE),'$',', $J_i=',num2str(JMAXI),'$'])
-saveas(fig,[SIMDIR,FIGNAME,'_vs_',param_name,'_',PARAMS,'.fig']);
+%% Space time
+    [YT,XT] = meshgrid(Ts3D,kx);
+    figure;
+%     pclr = surf(XT,YT,squeeze(abs(PHI_ST(1,:,:)))); set(pclr, 'edgecolor','none'); colorbar;
+%     pclr = pcolor(XT,YT,squeeze(abs(Ni00_ST(1,:,:)))); set(pclr, 'edgecolor','none'); colorbar;
+    semilogy(Ts3D(1:TMAX/SPS3D),squeeze(abs(PHI_ST(1,50:5:100,:))));
+end
 end
-%%
-end
\ No newline at end of file
diff --git a/wk/load_multiple_outputs_marconi.m b/wk/load_multiple_outputs_marconi.m
index 31afde74..4d2acbcb 100644
--- a/wk/load_multiple_outputs_marconi.m
+++ b/wk/load_multiple_outputs_marconi.m
@@ -4,3 +4,7 @@ load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HD_study/150x
 load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-01_SGGK_mu_3e-02/out.txt')
 load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_5e-01_SGGK_mu_0e+00/out.txt')
 load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-01_SGGK_mu_0e+00/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_2e-01_SGGK_mu_0e+00/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_7e-02_SGGK_mu_0e+00/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-02_SGGK_mu_3e-02/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e+00_SGGK_mu_0e+00/out.txt')
diff --git a/wk/local_run.m b/wk/local_run.m
index 5bb94669..542b5091 100644
--- a/wk/local_run.m
+++ b/wk/local_run.m
@@ -38,7 +38,7 @@ SIMID   = 'HD_study';  % Name of the simulation
 NON_LIN = -1;   % activate non-linearity (is cancelled if KXEQ0 = 1)
 % INIT options
 INIT_ZF = 0; ZF_AMP = 0.0;
-INIT_BLOB = 0; WIPE_TURB = 0;
+INIT_BLOB = 0; WIPE_TURB = 0; WIPE_ZF = 0;
 %% OUTPUTS
 W_DOUBLE = 0;
 W_GAMMA  = 1;
diff --git a/wk/marconi_run.m b/wk/marconi_run.m
index 1c4f2ac0..06111db2 100644
--- a/wk/marconi_run.m
+++ b/wk/marconi_run.m
@@ -2,7 +2,7 @@ clear all;
 addpath(genpath('../matlab')) % ... add
 SUBMIT = 1; % To submit the job automatically
 % EXECNAME = 'helaz_dbg';
-  EXECNAME = 'helaz_3.1';
+  EXECNAME = 'helaz_3.2';
 for ETAN = [1/0.6]
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% Set Up parameters
@@ -14,13 +14,13 @@ CLUSTER.TIME  = '24:00:00'; % allocation time hh:mm:ss
 if(strcmp(CLUSTER.PART,'dbg')); CLUSTER.TIME  = '00:30:00'; end;
 CLUSTER.MEM   = '128GB';     % Memory
 CLUSTER.JNAME = 'HeLaZ';% Job name
-NP_P          = 2;          % MPI processes along p  
+NP_P          = 2;          % MPI processes along p
 NP_KX         = 24;         % MPI processes along kx
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% PHYSICAL PARAMETERS
-NU      = 0.01;   % Collision frequency
+NU      = 1.0;   % Collision frequency
 ETAN    = 1.0/0.6;    % Density gradient drive (R/Ln)
-NU_HYP  = 1.0;
+NU_HYP  = 0.0;
 %% GRID PARAMETERS
 N       = 150;     % Frequency gridpoints (Nkx = N/2)
 L       = 100;     % Size of the squared frequency domain
@@ -32,7 +32,7 @@ P       = 4;
 J       = 2;
 %% TIME PARAMETERS
 TMAX    = 10000;  % Maximal time unit
-DT      = 1e-2;   % Time step
+DT      = 5e-3;   % Time step
 SPS0D   = 1;      % Sampling per time unit for profiler
 SPS2D   = 1;      % Sampling per time unit for 2D arrays
 SPS3D   = 1/2;      % Sampling per time unit for 3D arrays
@@ -52,7 +52,7 @@ SIMID   = 'HD_study';  % Name of the simulation
 NON_LIN = 1;   % activate non-linearity (is cancelled if KXEQ0 = 1)
 % INIT options
 INIT_ZF = 0; ZF_AMP = 0.0;
-INIT_BLOB = 0; WIPE_TURB = 0;
+INIT_BLOB = 0; WIPE_TURB = 0; WIPE_ZF = 0;
 %% OUTPUTS
 W_DOUBLE = 1;
 W_GAMMA  = 1;
@@ -88,7 +88,7 @@ MU_J    = 0.0;     % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
 % Compute processes distribution
 Ntot = NP_P * NP_KX;
 Nnodes = ceil(Ntot/48);
-Nppn   = Ntot/Nnodes; 
+Nppn   = Ntot/Nnodes;
 CLUSTER.NODES =  num2str(Nnodes);  % MPI process along p
 CLUSTER.NTPN  =  num2str(Nppn); % MPI process along kx
 CLUSTER.CPUPT = '1';        % CPU per task
@@ -103,4 +103,4 @@ if(SUBMIT)
     system('ssh ahoffman@login.marconi.cineca.it sh HeLaZ/wk/setup_and_run.sh');
 end
 disp('done');
-end
\ No newline at end of file
+end
diff --git a/wk/new_flux_results.m b/wk/new_flux_results.m
new file mode 100644
index 00000000..c52bec38
--- /dev/null
+++ b/wk/new_flux_results.m
@@ -0,0 +1,31 @@
+%% eta = 0.6
+%nu Gammainf mu
+SGGK_transport = [...
+    1.0e+0, 5.6e-1, 0.0e+0;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e+00_SGGK_mu_0e+00/
+    5.0e-1, 4.0e-1, 0.0e+0;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_5e-01_SGGK_mu_0e+00/ before wiping ZF
+    5.0e-1, 5.2e-1, 0.0e+0;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_5e-01_SGGK_mu_0e+00/ after wiping ZF
+    2.5e-1, 3.6e-1, 0.0e+0;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_2e-01_SGGK_mu_0e+00/
+    1.0e-1, 2.2e-1, 0.0e+0;... 
+    1.0e-1, 1.9e-1, 0.0e+0;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-01_SGGK_mu_0e+00
+    7.5e-2, 1.5e-1, 0.0e+0;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_7e-02_SGGK_mu_0e+00/
+    5.0e-2, 1.1e-1, 0.0e+0;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_7e-02_SGGK_mu_0e+00/
+    1.0e-2, 4.6e-2, 3.2e-2;...
+    ];
+
+DGGK_transport = [...
+    1.0e+0, 5.7e+00, 0.0e+0;... % HeLaZ 2.8 P,J=10,5
+    5.0e-1, 1.1e+01, 0.0e+0;... % simulation_B/cw_DGGK
+    2.5e-1, 6.4e+00, 0.0e+0;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_2e-01_SGGK_mu_0e+00
+    1.0e-1, 3.0e+00, 0.0e+0;... % HeLaZ 2.8 P,J=2,1
+    1.0e-1, 2.5e+00, 0.0e+0;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-01_SGGK_mu_0e+00
+    7.5e-2, 1.8e+00, 0.0e+0;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_7e-02_SGGK_mu_0e+00
+    1.0e-2, 3.3e-01, 0.0e+0;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-02_DGGK_mu_5e-04
+    1.0e-2, 2.2e-01, 1.0e-4;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-02_DGGK_mu_3e-03
+    1.0e-2, 1.4e-01, 3.0e-3;... % HeLaZ 2.8 P,J=6,3
+    ];
+
+figure;
+semilogy(SGGK_transport(:,1),SGGK_transport(:,2),'.','MarkerSize',30); hold on;
+semilogy(DGGK_transport(:,1),DGGK_transport(:,2),'.','MarkerSize',30);
+grid on; xlabel('$\nu$'); ylabel('$\Gamma_x$');
+legend(['SGGK';'DGGK'])
\ No newline at end of file
diff --git a/wk/open_figure_script.m b/wk/open_figure_script.m
index cc22dcf6..48416378 100644
--- a/wk/open_figure_script.m
+++ b/wk/open_figure_script.m
@@ -15,13 +15,13 @@ simname_ = fname_(54:end-8);
 % simname_ = '';
 % simname_ = '';
 % simname_ = '';
-% simname_ = '';
+simname_ = 'simulation_A/DGGK_damping_150x75_L_100_P_4_J_2_eta_0.6_nu_1e-01_SGGK_mu_0e+00';
 
 
 
 
 %%
-figname_ = '/fig/ZF_transport_drphi_';
+figname_ = '/fig/ZF_transport_drphi_*';
 % figname_ = '/fig/space_time_';
 % figname_ = '/fig/phi_shear_phase_space_';
 
@@ -31,5 +31,5 @@ path_    = '../results/';
 
 params_  = simname_(idx-3:end);
 
-
-openfig([path_,simname_,figname_,params_,'.fig']);
\ No newline at end of file
+[~,a] =  system(['ls ',[path_,simname_,figname_,'.fig']]);
+openfig(a(1:end-1))
diff --git a/wk/plot_cosol_mat.m b/wk/plot_cosol_mat.m
index d320b533..82a93066 100644
--- a/wk/plot_cosol_mat.m
+++ b/wk/plot_cosol_mat.m
@@ -94,7 +94,10 @@ subplot(224)
 
 % mat_file_name = '/home/ahoffman/HeLaZ/iCa/gk_coulomb_P_2_J_1_N_1_kpm_4.0_NFLR_5.h5';
 % mat_file_name = '/home/ahoffman/HeLaZ/iCa/gk_coulomb_P_6_J_3_N_150_kpm_8.0_NFLR_8.h5';
-mat_file_name = '/home/ahoffman/cosolver/gk.coulomb.NFLR/k.4/self.NFLR.8.0.h5';
+% mat_file_name = '/home/ahoffman/cosolver/gk.coulomb.NFLR/k.4/self.NFLR.8.0.h5';
+% mat_file_name = '/home/ahoffman/cosolver/gk.coulomb_6_3_NFLR_8_kp_4_npiflr_8_sigma_3e-3/scanfiles_00000/self.0.h5';
+% mat_file_name = '/home/ahoffman/cosolver/gk.coulomb_6_3_NFLR_8_kp_4_npiflr_0_sigma_5e-4/scanfiles_00000/self.0.h5';
+mat_file_name = '/home/ahoffman/cosolver/gk.coulomb_6_3_NFLR_8_kp_4_bjf/scanfiles_00000/self.0.h5';
 matidx = 74;
 
 matidx = sprintf('%5.5i',matidx);disp(matidx);
-- 
GitLab