diff --git a/wk/HD_study.m b/wk/HD_study.m
new file mode 100644
index 0000000000000000000000000000000000000000..f7769719529079480088e8127e82e3cfc549c739
--- /dev/null
+++ b/wk/HD_study.m
@@ -0,0 +1,48 @@
+% bursts simulations
+b_=[...
+    1e-1, 1e-1;...
+    1e-1, 3e-2;...
+    1e-1, 2e-2;...
+    1e-2, 3e-3;...
+    5e-2, 3e-3;...
+    1e-3, 2e-2;...
+    ];
+% converged plateau simulations
+cp_=[...
+    1e+0, 1e-2;...
+    1e+0, 3e-3;...
+    1e+0, 2e-3;...
+    1e-1, 2e-3;...
+    1e-1, 1e-3;...
+    ];
+% moving plateau
+dp_=[...
+    1e+0, 1e+0;...
+    1e-3, 2e-2;...
+    1e-3, 1e-2;...
+    1e-3, 5e-3;...
+    1e-3, 2.5e-3,...
+    ];
+% resolution issue plateau?
+rp_=[...
+    5e-2, 1e-3;...
+    1e-2, 1e-4;...
+    ];
+figure; set(gcf, 'Position',  [100, 100, 900, 400])
+grid on; xlim([5e-4,5e0]); ylim([5e-5,5e+0]);
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+xlabel('$\nu_{DGGK}$'); ylabel('$\mu_{HD}$'); hold on;
+ 
+scatter( b_(:,1), b_(:,2),'o',...
+    'MarkerFaceColor',[0.6350 0.0780 0.1840],'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,...
+    'DisplayName','Converged Plateau'); 
+scatter(dp_(:,1),dp_(:,2),'v',...
+    'MarkerFaceColor',[0.4660 0.6740 0.1880],'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','Moving Plateau'); 
+plot([1 1],[5e-5 5e-1],'--','Color',[0.4660 0.6740 0.1880]);
diff --git a/wk/analysis_2D.m b/wk/analysis_2D.m
index 4b0112d18644cc27014696ebf84b7a76ff6b013b..7278ebcdcc1f04d155a742cb9d977f3b5e2385d2 100644
--- a/wk/analysis_2D.m
+++ b/wk/analysis_2D.m
@@ -5,23 +5,23 @@ for i_ = 1
 if 1% Local results
 outfile ='';
 outfile ='';
-outfile ='HD_study/150x75_L_100_P_2_J_1_eta_0.6_nu_1e-02_DGGK_CLOS_0_mu_3e-03';
-% outfile ='HD_study/200x100_L_100_P_2_J_1_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_3e-02';
-% outfile ='HD_study/150x75_L_100_P_2_J_1_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_3e-02';
+outfile ='';
+outfile ='';
+% outfile ='HD_study/300x150_L_100_P_2_J_1_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_1e-03';
+% outfile ='HD_study/150x75_L_100_P_2_J_1_eta_0.6_nu_1e-02_DGGK_CLOS_0_mu_3e-03';
+outfile ='HD_study/200x100_L_200_P_2_J_1_eta_0.6_nu_1e+00_DGGK_CLOS_0_mu_0e+00';
     BASIC.RESDIR      = ['../results/',outfile,'/'];
     BASIC.MISCDIR     = ['/misc/HeLaZ_outputs/results/',outfile,'/'];
     CMD = ['cp ', BASIC.RESDIR,'outputs* ',BASIC.MISCDIR]; disp(CMD);
     system(CMD);
-end
-if 0% Marconi results
-outfile ='';
+else% Marconi results
 outfile ='';
 outfile ='';
 outfile ='';
 outfile ='';
 outfile ='';
 outfile ='';
-outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.7_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.6_nu_1e-02_DGGK_CLOS_0_mu_1e-02/out.txt';
+outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.8_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.6_nu_1e-03_SGGK_CLOS_0_mu_2e-02/out.txt';
 % outfile = outcl{i_};
 % load_marconi(outfile);
 BASIC.RESDIR      = ['../',outfile(46:end-8),'/'];
@@ -406,7 +406,7 @@ end
 
 if 1
 %% |phi_k|^2 spectra (Kobayashi 2015 fig 3)
-tstart = 2000; tend = 3000;
+tstart = 10000; tend = 14000;
 [~,itstart] = min(abs(Ts2D-tstart));
 [~,itend]   = min(abs(Ts2D-tend));
 trange = itstart:itend;
@@ -447,10 +447,10 @@ clear Z_rth phi_kp ni_kp Ti_kp
 end
 
 %%
-t0    =3000;
+t0    =800;
 [~, it02D] = min(abs(Ts2D-t0));
 [~, it05D] = min(abs(Ts5D-t0));
-skip_ = 1; 
+skip_ = 5; 
 DELAY = 1e-3*skip_;
 FRAMES_2D = it02D:skip_:numel(Ts2D);
 FRAMES_5D = it05D:skip_:numel(Ts5D);
@@ -492,8 +492,8 @@ if 0
 GIFNAME = ['ni00',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0;
 FIELD = real(ni00); X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D;
 FIELDNAME = '$n_i^{00}$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$';
-create_gif
-% create_mov
+% create_gif
+create_mov
 end
 if 0
 %% GC Density electrons
@@ -589,5 +589,5 @@ FIELDNAME = 'N_i'; XNAME = '$k_z\rho_s$'; YNAME = '$P$, $\sum_r$';
 create_gif_5D
 end
 %%
-ZF_fourier_analysis
+% ZF_fourier_analysis
 end
\ No newline at end of file
diff --git a/wk/continue_multiple_runs_marconi.m b/wk/continue_multiple_runs_marconi.m
index cc2b930440d80c19b0cf337cd5fac5652e318d14..f6f8a27267157d0574a7a2df1ad8b579c094f33f 100644
--- a/wk/continue_multiple_runs_marconi.m
+++ b/wk/continue_multiple_runs_marconi.m
@@ -1,5 +1,5 @@
 %% Paste the list of continue_run calls
-continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.7_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.6_nu_1e-02_DGGK_CLOS_0_mu_1e-02/out.txt')
+continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.8_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.6_nu_1e-03_SGGK_CLOS_0_mu_2e-02/out.txt')
 
 %% Functions to modify preexisting fort.90 input file and launch on marconi
 function [] = continue_run(outfilename)
@@ -51,15 +51,16 @@ function [] = continue_run(outfilename)
     disp(A{35})
     % Change time step
     line_= A{3}; dt_old = str2num(line_(12:end));    
-    A{3} = ['  dt     = ',num2str(dt_old)];
+    A{3} = ['  dt     = ',num2str(1.5*dt_old)];
     % Increase endtime
     A{4} = ['  tmax      = 20000'];
     % Put non linear term back
     line_= A{41}; NL_old = str2num(line_(13:end));    
     A{41} = ['  NL_CLOS = ',num2str(NL_old)];
+%     A{41} = ['  NL_CLOS = -1'];
     % change HD
     line_= A{43}; mu_old = str2num(line_(13:end));
-    A{43} = ['  mu      = ',num2str(mu_old)];
+    A{43} = ['  mu      = ',num2str(0*mu_old)];
     % change N
     line_= A{13}; N_old = str2num(line_(10:end));
     A{13} = ['  Nr = ',num2str(N_old)];
diff --git a/wk/load_multiple_outputs_marconi.m b/wk/load_multiple_outputs_marconi.m
index fc48c1e5bea030b81167d92d47461e9cb9a24028..f647156162f0e6611d726924d46cec265d2d4fb5 100644
--- a/wk/load_multiple_outputs_marconi.m
+++ b/wk/load_multiple_outputs_marconi.m
@@ -1,7 +1,5 @@
 addpath(genpath('../matlab')) % ... add
 %% Paste the list of simulation results to load
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/kobayashi/100x50_L_50_P_6_J_3_eta_0.71429_nu_1e-02_PAGK_CLOS_0_mu_0e+00/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/kobayashi/300x150_L_100_P_10_J_5_eta_0.71429_nu_1e-02_PAGK_CLOS_0_mu_0e+00/out.txt')
 load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.8_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.6_nu_1e-03_SGGK_CLOS_0_mu_2e-02/out.txt')
-load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.7_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.6_nu_1e-02_DGGK_CLOS_0_mu_1e-02/out.txt')
-load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/kobayashi/300x150_L_100_P_6_J_3_eta_0.71429_nu_1e-02_PAGK_CLOS_0_mu_0e+00/out.txt')
-load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/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/out.txt')
-load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.7_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.6_nu_1e-02_SGGK_CLOS_0_mu_1e-02/out.txt')
diff --git a/wk/local_run.m b/wk/local_run.m
index 0e002fc07e52a0b70c4a5e7f1b00da5b93a4ec0a..fa0005c5a6e7c424e98b0bec60a06fea5b72bb97 100644
--- a/wk/local_run.m
+++ b/wk/local_run.m
@@ -4,22 +4,22 @@ addpath(genpath('../matlab')) % ... add
 CLUSTER.TIME  = '99:00:00'; % allocation time hh:mm:ss
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% PHYSICAL PARAMETERS
-NU      = 0.1;   % Collision frequency
+NU      = 1.0;   % Collision frequency
 ETAB    = 0.6;    % Magnetic gradient
 ETAN    = 1.0;    % Density gradient
-NU_HYP  = 0.5;
+NU_HYP  = 0.0;
 %% GRID PARAMETERS
-N       = 300;     % Frequency gridpoints (Nkr = N/2)
+N       = 100;     % Frequency gridpoints (Nkr = N/2)
 L       = 100;     % Size of the squared frequency domain
 P       = 2;
 J       = 1;
 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 PARAMETERS
-TMAX    = 1000;  % Maximal time unit
+TMAX    = 500;  % Maximal time unit
 DT      = 2e-2;   % Time step
-SPS0D   = 1;      % Sampling per time unit for profiler
-SPS2D   = 1;      % Sampling per time unit for 2D arrays
+SPS0D   = 2;      % Sampling per time unit for profiler
+SPS2D   = 2;      % Sampling per time unit for 2D arrays
 SPS5D   = 1/200;  % Sampling per time unit for 5D arrays
 SPSCP   = 0;    % Sampling per time unit for checkpoints/10
 RESTART = 0;      % To restart from last checkpoint
diff --git a/wk/photomaton.m b/wk/photomaton.m
index b8801d6ef1ebbff96a0f7890dbb4e45a308b9432..34474be94a4e862e00794ed437a7b046579aada2 100644
--- a/wk/photomaton.m
+++ b/wk/photomaton.m
@@ -1,21 +1,21 @@
 if 0
 %% Photomaton : real space
 % FIELD = ni00;   FNAME = 'ni00'; FIELDLTX = '$n_i^{00}$'; XX = RR; YY = ZZ;
-% FIELD = ne00;   FNAME = 'ne00'; FIELDLTX = '$n_e^{00}$'; XX = RR; YY = ZZ;
+FIELD = ne00;   FNAME = 'ne00'; FIELDLTX = '$n_e^{00}$'; XX = RR; YY = ZZ;
 % FIELD = dens_i; FNAME = 'ni';   FIELDLTX = '$n_i$'; XX = RR; YY = ZZ;
 % FIELD = dens_e; FNAME = 'ne';   FIELDLTX = '$n_e$'; XX = RR; YY = ZZ;
 % FIELD = temp_i; FNAME = 'Ti';   FIELDLTX = '$T_i$'; XX = RR; YY = ZZ;
 % FIELD = temp_e; FNAME = 'Te';   FIELDLTX = '$T_e$'; XX = RR; YY = ZZ;
-FIELD = phi; FNAME = 'phi'; FIELDLTX = '$\phi$'; XX = RR; YY = ZZ;
+% FIELD = phi; FNAME = 'phi'; FIELDLTX = '$\phi$'; XX = RR; YY = ZZ;
 % FIELD = drphi; FNAME = 'ZF'; FIELDLTX = '$u^{ZF}_z$'; XX = RR; YY = ZZ;
 % FIELD = -dr2phi; FNAME = 'shear'; FIELDLTX = '$s$'; XX = RR; YY = ZZ;
+ plt = @(x) x./max(max(x));
 TNAME = [];
-tf = 2500; [~,it1] = min(abs(Ts2D-tf)); TNAME = [TNAME,'_',num2str(tf)];
-tf = 3000; [~,it2] = min(abs(Ts2D-tf)); TNAME = [TNAME,'_',num2str(tf)];
-tf = 3500; [~,it3] = min(abs(Ts2D-tf)); TNAME = [TNAME,'_',num2str(tf)];
-tf = 4000; [~,it4] = min(abs(Ts2D-tf)); TNAME = [TNAME,'_',num2str(tf)];
+tf = 500; [~,it1] = min(abs(Ts2D-tf)); TNAME = [TNAME,'_',num2str(tf)];
+tf = 1000; [~,it2] = min(abs(Ts2D-tf)); TNAME = [TNAME,'_',num2str(tf)];
+tf = 1250; [~,it3] = min(abs(Ts2D-tf)); TNAME = [TNAME,'_',num2str(tf)];
+tf = 1750; [~,it4] = min(abs(Ts2D-tf)); TNAME = [TNAME,'_',num2str(tf)];
 fig = figure; FIGNAME = [FNAME,TNAME,'_snaps','_',PARAMS]; set(gcf, 'Position',  [100, 100, 1500, 350])
-plt = @(x) x./max(max(x));
     subplot(141)
         DATA = plt(FIELD(:,:,it1));
         pclr = pcolor((XX),(YY),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
@@ -65,26 +65,44 @@ end
 
 %%
 if 0
-%% Show frame in kspace
-tf = 800; [~,it2] = min(abs(Ts2D-tf)); [~,it5] = min(abs(Ts5D-tf));
-fig = figure; FIGNAME = ['krkz_',sprintf('t=%.0f',Ts2D(it2)),'_',PARAMS];set(gcf, 'Position',  [100, 100, 700, 600])
-CLIM = [0,1];
-    subplot(223); plt = @(x) fftshift((abs(x)),2)/max(max(abs(x)));
-        pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(PHI(:,:,it2))); set(pclr, 'edgecolor','none');
-        caxis(CLIM); colormap hot
-        xlabel('$k_r$'); ylabel('$k_z$'); title(sprintf('$t c_s/R=%.0f$',Ts2D(it2))); legend('$|\hat\phi|$');
-    subplot(222);
-        pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(Ni00(:,:,it2))); set(pclr, 'edgecolor','none');
-        caxis(CLIM); colormap hot
-        xlabel('$k_r$'); ylabel('$k_z$'); legend('$|\hat n_i^{00}|$');
-    subplot(221);
-        pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(Ne00(:,:,it2))); set(pclr, 'edgecolor','none');
-        caxis(CLIM); colormap hot
-        xlabel('$k_r$'); ylabel('$k_z$'); legend('$|\hat n_e^{00}|$');
-    subplot(224);
-    colorbar;
-        caxis(CLIM); colormap hot
-%         pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(Si00(:,:,it5))); set(pclr, 'edgecolor','none'); colorbar;
-%         xlabel('$k_r$'); ylabel('$k_z$');legend('$\hat S_i^{00}$');
+%% Photomaton : k space
+% FIELD = Ni00;   FNAME = 'Ni00'; FIELDLTX = '$N_i^{00}$'; XX = KR; YY = KZ;
+FIELD = Ne00;   FNAME = 'Ne00'; FIELDLTX = '$N_e^{00}$'; XX = KR; YY = KZ;
+FIELD = PHI;   FNAME = 'PHI'; FIELDLTX = '$\tilde\phi$'; XX = KR; YY = KZ;
+FIELD = ifftshift((abs(FIELD)),2); XX = fftshift(XX,2); YY = fftshift(YY,2);
+plt = @(x) x./max(max(x));
+TNAME = [];
+tf = 500; [~,it1] = min(abs(Ts2D-tf)); TNAME = [TNAME,'_',num2str(tf)];
+tf = 1000; [~,it2] = min(abs(Ts2D-tf)); TNAME = [TNAME,'_',num2str(tf)];
+tf = 1250; [~,it3] = min(abs(Ts2D-tf)); TNAME = [TNAME,'_',num2str(tf)];
+tf = 1750; [~,it4] = min(abs(Ts2D-tf)); TNAME = [TNAME,'_',num2str(tf)];
+fig = figure; FIGNAME = [FNAME,TNAME,'_snaps','_',PARAMS]; set(gcf, 'Position',  [100, 100, 1500, 350])
+    subplot(141)
+        DATA = plt(FIELD(:,:,it1));
+        pclr = pcolor((XX),(YY),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
+        colormap gray;
+        xlabel('$r/\rho_s$'); ylabel('$z/\rho_s$');set(gca,'ytick',[]); 
+        title(sprintf('$t c_s/R=%.0f$',Ts2D(it1)));
+    subplot(142)
+        DATA = plt(FIELD(:,:,it2));
+        pclr = pcolor((XX),(YY),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
+        colormap(bluewhitered);
+        xlabel('$r/\rho_s$');ylabel('$z/\rho_s$'); set(gca,'ytick',[]); 
+        title(sprintf('$t c_s/R=%.0f$',Ts2D(it2)));
+    subplot(143)
+        DATA = plt(FIELD(:,:,it3));
+        pclr = pcolor((XX),(YY),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
+        colormap(bluewhitered);
+        xlabel('$r/\rho_s$');ylabel('$z/\rho_s$');set(gca,'ytick',[]); 
+        title(sprintf('$t c_s/R=%.0f$',Ts2D(it3)));
+    subplot(144)
+        DATA = plt(FIELD(:,:,it4));
+        pclr = pcolor((XX),(YY),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
+        colormap(bluewhitered);
+        xlabel('$r/\rho_s$');ylabel('$z/\rho_s$'); set(gca,'ytick',[]); 
+        title(sprintf('$t c_s/R=%.0f$',Ts2D(it4)));
+    suptitle([FIELDLTX,', $\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB),...
+        ', $L=',num2str(L),'$, $N=',num2str(Nr),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',...
+        ' $\mu_{hd}=$',num2str(MU)]);
 save_figure
 end
\ No newline at end of file