diff --git a/wk/lin_EPY.m b/wk/lin_EPY.m
new file mode 100644
index 0000000000000000000000000000000000000000..41abd8aa00081b725df41e4b2ed803da719774ea
--- /dev/null
+++ b/wk/lin_EPY.m
@@ -0,0 +1,179 @@
+%% QUICK RUN SCRIPT for linear entropy mode in a Zpinch
+% This script create a directory in /results and run a simulation directly
+% from matlab framework. It is meant to run only small problems in linear
+% for benchmark and debugging purpose since it makes matlab "busy"
+%
+SIMID = 'lin_EPY'; % Name of the simulation
+RUN = 1; % To run or just to load
+addpath(genpath('../matlab')) % ... add
+default_plots_options
+PROGDIR = '/home/ahoffman/gyacomo/';
+% EXECNAME = 'gyacomo_1.0';
+EXECNAME = 'gyacomo';
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Set Up parameters
+CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% PHYSICAL PARAMETERS
+NU = 0.1; % Collision frequency
+TAU = 1.0; % e/i temperature ratio
+K_Ne = 2.2; % ele Density '''
+K_Te = K_Ne/4; % ele Temperature '''
+K_Ni = K_Ne; % ion Density gradient drive
+K_Ti = K_Ni/4; % ion Temperature '''
+SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
+KIN_E = 1; % 1: kinetic electrons, 2: adiabatic electrons
+BETA = 0.0; % electron plasma beta
+%% GRID PARAMETERS
+P = 4;
+J = P/2;
+PMAXE = P; % Hermite basis size of electrons
+JMAXE = J; % Laguerre "
+PMAXI = P; % " ions
+JMAXI = J; % "
+NX = 2; % real space x-gridpoints
+NY = 100; % '' y-gridpoints
+LX = 2*pi/0.8; % Size of the squared frequency domain
+LY = 120;%2*pi/0.05; % Size of the squared frequency domain
+NZ = 1; % number of perpendicular planes (parallel grid)
+NPOL = 1;
+SG = 0; % Staggered z grids option
+%% GEOMETRY
+GEOMETRY= 'Z-pinch'; % Z-pinch overwrites q0, shear and eps
+Q0 = 1.4; % safety factor
+SHEAR = 0.8; % magnetic shear
+NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
+EPS = 0.18; % inverse aspect ratio
+%% TIME PARMETERS
+TMAX = 50; % Maximal time unit
+DT = 1e-2; % Time step
+SPS0D = 1; % Sampling per time unit for 2D arrays
+SPS2D = 0; % Sampling per time unit for 2D arrays
+SPS3D = 2; % Sampling per time unit for 2D arrays
+SPS5D = 1/5; % Sampling per time unit for 5D arrays
+SPSCP = 0; % Sampling per time unit for checkpoints
+JOB2LOAD= -1;
+%% OPTIONS
+LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
+% Collision operator
+% (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
+CO = 'DG';
+GKCO = 1; % gyrokinetic operator
+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= 'mom00'; % Start simulation with a noisy mom00/phi/allmom
+%% OUTPUTS
+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;
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% unused
+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;
+%%-------------------------------------------------------------------------
+%% RUN
+setup
+% system(['rm fort*.90']);
+% Run linear simulation
+if RUN
+ system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',PROGDIR,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk'])
+end
+
+%% Load results
+%%
+filename = [SIMID,'/',PARAMS,'/'];
+LOCALDIR = [PROGDIR,'results/',filename,'/'];
+% Load outputs from jobnummin up to jobnummax
+JOBNUMMIN = 00; JOBNUMMAX = 00;
+data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
+
+%% Short analysis
+if 1
+%% linear growth rate (adapted for 2D zpinch and fluxtube)
+trange = [0.5 1]*data.Ts3D(end);
+nplots = 1; % 1 for only growth rate and error, 2 for omega local evolution, 3 for plot according to z
+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);
+end
+
+if 0
+%% Ballooning plot
+options.time_2_plot = [120];
+options.kymodes = [0.9];
+options.normalized = 1;
+% options.field = 'phi';
+fig = plot_ballooning(data,options);
+end
+
+if 0
+%% Hermite-Laguerre spectrum
+% options.TIME = 'avg';
+options.P2J = 1;
+options.ST = 1;
+options.PLOT_TYPE = 'space-time';
+% options.PLOT_TYPE = 'Tavg-1D';
+% options.PLOT_TYPE = 'Tavg-2D';
+options.NORMALIZED = 0;
+options.JOBNUM = 0;
+options.TIME = [0 50];
+options.specie = 'i';
+options.compz = 'avg';
+fig = show_moments_spectrum(data,options);
+% fig = show_napjz(data,options);
+save_figure(data,fig)
+end
+
+if 0
+%% linear growth rate for 3D Zpinch (kz fourier transform)
+trange = [0.5 1]*data.Ts3D(end);
+options.keq0 = 1; % chose to plot planes at k=0 or max
+options.kxky = 1;
+options.kzkx = 0;
+options.kzky = 0;
+[lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options);
+save_figure(data,fig)
+end
+if 0
+%% Mode evolution
+options.NORMALIZED = 0;
+options.K2PLOT = 1;
+options.TIME = [0:1000];
+options.NMA = 1;
+options.NMODES = 1;
+options.iz = 'avg';
+fig = mode_growth_meter(data,options);
+save_figure(data,fig,'.png')
+end
+
+
+if 0
+%% RH TEST
+ikx = 2; iky = 2; t0 = 0; t1 = data.Ts3D(end);
+[~, it0] = min(abs(t0-data.Ts3D));[~, it1] = min(abs(t1-data.Ts3D));
+plt = @(x) squeeze(mean(real(x(iky,ikx,:,it0:it1)),3))./squeeze(mean(real(x(iky,ikx,:,it0)),3));
+figure
+plot(data.Ts3D(it0:it1), plt(data.PHI));
+xlabel('$t$'); ylabel('$\phi_z(t)/\phi_z(0)$')
+title(sprintf('$k_x=$%2.2f, $k_y=$%2.2f',data.kx(ikx),data.ky(iky)))
+end