diff --git a/matlab/read_flux_out_XX.m b/matlab/read_flux_out_XX.m
index 22f4ce581ac07712ec5492602e0ff5f354c28395..e442f5235b852bbaa0cece0e5c21ee66b968dd04 100644
--- a/matlab/read_flux_out_XX.m
+++ b/matlab/read_flux_out_XX.m
@@ -1,4 +1,4 @@
-function [out] = read_flux_out_XX(folderPath,PLOT)
+function [out] = read_flux_out_XX(folderPath,PLOT,nmvm)
     % Check if the prompt_string is provided as an argument
     if nargin < 1
         % If not provided, prompt the user for input
@@ -6,8 +6,12 @@ function [out] = read_flux_out_XX(folderPath,PLOT)
         % Get the input from the user
         folderPath = input(prompt_string, 's');   
         PLOT = 1;
-    else if nargin == 1
+        nmvm = 1;
+    elseif nargin == 1
             PLOT = 0;
+            nmvm = 1;
+    elseif nargin == 2
+            nmvm = 1;
     end
 
     % Initialize empty arrays to store the values
@@ -74,17 +78,23 @@ function [out] = read_flux_out_XX(folderPath,PLOT)
     if PLOT
     figure
     subplot(211)
-    plot(t_all,Pxi_all,'r','DisplayName','ions'); hold on;
+    x_ = movmean(t_all,nmvm);
+    y_ = movmean(Pxi_all,nmvm);
+    plot(x_,y_,'r','DisplayName','ions'); hold on;
     if(numel(t_all)==numel(Pxe_all))
-        plot(t_all,Pxe_all,'-.b','DisplayName','electrons'); hold on;
+        y_ = movmean(Pxe_all,nmvm);
+        plot(x_,y_,'-.b','DisplayName','electrons'); hold on;
     end
     xlabel('$tc_s/R$'); ylabel('$\Gamma_{x}$');
     legend('show')
     title('Radial particle flux')
     subplot(212)
-    plot(t_all,Qxi_all,'r','DisplayName','ions'); hold on;
+    x_ = movmean(t_all,nmvm);
+    y_ = movmean(Qxi_all,nmvm);
+    plot(x_,y_,'r','DisplayName','ions'); hold on;
     if(numel(t_all)==numel(Qxe_all))
-        plot(t_all,Qxe_all,'-.b','DisplayName','electrons'); hold on;
+        y_ = movmean(Qxe_all,nmvm);
+        plot(x_,y_,'-.b','DisplayName','electrons'); hold on;
     end
     xlabel('$tc_s/R$'); ylabel('$Q_{x}$');
     legend('show');
diff --git a/wk/ZF_energy_t_evolution.m b/wk/ZF_energy_t_evolution.m
new file mode 100644
index 0000000000000000000000000000000000000000..6075cd1aa1e70570ea685243c873920ceee7fb00
--- /dev/null
+++ b/wk/ZF_energy_t_evolution.m
@@ -0,0 +1,10 @@
+sz_ = size(data.PHI);
+iz  = sz_(3)/2+1;
+Ezf = sum(abs(data.PHI(:,:,:,:)),2);
+Ezf = squeeze(Ezf(1,:,13,:));
+Enz = sum(abs(data.PHI(:,:,:,:)),2);
+Enz = squeeze(sum(Enz(2:end,:,13,:),1));
+
+figure; 
+plot(data.Ts3D,Ezf,'DisplayName','Zonal energy'); hold on;
+plot(data.Ts3D,Enz,'DisplayName','Non-zonal energy');
diff --git a/wk/fast_analysis.m b/wk/fast_analysis.m
index 4e0e5814654e347f9d832e3af220307ef7567f15..4bc0c79a1248a38a6d7eda8b97476f2de220fbd1 100644
--- a/wk/fast_analysis.m
+++ b/wk/fast_analysis.m
@@ -6,97 +6,29 @@ addpath(genpath([gyacomodir,'matlab/load'])) % ... add
 default_plots_options
 % Partition of the computer where the data have to be searched
 % PARTITION='/Users/ahoffmann/gyacomo/results/paper_3/';
-PARTITION='/misc/gyacomo23_outputs/paper_3/';
-% PARTITION = '../results/paper_3/';
-%% Paper 3
-% resdir = 'DTT_rho85/3x2x192x48x32';
-% resdir = 'DTT_rho85/3x2x192x48x32_NT';
-% resdir = 'DTT_rho98/3x2x192x48x32';
-% resdir = 'DTT_rho98/3x2x192x48x32_0.25grad';
-% resdir = 'LM_DIIID_rho95/5x3x512x92x32';
-% resdir = 'LM_DIIID_rho95/3x2x512x92x32';
-% resdir = 'DIIID_LM_rho90/3x2x256x128x32';
-% resdir = 'DTT_rho85_geom_scan/P8_J4_delta_nuDGGK_conv_test/delta_-0.3_nu_0.9';
-% resdir = 'NT_DIIID_Austin2019_rho95/3x2x256x64x32';
-
-% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_2500/largerbox_moremodes';
-% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_2500/PT/3x2x128x32x32';
-% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_2500/NT/3x2x128x32x32';
-% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_2500/VNT/3x2x128x32x32';
-% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_2500/VPT/3x2x128x32x32';
-% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_2500/CIRC/3x2x128x32x32';
-% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_2500/ELONG/3x2x128x32x32';
-% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_2500/PT/lin_3x2x128x32x32';
-% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_2500/NT/lin_3x2x128x32x32';
-
-% resdir = 'DIIID_rho95_cold_ions_tau1e-3/huge';
-
-% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_1000/PT/3x2x128x32x32';
-% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_1000/NT/3x2x128x32x32';
-% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_1000/PT/lin_3x2x128x32x32';
-% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_1000/NT/lin_3x2x128x32x32';
-% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_750/PT/lin_3x2x128x32x32';
-% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_750/NT/lin_3x2x128x32x32';
-% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_250/PT/lin_3x2x128x32x32';
-% resdir = 'DIIID_rho95_cold_ions_tau1e-3/RT_250/NT/lin_3x2x128x32x32';
-% resdir = 'DIIID_rho95_cold_ions_tau1e-3/huge';
-
-% resdir = 'DIIID_rho95_cold_ions_tau1e-2/PT/3x2x128x32x32';
-% resdir = 'DIIID_rho95_cold_ions_tau1e-2/NT/3x2x128x32x32';
-% resdir = 'DIIID_rho95_cold_ions_tau1e0/PT/3x2x128x32x32';
-% resdir = 'DIIID_rho95_cold_ions_tau1e0/PT/5x3x192x48x32';
-% resdir = 'DIIID_rho95_cold_ions_tau1e0/NT/3x2x128x32x32';
-% resdir = 'DIIID_rho95_cold_ions_tau1e0/NT/5x3x192x48x32';
-% resdir = 'DTT_rho85_cold_ions_tau1e-3/NT/3x2x128x32x32';
-% % resdir = '../testcases/cyclone_example';
-% resdir = '../testcases/CBC_ExBshear/std';
-% resdir = '../results/paper_3/HM_DTT_rho98/3x2x128x64x64';
-
-% resdir = 'DIIID_cold_ions_rho95_geom_scan/3x2x192x48x32_RT_1000_eps_q0_scan/NT/eps_0.35_q0_4.0';
-% resdir = 'DTT_rho85_geom_scan/P2_J1_PT_sfact_shear_scan/shear_2.7_q0_-2.9';
-% PARTITION = ''; resdir ='../results/HEL_CBC/tau_1e-3_kT_3500/128x32x24';
-% PARTITION = ''; resdir ='../results/HEL_CBC/tau_1e-3_kT_3500/256x64x24';
-% PARTITION = ''; resdir ='../results/HEL_CBC/tau_1e-3_colless';
-% PARTITION = ''; resdir ='../results/HEL_CBC/tau_1e-3_kT_2000/lin_128x32x24';
-% PARTITION = ''; resdir ='../results/HEL_CBC/tau_1e-3_kT_2000/128x32x24';
-% PARTITION = ''; resdir ='../results/HEL_CBC/CBC_21/128x32x24';
-% PARTITION = ''; resdir ='../results/HEL_CBC/192x48x24';
-% PARTITION = ''; resdir ='../testcases/Ivanov_2020';
-% PARTITION = ''; resdir ='../testcases/Hasegawa_Wakatani';
-% resdir = 'HEL_CBC/256x92x24_max_trunc';
-
-% resdir ='HEL_CBC/192x48x24';
-% resdir ='HEL_CBC/256x92x24'3;
-% resdir ='HEL_CBC/256x256x32/k_N__0.0_k_T__1750';
-
-PARTITION = '/misc/gyacomo23_outputs/paper_3/DIIID_rho_95_Tstudy/';
-% resdir = 'multi_scale_3x2x512x128x24';
-% resdir = 'multi_scale_3x2x512x128x24_larger_box';
-% resdir = 'multi_scale_3x2x768x192x24/continue_with_gradN_and_tau';
-% resdir = 'multi_scale/multi_scale_5x2x768x192x24/PT';
-% resdir = 'NT';
-% resdir = 'electron_scale_3x2x256x64x24';
-% resdir = 'electron_scale_3x2x128x64x24';
-% resdir = 'ion_scale_3x2x192x48x32_larger_box';
-% resdir = 'ion_scale_3x2x256x64x32'; % PT and NT also
-% resdir = 'ion_scale/ion_scale_5x2x256x64x32_tau_1_RN_0/NT';
-% resdir = 'adiab_e/5x2x256x64x32_tau_1_RN_0/NT';
-
-resdir = 'multi_scale/multi_scale_3x2x768x192x24/NT';
-% resdir = 'multi_scale/multi_scale_3x2x512x128x24_larger_box/NT';
-% resdir = 'multi_scale/multi_scale_3x2x768x192x24/continue_with_gradN_and_tau';
-% resdir = 'ion_scale/ion_scale_5x2x256x64x32_tau_1_RN_0/NT';
-% resdir = 'adiab_e/5x2x256x64x32_tau_1_RN_0/0T';
-% resdir = 'hot_electrons/hot_electrons_256x64x24/0T';
-
+% PARTITION='/misc/gyacomo23_outputs/paper_3/';
+% resdir = 'DIIID_rho_95_Tstudy/ion_scale/ion_scale_5x2x256x64x32_tau_1_RN_0/0T/';
+% resdir = 'DIIID_rho_95_Tstudy/adiab_e/5x2x256x64x32_tau_1_RN_0/0T/';
+%%
 
+PARTITION = '/misc/gyacomo23_outputs/triangularity_paper/';
 
+% resdir = 'ion_scale/3x2x256x64x32/PT';
+% resdir = 'ion_scale/5x2x256x64x32/0T';
+% resdir = 'adiabatic_electrons/5x2x256x64x32/0T';
+% resdir = 'adiabatic_electrons/5x2x192x48x24/NT';
+% resdir = 'hot_electrons/256x64x32/0T';
+% resdir = 'hot_electrons/L_300/256x64x32/0T';
+resdir = 'hot_electrons/L_300_gradN_scaled/256x64x32/0T';
+% resdir = 'hot_electrons/256x64x32/0T';
+% resdir = 'hot_electrons/512x64x32/0T';
 % PARTITION = '/misc/gyacomo23_outputs/paper_3/';
 % resdir = 'DIIID_HEL_rho95/PT';
 
 DATADIR = [PARTITION,resdir,'/'];
+% read_flux_out_XX(DATADIR,1,5)
 %%
-J0 = 00; J1 = 20;
+J0 = 00; J1 = 10;
 
 % Load basic info (grids and time traces)
 data    = {};
@@ -132,11 +64,11 @@ options.TAVG_1   = data.Ts3D(end);
 options.NCUT     = 5;              % 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 = 'n_e';          % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
+% options.ST_FIELD = 'n_i';          % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
 % options.ST_FIELD = 'u_i';          % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
-% options.ST_FIELD = 'T_i';          % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
+options.ST_FIELD = 'T_i';          % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
 % options.ST_FIELD = 'n_i T_i';          % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
-options.ST_FIELD = '\phi';          % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
+% options.ST_FIELD = '\phi';          % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
 % options.ST_FIELD = 'Q_{xi}';          % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
 % options.ST_FIELD = 'G_x';          % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
 % options.ST_FIELD = 'w_{Ez}';          % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
@@ -157,13 +89,13 @@ options.NORMALIZE = 0;
 options.LOGSCALE  = 0;
 options.CLIMAUTO  = 1;
 options.TAVG      = 1;
-% options.NAME      = ['N_e^{00}'];
+% options.NAME      = ['N_i^{00}'];
 % options.NAME      = 'n_e';
-% % options.NAME      = 'u_i';
-% options.NAME      = 'T_i';
-% options.NAME      = 'Q_{xe}';
-options.NAME      = 'v_{Ey}';
-% options.NAME      = 'w_{Ez}';
+% options.NAME      = 'u_i';
+% options.NAME      = 'n_i';
+% options.NAME      = 'Q_{xi}';
+% options.NAME      = 'v_{Ey}';
+options.NAME      = 'w_{Ez}';
 % options.NAME      = '\omega_z';
 % options.NAME      = '\phi';
 % options.NAME      = 'n_i-n_e';
@@ -172,10 +104,11 @@ loc =11;
 options.COMP =i_;
 % options.PLAN      = '3D';  
 options.PLAN      = 'xy'; options.COMP =floor(data.grids.Nz/2)+1; 
-% options.PLAN      = 'yz'; options.COMP ='avg';
+% options.PLAN      = 'xz'; options.COMP ='avg';
 % options.COMP ='avg'; 
 options.XYZ  =[-11 20 0]; 
-options.TIME      = [90:150];
+options.TIME = [100 250 350 500]; options.TAVG = 0;
+% options.TIME = [100:250]; options.TAVG = 1;
 options.RESOLUTION = 256;
 fig = photomaton(data,options);
 % colormap(gray)
@@ -209,40 +142,42 @@ options.FIELD  = 'Ni00';
 options.GOK2   = 0;
 options.SHOWFIG = 1;
 [fig, wkykx, ekykx] = mode_growth_meter(data,options);
-%%
-kx = (1:data.grids.Nx/2)'*2*pi/data.fort_00.GRID.Lx;
-ky = (1:data.grids.Ny/2)'*2*pi/data.fort_00.GRID.Ly;
-gkxky = real(wkykx(2:end,1:data.grids.Nx/2))';
-gkxky(isnan(gkxky)) =0;
-gkxky(isinf(gkxky)) =0;
-% gkxky(gkxky<0)      =0;
-% gkxky = imgaussfilt(gkxky,1);
-%
-wkxky = imag(wkykx(2:end,1:data.grids.Nx/2))';
-wkxky(isnan(wkxky)) =0;
-wkxky(isinf(wkxky)) =0;
-% wkxky(wkxky<0)      =0;
-% wkxky = imgaussfilt(wkxky,1.5);
-%
-figure; 
-subplot(121)
-    contourf(kx,ky,gkxky',10)
-    % clim(0.5*[0 1]); 
-    % colormap(bluewhitered); colorbar;
-    xlim([0.025 1]);
-    xlabel('$k_x\rho_s$'); ylabel('$k_y\rho_s$')
-subplot(122)
-    contourf(kx,ky,wkxky',10)
-    % clim(1*[0 1]); 
-    % colormap(bluewhitered); colorbar 
-    xlim([0.025 1]);
-    xlabel('$k_x\rho_s$'); ylabel('$k_y\rho_s$')
-% save_figure(data,fig,'.png')
+% %%
+% kx = (1:data.grids.Nx/2)'*2*pi/data.fort_00.GRID.Lx;
+% ky = (1:data.grids.Ny/2)'*2*pi/data.fort_00.GRID.Ly;
+% gkxky = real(wkykx(2:end,1:data.grids.Nx/2))';
+% gkxky(isnan(gkxky)) =0;
+% gkxky(isinf(gkxky)) =0;
+% % gkxky(gkxky<0)      =0;
+% % gkxky = imgaussfilt(gkxky,1);
+% %
+% wkxky = imag(wkykx(2:end,1:data.grids.Nx/2))';
+% wkxky(isnan(wkxky)) =0;
+% wkxky(isinf(wkxky)) =0;
+% % wkxky(wkxky<0)      =0;
+% % wkxky = imgaussfilt(wkxky,1.5);
+% %
+% figure; 
+% subplot(121)
+%     contourf(kx,ky,gkxky',10)
+%     % clim(0.5*[0 1]); 
+%     % colormap(bluewhitered); colorbar;
+%     xlim([0.025 1]);
+%     xlabel('$k_x\rho_s$'); ylabel('$k_y\rho_s$')
+% subplot(122)
+%     contourf(kx,ky,wkxky',10)
+%     % clim(1*[0 1]); 
+%     % colormap(bluewhitered); colorbar 
+%     xlim([0.025 1]);
+%     xlabel('$k_x\rho_s$'); ylabel('$k_y\rho_s$')
+% % save_figure(data,fig,'.png')
 end
 
 if 1
 %% Hermite-Laguerre spectrum
-[data.Napjz, data.Ts3D] = compile_results_3Da(DATADIR,0,10,'Napjz');
+[data.Napjz, data.Ts3D] = compile_results_3Da(DATADIR,J0,J1,'Napjz');
+data.Napjz(1,3,1,:,:) = data.Napjz(1,3,1,:,:)*data.inputs.tau;
+data.Napjz(1,1,2,:,:) = data.Napjz(1,1,2,:,:)*data.inputs.tau;
 % [data.Napjz, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'Nipjz');
 options.ST         = 1;
 options.NORMALIZED = 0;
@@ -271,13 +206,14 @@ end
 
 if 1
 %% 1D spectral plot
-options.TIME  = [30 80]; % averaging time window
+options.TIME  = [100 300]; % averaging time window
 % options.NAME      = ['N_i^{00}'];
 % options.NAME      = 'n_i';
 % options.NAME      = 'T_i';
 % options.NAME      = 'Q_{xi}';
 % options.NAME      = 's_{Ey}';
-options.NAME      = '\psi';
+options.NAME      = '\phi';
+% options.NAME      = '\psi';
 options.NORMALIZE = 0;
 [fig] = plot_spectrum(data,options);
 end
@@ -302,15 +238,14 @@ options.CLIMAUTO  = 0; % adjust the colormap auto
 % options.NAME      = '\phi';
 % options.NAME      = 'w_{Ez}';
 % options.NAME      = '\psi';
-options.NAME      = 'T_i';
+% options.NAME      = 'n_i';
 % options.NAME      = '\phi^{NZ}';
-% options.NAME     = ['N_e^{00}'];
+options.NAME     = ['N_i^{00}'];
 % options.NAME     = ['N_i^{00}'];
-options.PLAN      = 'xy';
+options.PLAN      = 'xy'; options.COMP =floor(data.grids.Nz/2)+1; 
+% options.PLAN      = 'xz'; options.COMP ='avg';
 % options.PLAN      = '3D';  
-% options.XYZ  =[-11 20 0]; 
-% options.COMP      = 'avg';
-options.COMP      = floor(data.grids.Nz/2+1);
+options.XYZ  =[-11 20 0]; 
 options.TIME      =  data.Ts3D(1:1:end);
 % options.TIME      = [0:1500];
 data.EPS          = 0.1;
diff --git a/wk/lin_run_script.m b/wk/lin_run_script.m
index 7e0eb7271fab42b84a51e3fea2c12c986b5eee4b..09ace133bdbbe5aeaa9e1d5b6e4ed030d48aeb0f 100644
--- a/wk/lin_run_script.m
+++ b/wk/lin_run_script.m
@@ -26,7 +26,7 @@ EXECNAME = 'gyacomo23_dp'; % double precision
 % run lin_DTT_HM_rho85
 % run lin_DTT_HM_rho98
 % run lin_DIIID_LM_rho90
-% run lin_DIIID_LM_rho95
+run lin_DIIID_LM_rho95
 % run lin_DIIID_LM_rho95_HEL
 % run lin_JET_rho97
 % run lin_Entropy
@@ -34,30 +34,31 @@ EXECNAME = 'gyacomo23_dp'; % double precision
 % run lin_KBM
 % run lin_RHT
 % run lin_Ivanov
-rho  = 0.95; TRIANG = 'NT'; READPROF = 1; 
+% rho  = 0.95; TRIANG = 'NT'; READPROF = 1; 
 % prof_folder = ['parameters/profiles/DIIID_Austin_et_al_2019/',TRIANG,'/'];
 % prof_folder = ['parameters/profiles/DIIID_Oak_Nelson/',TRIANG,'/'];
-prof_folder = ['parameters/profiles/DIIID_Oak_Nelson_high_density/',TRIANG,'/'];
-run lin_DIIID_data
+% prof_folder = ['parameters/profiles/DIIID_Oak_Nelson_high_density/',TRIANG,'/'];
+% run lin_DIIID_data
 % run lin_STEP_EC_HD_psi71
 % run lin_STEP_EC_HD_psi49
-if 1
+if 0
 % Plot the profiles
  plot_params_vs_rho(geom,prof_folder,rho,0.5,1.1,Lref,mref,Bref,READPROF);
 end
 % SIMID = ['rho_scan_DIIID_AUSTIN_2019/3x2x192x96x32/rho',num2str(rho)];
 %% Change parameters
 % GEOMETRY = 's-alpha';
+PMAX  = 2; JMAX = 1;
 DELTA =0.0; 
-K_Ni = 0; K_Ne = 0;
+% K_Ni = 0; K_Ne = 0;
 % DELTA = 0.0; 
 % DELTA = 0.2; 
 S_DELTA = DELTA/2;
-LY   = 2*pi/0.25;
-TMAX = 20;
+LY   = 2*pi/0.05;
+TMAX = 40;
 NY   = 2;
 DT   = 0.01;
-TAU  = 1; NU = 0.05;
+% TAU  = 1; NU = 0.05;
 % TAU = 1e-3; K_Ti = K_Ti/2/TAU; NU = 3*NU/8/TAU; ADIAB_E = 1; NA = 1;
 % MU_X = 1; MU_Y = 1;
 %% RUN
diff --git a/wk/lin_scan_script.m b/wk/lin_scan_script.m
index 3883f0e9b4fff3be942ca2dc74e8f14c0b89dbbc..f95c86bb848fcafe6a53f647c7f219270d67e5e5 100644
--- a/wk/lin_scan_script.m
+++ b/wk/lin_scan_script.m
@@ -28,28 +28,31 @@ EXECNAME = 'gyacomo23_dp'; % double precision
 % run lin_Entropy
 % run lin_ITG
 % run lin_RHT
-rho  = 0.95; TRIANG = 'PT'; READPROF = 1; 
+% rho  = 0.95; TRIANG = 'PT'; READPROF = 1; 
 % prof_folder = ['parameters/profiles/DIIID_Austin_et_al_2019/',TRIANG,'/'];
 % prof_folder = ['parameters/profiles/DIIID_Oak_Nelson/',TRIANG,'/'];
-prof_folder = ['parameters/profiles/DIIID_Oak_Nelson_high_density/',TRIANG,'/'];
-run lin_DIIID_data
+% prof_folder = ['parameters/profiles/DIIID_Oak_Nelson_high_density/',TRIANG,'/'];
+% run lin_DIIID_data
+run lin_DIIID_LM_rho95
 
 %% Change parameters
 % NU   = 1;
 % TAU  = 1;
 NY   = 2;
-EXBRATE = 0;
-S_DELTA = min(2.0,S_DELTA);
-SIGMA_E  = 0.023;
-NEXC = 0;
+% EXBRATE = 0;
+% S_DELTA = min(2.0,S_DELTA);
+% SIGMA_E  = 0.023;
+% NEXC = 0;
 LX   = 120;
 %% Scan parameters
 SIMID = [SIMID,'_scan'];
-P_a   = [2 4];
+P_a   = [2 4]; J_a = [1 1];
 % P_a   = 2;
-ky_a  = [0.01 0.02 0.05 0.1  0.2  0.5  1.0  2.0  5.0  10.0];
+% ky_a  = [0.01 0.02 0.05 0.1  0.2  0.5  1.0  2.0  5.0  10.0];
+ky_a  = [0.025:0.025:1.1];
 % ky_a  = 4.0;
-dt_a  = logspace(-2,-3,numel(ky_a));
+% dt_a  = logspace(-2,-3,numel(ky_a));
+dt_a  = linspace(0.01,0.01,numel(ky_a));
 CO    = 'DG';
 %% Scan loop
 % arrays for the result
@@ -59,7 +62,7 @@ w_ky = g_ky*0;
 w_std= g_ky*0;
 j = 1;
 for PMAX = P_a
-    JMAX = P/2;
+    JMAX = J_a(j);
     i = 1;
     for ky = ky_a
         LY   = 2*pi/ky;
diff --git a/wk/parameters/lin_DIIID_LM_rho95.m b/wk/parameters/lin_DIIID_LM_rho95.m
index 17c6cbf8552518eb05c6cbd2b62d9671e59d0730..3dfccf17b299e2bf4e454396b26a87b0f210dfc0 100644
--- a/wk/parameters/lin_DIIID_LM_rho95.m
+++ b/wk/parameters/lin_DIIID_LM_rho95.m
@@ -1,25 +1,23 @@
 %% Reference values
 % See Neiser et al. 2019 Gyrokinetic GENE simulations of DIII-D near-edge L-mode plasmas
 %% Set simulation parameters
-SIMID   = 'lin_DIIID_HM_rho90';  % Name of the simulation
+SIMID   = 'lin_DIIID_LM_rho95';  % Name of the simulation
 %% Set up physical parameters
 CLUSTER.TIME = '99:00:00';  % Allocation time hh:mm:ss
-NU = 1.0; %(0.00235 in GENE)
-TAU = 0.281/0.0831;               % i/e temperature ratio
-K_Ne    = 7.05;             % ele Density '''
-K_Te    = 13.5;             % ele Temperature '''
+NU = 0.02; %(2x)
+TAU = 1.6;               % i/e temperature ratio
+K_Ne    = 1.7;             % ele Density '''
+K_Te    = 6.0;             % ele Temperature '''
 K_Ni    = K_Ne;             % ion Density gradient drive
-K_Ti    = 2.32;             % ion Temperature '''
-SIGMA_E = 0.0233380/sqrt(2);        % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
+K_Ti    = 5.2;             % ion Temperature '''
+SIGMA_E = sqrt(1/1000);%0.0233380/sqrt(2) % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
 NA      = 2;          % number of kinetic species
 ADIAB_E = (NA==1);          % adiabatic electron model
-BETA    = 1.32e-2*0.01;           % electron plasma beta in prct
+BETA    = 7.59e-4;           % electron plasma beta in prct
 MHD_PD  = 1;
 %% Set up grid parameters
-P = 2;
-J = P/2;%P/2;
-PMAX = P;                   % Hermite basis size
-JMAX = J;                   % Laguerre basis size
+PMAX = 4;                   % Hermite basis size
+JMAX = 1;                   % Laguerre basis size
 NX = 6;                    % real space x-gridpoints
 NY = 2;                     % real space y-gridpoints
 LX = 2*pi/0.1;              % Size of the squared frequency domain in x direction
@@ -31,15 +29,15 @@ NEXC = 1;                   % To extend Lx if needed (Lx = Nexc/(kymin*shear))
 %% GEOMETRY
 % GEOMETRY= 's-alpha';
 GEOMETRY= 'miller';
-Q0      = 5.18;    % safety factor
-SHEAR   = 4.47;    % magnetic shear
-EPS     = 0.28;    % inverse aspect ratio
-KAPPA   = 1.53;    % elongation
-S_KAPPA = 0.77;
-DELTA   = 0.23;    % triangularity
-S_DELTA = 1.05;
-ZETA    =-0.01;    % squareness
-S_ZETA  =-0.17;
+Q0      = 4.8;    % safety factor
+SHEAR   = 2.5;    % magnetic shear
+EPS     = 0.3;    % inverse aspect ratio
+KAPPA   = 1.6;    % elongation
+S_KAPPA = 0.5;
+DELTA   = 0.0;    % triangularity
+S_DELTA = 0.0;
+ZETA    = 0.0;    % squareness
+S_ZETA  = 0.0;
 PARALLEL_BC = 'dirichlet'; % Boundary condition for parallel direction ('dirichlet','periodic','shearless','disconnected')
 SHIFT_Y = 0.0;    % Shift in the periodic BC in z
 NPOL   = 1;       % Number of poloidal turns
diff --git a/wk/plot_dist_func_and_trap_pass_el.m b/wk/plot_dist_func_and_trap_pass_el.m
new file mode 100644
index 0000000000000000000000000000000000000000..b07ad0bfbb375aed9f17c7e3961cf05fa13514e1
--- /dev/null
+++ b/wk/plot_dist_func_and_trap_pass_el.m
@@ -0,0 +1,19 @@
+options.SHOW_FLUXSURF = 0;
+options.SHOW_METRICS  = 0;
+options.SHOW_CURVOP   = 0;
+[fig, geo_arrays] = plot_metric(data,options);
+
+eps_eff = 1/max((geo_arrays.hatB));
+T = [10];
+s = linspace(-3,3,64);
+x = linspace(0,8,48);
+[SS,XX,FF] = compute_fa_2D(data,'i',s,x,T);
+
+figure
+% contour(SS,XX,log10(FF))
+pca = pcolor(SS,XX,(FF)); set(pca,'EdgeColor','None');
+set(gca,'ColorScale','log'); %shading interp;
+% clim([1e-4 1]);
+hold on
+plot(s,(s.^2/(1+data.fort_00.GEOMETRY.eps)),'--w');
+plot(s,eps_eff*(s.^2),'--k');
\ No newline at end of file