From 1fda930a40c53a5d30079429309058d4f619c59e Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Wed, 11 May 2022 10:01:57 +0200
Subject: [PATCH] scripts update

---
 matlab/compute/compute_fa_1D.m |  10 +--
 matlab/load/load_pjz_data.m    |  16 ++++
 matlab/plot/plot_ballooning.m  |  33 ++++----
 matlab/plot/plot_fa.m          |   6 +-
 matlab/plot/plot_fa_gene.m     |  12 +--
 matlab/plot/show_napjz.m       | 136 +++++++++++++++++++++++++++++++++
 wk/analysis_3D.m               |  42 +++++-----
 wk/gene_analysis_3D.m          |  96 ++---------------------
 wk/quick_run.m                 |  37 ++++-----
 9 files changed, 234 insertions(+), 154 deletions(-)
 create mode 100644 matlab/load/load_pjz_data.m
 create mode 100644 matlab/plot/show_napjz.m

diff --git a/matlab/compute/compute_fa_1D.m b/matlab/compute/compute_fa_1D.m
index d8ab8d74..79fe35d8 100644
--- a/matlab/compute/compute_fa_1D.m
+++ b/matlab/compute/compute_fa_1D.m
@@ -54,7 +54,7 @@ if options.non_adiab
         for ikx = 1:data.Nkx
             for iky = 1:data.Nky    
                 kp_ = sqrt(kx_(ikx)^2 + ky_(iky)^2);
-                Napj_(1,ij_,ikx,iky) = Napj_(1,ij_,ikx,iky) + kernel(ij_,kp_)*phi_(ikx,iky);
+                Napj_(1,ij_,iky,ikx) = Napj_(1,ij_,iky,ikx) + kernel(ij_,kp_)*phi_(iky,ikx);
             end
         end
     end
@@ -63,7 +63,7 @@ end
 
 % x = 0
 if options.RMS
-    Fs = zeros(data.Nkx,data.Nky,numel(s));
+    Fs = zeros(data.Nky,data.Nkx,numel(s));
     FAM = FaM(s,xmin);
     for ip_ = 1:Np
         p_ = parray(ip_);
@@ -74,7 +74,7 @@ if options.RMS
             HLF = HH.*LL.*FAM;
             for ikx = 1:data.Nkx
                 for iky = 1:data.Nky
-                    Fs(ikx,iky,:) = squeeze(Fs(ikx,iky,:))' + Napj_(ip_,ij_,ikx,iky)*HLF;
+                    Fs(iky,ikx,:) = squeeze(Fs(iky,ikx,:))' + Napj_(ip_,ij_,iky,ikx)*HLF;
                 end
             end
        end
@@ -95,7 +95,7 @@ end
 
 % s = 0
 if options.RMS
-    Fx = zeros(data.Nkx,data.Nky,numel(x));
+    Fx = zeros(data.Nky,data.Nkx,numel(x));
     FAM = FaM(x,smin);
     for ip_ = 1:Np
         p_ = parray(ip_);
@@ -106,7 +106,7 @@ if options.RMS
             HLF = HH.*LL.*FAM;
             for ikx = 1:data.Nkx
                 for iky = 1:data.Nky
-                    Fx(ikx,iky,:) = squeeze(Fx(ikx,iky,:))' + Napj_(ip_,ij_,ikx,iky)*HLF;
+                    Fx(iky,ikx,:) = squeeze(Fx(iky,ikx,:))' + Napj_(ip_,ij_,iky,ikx)*HLF;
                 end
             end
        end
diff --git a/matlab/load/load_pjz_data.m b/matlab/load/load_pjz_data.m
new file mode 100644
index 00000000..81eea09d
--- /dev/null
+++ b/matlab/load/load_pjz_data.m
@@ -0,0 +1,16 @@
+function [ data, time, dt ] = load_pjz_data( filename,variablename,specie)
+    time     = h5read(filename,'/data/var3d/time');
+    p        = h5read(filename,['/data/grid/coordp_',specie]);
+    j        = h5read(filename,['/data/grid/coordj_',specie]);
+    z        = h5read(filename,'/data/grid/coordz');
+    dt    = h5readatt(filename,'/data/input','dt');
+    cstart= h5readatt(filename,'/data/input','start_iframe3d'); 
+    
+    data     = zeros(numel(p),numel(j),numel(z),numel(time));
+    for it = 1:numel(time)
+        tmp         = h5read(filename,['/data/var3d/',variablename,'/', num2str(cstart+it,'%06d')]);
+        data(:,:,:,it) = tmp;
+    end
+
+end
+
diff --git a/matlab/plot/plot_ballooning.m b/matlab/plot/plot_ballooning.m
index f209aea5..49d59248 100644
--- a/matlab/plot/plot_ballooning.m
+++ b/matlab/plot/plot_ballooning.m
@@ -1,30 +1,30 @@
 function [FIG] = plot_ballooning(data,options)
-    
+    FIG.fig = figure; iplot = 1;
     [~,it] = min(abs(data.Ts3D - options.time_2_plot));
+    [~,ikyarray] = min(abs(data.ky - options.kymodes));
     phi_real=(real(data.PHI(:,:,:,it)));
     phi_imag=(imag(data.PHI(:,:,:,it)));
     % Apply baollooning tranform
-    for iky=2
+    for iky=ikyarray
         dims = size(phi_real);
-        phib_real = zeros(  dims(1)*dims(3)  ,1);
+        phib_real = zeros(  dims(2)*dims(3)  ,1);
         phib_imag= phib_real;
         b_angle = phib_real;
 
-        midpoint = floor((dims(1)*dims(3) )/2)+1;
-        for ip =1: dims(1)
-            start_ =  (ip -1)*dims(3) +1;
-            end_ =  ip*dims(3);
-            phib_real(start_:end_)  = phi_real(ip,iky,:);
-            phib_imag(start_:end_)  = phi_imag(ip,iky, :);
+        for ikx =1: dims(2)
+            start_ =  (ikx -1)*dims(3) +1;
+            end_ =  ikx*dims(3);
+            phib_real(start_:end_)  = phi_real(iky,ikx,:);
+            phib_imag(start_:end_)  = phi_imag(iky,ikx,:);
         end
 
         % Define ballooning angle
         Nkx = numel(data.kx)-1; coordz = data.z;
         idx = -Nkx:1:Nkx;
-        for ip =1: dims(1)
+        for ikx =1: dims(2)
             for iz=1:dims(3)
-                ii = dims(3)*(ip -1) + iz;
-                b_angle(ii) =coordz(iz) + 2*pi*idx(ip);
+                ii = dims(3)*(ikx -1) + iz;
+                b_angle(ii) =coordz(iz) + 2*pi*idx(ikx);
             end
         end
         
@@ -40,14 +40,15 @@ function [FIG] = plot_ballooning(data,options)
         phib_real_norm  = real( phib_norm);%phib_real(:)/phib_real(idxLFS);
         phib_imag_norm  = imag( phib_norm);%phib_imag(:)/ phib_imag(idxLFS);
 
-
-        FIG.fig = figure; hold all;
-        plot(b_angle/pi, phib_real_norm,'-b');
+        subplot(numel(ikyarray),1,iplot)
+        plot(b_angle/pi, phib_real_norm,'-b'); hold on;
         plot(b_angle/pi, phib_imag_norm ,'-r');
         plot(b_angle/pi, sqrt(phib_real_norm .^2 + phib_imag_norm.^2),'-k');
         legend('real','imag','norm')
         xlabel('$\chi / \pi$')
         ylabel('$\phi_B (\chi)$');
-        title(['HeLaZ ballooning, t=',num2str(data.Ts3D(it))]);
+        title(['ky=',sprintf('%1.1f',data.ky(iky)),...
+               ',t=',sprintf('%1.1f',data.Ts3D(it))]);
+        iplot = iplot + 1;
     end
 end
diff --git a/matlab/plot/plot_fa.m b/matlab/plot/plot_fa.m
index ab93497e..00a2ca2f 100644
--- a/matlab/plot/plot_fa.m
+++ b/matlab/plot/plot_fa.m
@@ -26,13 +26,15 @@ if OPTIONS.ONED
         end
 
 else
-    [SS,XX,FFa] = compute_fa_2D(DATA, OPTIONS); 
+    [SS,XX,FFa] = compute_fa_2D(DATA, OPTIONS);  sz = size(SS);
     [~,it] = min(abs(OPTIONS.T-DATA.Ts5D)); 
         switch OPTIONS.PLT_FCT
             case 'contour'
-            contour(SS,XX,FFa',128);
+            contour(SS,XX,FFa',sum(sz)/2);
             case 'pcolor'
             pclr = pcolor(SS,XX,FFa'); set(pclr, 'edgecolor','none'); shading interp
+            case 'contourf'
+            contourf(SS,XX,FFa',sum(sz)/2);                
         end
         xlabel('$v_\parallel$'); ylabel('$\mu$');
         legend(['$\langle |f_',OPTIONS.SPECIE,'|^2\rangle_{xy}^{1/2}$',zcomp])
diff --git a/matlab/plot/plot_fa_gene.m b/matlab/plot/plot_fa_gene.m
index 5a228a1f..997cfe1f 100644
--- a/matlab/plot/plot_fa_gene.m
+++ b/matlab/plot/plot_fa_gene.m
@@ -5,8 +5,8 @@ specie = OPTIONS.specie;
 PLT_FCT= OPTIONS.PLT_FCT;
 
 file = 'coord.dat.h5';
-vp = h5read([folder,file],'/coord/vp');
-mu = h5read([folder,file],'/coord/mu');
+vp = h5read([folder,file],'/coord/vp'); nvp = numel(vp);
+mu = h5read([folder,file],'/coord/mu'); nmu = numel(mu);
 z  = h5read([folder,file],'/coord/z');
 [XX,SS] = meshgrid(mu,vp);
 
@@ -68,9 +68,9 @@ switch specie
     FFa    = abs(FFa)./max(max(abs(FFa)));
     switch PLT_FCT
         case 'contour'
-            contour(SS,XX,FFa,128);
+            contour(SS,XX,FFa);
         case 'contourf'
-            pclr = contourf(SS,XX,FFa,128); set(pclr, 'edgecolor','none')
+            pclr = contourf(SS,XX,FFa);
         case 'pcolor'
             pclr = pcolor(SS,XX,FFa); set(pclr, 'edgecolor','none'); shading interp
     end
@@ -80,9 +80,9 @@ switch specie
     FFa    = abs(FFa)./max(max(abs(FFa)));
     switch PLT_FCT
         case 'contour'
-            contour(SS,XX,FFa,128);
+            contour(SS,XX,FFa,(nvp+nmu)/2);
         case 'contourf'
-            pclr = contourf(SS,XX,FFa,128); set(pclr, 'edgecolor','none')
+            pclr = contourf(SS,XX,FFa,(nvp+nmu)/2);
         case 'pcolor'
             pclr = pcolor(SS,XX,FFa); set(pclr, 'edgecolor','none'); shading interp
     end
diff --git a/matlab/plot/show_napjz.m b/matlab/plot/show_napjz.m
new file mode 100644
index 00000000..e521b944
--- /dev/null
+++ b/matlab/plot/show_napjz.m
@@ -0,0 +1,136 @@
+function [ FIGURE ] = show_napjz( DATA, OPTIONS )
+fname = DATA.outfilenames{OPTIONS.JOBNUM+1};
+switch OPTIONS.specie
+    case 'i'
+        Pa = DATA.Pi;
+        Ja = DATA.Ji;
+        Napjz = load_pjz_data(fname,'Nipjz','i');
+        name = 'N_i^{pj}';
+        FIGURE.FIGNAME = ['Nipj_spectrum_',DATA.PARAMS];
+    case 'e'
+        Pa = DATA.Pe;
+        Ja = DATA.Je;
+        Napjz = load_pjz_data(fname,'Nipjz','i');
+        name = 'N_e^{pj}';
+        FIGURE.FIGNAME = ['Nepj_spectrum_',DATA.PARAMS];
+end
+
+switch OPTIONS.compz
+    case 'avg'
+        Napjz = squeeze(mean(Napjz,3));
+    otherwise
+        Napjz = squeeze(Napjz(:,:,OPTIONS.compz));
+end
+
+FIGURE.fig = figure; 
+set(gcf, 'Position',  [100 50 1000 400])
+
+switch OPTIONS.PLOT_TYPE
+    case 'space-time'
+    [JJ,PP] = meshgrid(Ja,Pa);
+    P2Ja = PP + 2*JJ;
+    % form an axis of velocity ordered moments
+    p2ja = unique(P2Ja);
+    % weights to average
+    weights = zeros(size(p2ja));
+    % space time of moments amplitude wrt to p+2j degree
+    Na_ST = zeros(numel(p2ja),numel(DATA.Ts5D));
+    % fill the st diagramm by averaging every p+2j=n moments together
+    for ip = 1:numel(Pa)
+        for ij = 1:numel(Ja)
+            [~,ip2j] = min(abs(p2ja-(Pa(ip)+2*Ja(ij))));
+            Na_ST(ip2j,:) = Na_ST(ip2j,:) + transpose(squeeze(Napjz(ip,ij,:)));
+            weights(ip2j) = weights(ip2j) + 1;
+        end
+    end
+    % doing the average
+    for ip2j = 1:numel(p2ja)
+        Na_ST(ip2j,:) = Na_ST(ip2j,:)/weights(ip2j);
+    end
+    % plots
+    if OPTIONS.NORMALIZED
+    plt = @(x,ip2j) x(ip2j,:)./max(x(ip2j,:));
+    else
+    plt = @(x,ip2j) x;
+    end
+    imagesc(DATA.Ts5D,p2ja,plt(Na_ST,1:numel(p2ja))); 
+    set(gca,'YDir','normal')        
+    xlabel('$t$'); ylabel('$p+2j$')
+    title('$\langle\sum_k |',name,'|\rangle_{p+2j=const}$')
+    if DATA.K_E
+    subplot(2,1,2)
+        imagesc(DATA.Ts5D,p2je,plt(Ne_ST,1:numel(p2ja))); 
+        set(gca,'YDir','normal')
+    xlabel('$t$'); ylabel('$p+2j$')
+    title('$\langle\sum_k |N_e^{pj}|\rangle_{p+2j=const}$')
+    suptitle(DATA.param_title);
+    end
+    
+    case 'Tavg-1D'
+    t0 = OPTIONS.TIME(1);
+    t1 = OPTIONS.TIME(end);
+    [~,it0] = min(abs(t0-DATA.Ts5D));
+    [~,it1] = min(abs(t1-DATA.Ts5D));
+    Napjz = mean(Napjz(:,:,it0:it1),3);
+    if DATA.K_E
+    Napjz = mean(Napjz(:,:,it0:it1),3);
+    end
+    if numel(OPTIONS.TIME) == 1
+        TITLE=['$t=',num2str(OPTIONS.TIME),'$'];
+    else
+        TITLE=['$t\in[',num2str(t0),',',num2str(t1),']$'];
+    end 
+    Napjz = squeeze(Napjz);
+
+    ymini = min(Napjz); ymaxi = max(Napjz);
+
+    if DATA.K_E
+    Napjz = squeeze(Napjz);
+    ymine = min(Napjz); ymaxe = max(Napjz);
+    ymax  = max([ymaxi ymaxe]);
+    ymin  = min([ymini ymine]);
+    end
+    if DATA.K_E
+    subplot(121)
+    end
+    if ~OPTIONS.P2J
+        for ij = 1:numel(Ja)
+            name = ['$j=',num2str(Ja(ij)),'$'];
+            semilogy(Pa,Napjz(:,ij),'o-','DisplayName',name); hold on;
+        end
+        xlabel('$p$'); 
+    else
+        for ij = 1:numel(Ja)
+            name = ['$j=',num2str(Ja(ij)),'$'];
+            semilogy(Pa+2*Ja(ij),Napjz(:,ij),'o-','DisplayName',name); hold on;
+        end
+        xlabel('$p+2j$')
+    end
+    ylabel(['$\sum_{kx,ky}|N_i^{pj}|$']);
+    legend('show');
+    title([TITLE,' He-La ion spectrum']);
+    
+    if DATA.K_E
+    subplot(122)
+    if ~OPTIONS.P2J
+        for ij = 1:numel(Ja)
+            name = ['$j=',num2str(Ja(ij)),'$'];
+            semilogy(Pa,Napjz(:,ij),'o-','DisplayName',name); hold on;
+        end
+        xlabel('p'); 
+    else
+    for ij = 1:numel(Ja)
+        name = ['$j=',num2str(Ja(ij)),'$'];
+        semilogy(Pa+2*Ja(ij),Napjz(:,ij),'o-','DisplayName',name); hold on;
+    end
+    xlabel('$p+2j$')
+    end
+    ylabel(['$\sum_{kx,ky}|N_e^{pj}|$']);
+    legend('show');
+    title([TITLE,' He-La elec. spectrum']);
+    end
+
+end
+
+end
+
diff --git a/wk/analysis_3D.m b/wk/analysis_3D.m
index 34dbcadc..f8ea45ee 100644
--- a/wk/analysis_3D.m
+++ b/wk/analysis_3D.m
@@ -40,20 +40,20 @@ end
 if 0
 %% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Options
-options.INTERP    = 1;
+options.INTERP    = 0;
 options.POLARPLOT = 0;
-options.NAME      = '\phi';
-% options.NAME      = 'N_i^{00}';
+% options.NAME      = '\phi';
+options.NAME      = 'N_i^{00}';
 % options.NAME      = 'v_y';
 % options.NAME      = 'n_i^{NZ}';
 % options.NAME      = '\Gamma_x';
 % options.NAME      = 'n_i';
-options.PLAN      = 'xy';
+options.PLAN      = 'kxky';
 % options.NAME      = 'f_e';
 % options.PLAN      = 'sx';
-options.COMP      = 'avg';
+options.COMP      = 9;
 % options.TIME      = dat.Ts5D;
-options.TIME      = 00:1:250;
+options.TIME      = 200:1:600;
 data.EPS          = 0.1;
 data.a = data.EPS * 2000;
 create_film(data,options,'.gif')
@@ -62,7 +62,7 @@ end
 if 0
 %% 2D snapshots
 % Options
-options.INTERP    = 0;
+options.INTERP    = 1;
 options.POLARPLOT = 0;
 options.AXISEQUAL = 1;
 options.NAME      = '\phi';
@@ -71,12 +71,12 @@ options.NAME      = '\phi';
 % options.NAME      = 'T_i';
 % options.NAME      = '\Gamma_x';
 % options.NAME      = 'k^2n_e';
-options.PLAN      = 'xy';
+options.PLAN      = 'xz';
 % options.NAME      = 'f_e';
 % options.PLAN      = 'sx';
-options.COMP      = 1;
-options.TIME      = [150];
-data.a = data.EPS * 1000;
+options.COMP      = 'avg';
+options.TIME      = [260:265];
+data.a = data.EPS * 2e3;
 fig = photomaton(data,options);
 save_figure(data,fig)
 end
@@ -85,8 +85,8 @@ if 0
 %% 3D plot on the geometry
 options.INTERP    = 1;
 options.NAME      = 'n_i';
-options.PLANES    = 1;
-options.TIME      = [100 200];
+options.PLANES    = [1:2:12];
+options.TIME      = [60];
 options.PLT_MTOPO = 1;
 data.rho_o_R      = 2e-3; % Sound larmor radius over Machine size ratio
 fig = show_geometry(data,options);
@@ -100,8 +100,8 @@ if 0
 options.SPAR      = gene_data.vp';
 options.XPERP     = gene_data.mu';
 options.Z         = 'avg';
-options.T         = 200;
-options.PLT_FCT   = 'pcolor';
+options.T         = 340;
+options.PLT_FCT   = 'contour';
 options.ONED      = 0;
 options.non_adiab = 1;
 options.SPECIE    = 'i';
@@ -113,10 +113,16 @@ end
 if 0
 %% Hermite-Laguerre spectrum
 % options.TIME = 'avg';
-options.P2J  = 1;
-options.ST   = 0;
+options.P2J        = 1;
+options.ST         = 0;
+options.PLOT_TYPE  = 'space-time';
 options.NORMALIZED = 0;
-fig = show_moments_spectrum(data,options);
+options.JOBNUM     = 03;
+options.TIME       = [200 500];
+options.specie     = 'i';
+options.compz      = 'avg';
+% fig = show_moments_spectrum(data,options);
+fig = show_napjz(data,options);
 save_figure(data,fig)
 end
 
diff --git a/wk/gene_analysis_3D.m b/wk/gene_analysis_3D.m
index 7479f55e..5716569d 100644
--- a/wk/gene_analysis_3D.m
+++ b/wk/gene_analysis_3D.m
@@ -7,7 +7,6 @@ folder = '/misc/gene_results/shearless_cyclone/s_alpha_output_1.0/';
 % folder = '/misc/gene_results/HP_fig_2b_mu_5e-2/';
 % folder = '/misc/gene_results/HP_fig_2c_mu_5e-2/';
 gene_data = load_gene_data(folder);
-% gene_data = rotate_c_plane_nxnky_to_nkxny(gene_data);
 gene_data = invert_kxky_to_kykx_gene_results(gene_data);
 if 1
 %% Space time diagramm (fig 11 Ivanov 2020)
@@ -23,100 +22,19 @@ end
 if 0
 %% 2D snapshots
 % Options
-options.INTERP    = 1;
-options.POLARPLOT = 0;
-options.AXISEQUAL = 1;
-options.NAME      = '\phi';
-% options.NAME      = 'n_i';
-% options.NAME      = 'T_i';
-% options.NAME      = '\Gamma_x';
-% options.NAME      = 'k^2n_e';
-options.PLAN      = 'xz';
-% options.NAME      ='f_e';
-% options.PLAN      = 'sx';
-options.COMP      = 'avg';
-options.TIME      = [100 300 900];
-gene_data.a = data.EPS * 2000;
-fig = photomaton(gene_data,options);
-save_figure(gene_data,fig)
-end
-
-if 0
-%% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% Options
 options.INTERP    = 0;
 options.POLARPLOT = 0;
-options.NAME      = '\phi';
-% options.NAME      = 'v_y';
-% options.NAME      = 'n_i^{NZ}';
-% options.NAME      = '\Gamma_x';
-% options.NAME      = 'n_i';
-options.PLAN      = 'xy';
-% options.NAME      = 'f_e';
-% options.PLAN      = 'sx';
-options.COMP      = 'avg';
-options.TIME      = 000:200;
-gene_data.a = data.EPS * 2000;
-create_film(gene_data,options,'.gif')
-end
-
-if 0
-%% Geometry
-names = {'$g^{xx}$','$g^{xy}$','$g^{xz}$','$g^{yy}$','$g^{yz}$','$g^{zz}$',...
-         '$B_0$','$\partial_x B_0$','$\partial_y B_0$','$\partial_z B_0$',...
-         '$J$','$R$','$\phi$','$Z$','$\partial_R x$','$\partial_Z x$'};
-figure;
-subplot(311)
-    for i = 1:6
-    plot(gene_data.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on;
-    end
-    xlim([min(gene_data.z),max(gene_data.z)]); legend('show'); title('GENE geometry');
-    
-subplot(312)
-    for i = 7:10
-    plot(gene_data.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on;
-    end
-    xlim([min(gene_data.z),max(gene_data.z)]); legend('show');
-
-subplot(313)
-    for i = 11:16
-    plot(gene_data.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on;
-    end
-    xlim([min(gene_data.z),max(gene_data.z)]); legend('show');
-
-end
-% folder = '/misc/gene_results/HP_fig_2b_mu_5e-2/';
-folder = '/misc/gene_results/shearless_cyclone/s_alpha_output_1.0/';
-% folder = '/misc/gene_results/cyclone/s_alpha_output_1.0/';
-gene_data = load_gene_data(folder);
-gene_data = rotate_c_plane_nxnky_to_nkxny(gene_data);
-if 1
-%% Space time diagramm (fig 11 Ivanov 2020)
-TAVG_0 = 0.8*gene_data.Ts3D(end); TAVG_1 = gene_data.Ts3D(end); % Averaging times duration
-options.NMVA     = 1;              % Moving average for time traces
-% options.ST_FIELD = '\Gamma_x';          % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
-options.ST_FIELD = '\phi';          % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
-options.INTERP   = 1;
-fig = plot_radial_transport_and_spacetime(gene_data,options);
-% save_figure(data,fig)
-end
-
-if 0
-%% 2D snapshots
-% Options
-options.INTERP    = 1;
-options.POLARPLOT = 0;
 options.AXISEQUAL = 1;
-options.NAME      = '\phi';
-% options.NAME      = 'n_i';
+% options.NAME      = 'Q_x';
+options.NAME      = 'n_i';
 % options.NAME      = 'T_i';
 % options.NAME      = '\Gamma_x';
 % options.NAME      = 'k^2n_e';
-options.PLAN      = 'xz';
+options.PLAN      = 'xy';
 % options.NAME      ='f_e';
 % options.PLAN      = 'sx';
 options.COMP      = 'avg';
-options.TIME      = [100 300 900];
+options.TIME      = [200];
 gene_data.a = data.EPS * 2000;
 fig = photomaton(gene_data,options);
 save_figure(gene_data,fig)
@@ -132,11 +50,11 @@ options.NAME      = '\phi';
 % options.NAME      = 'n_i^{NZ}';
 % options.NAME      = '\Gamma_x';
 % options.NAME      = 'n_i';
-options.PLAN      = 'xy';
+options.PLAN      = 'kxky';
 % options.NAME      = 'f_e';
 % options.PLAN      = 'sx';
 options.COMP      = 'avg';
-options.TIME      = 000:170;
+options.TIME      = 000:700;
 gene_data.a = data.EPS * 2000;
 create_film(gene_data,options,'.gif')
 end
@@ -175,7 +93,7 @@ options.PLT_FCT = 'pcolor';
 options.folder  = folder;
 options.Z       = 'avg';
 options.FIELD   = '<f_>';
-options.ONED    = 0;
+options.ONED    = 1;
 % options.FIELD   = 'Q_es';
 plot_fa_gene(options);
 end
diff --git a/wk/quick_run.m b/wk/quick_run.m
index e7525829..b3f5829b 100644
--- a/wk/quick_run.m
+++ b/wk/quick_run.m
@@ -13,33 +13,33 @@ EXECNAME = 'helaz3.12';
 CLUSTER.TIME  = '99:00:00'; % allocation time hh:mm:ss
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% PHYSICAL PARAMETERS
-NU      = 0.1;   % Collision frequency
+NU      = 0.0;   % Collision frequency
 TAU     = 1.0;    % e/i temperature ratio
-K_N     = 1.9;%2.22;   % Density gradient drive
-K_T     = 0.25*K_N;   % Temperature '''
+K_N     = 2.22;   % Density gradient drive
+K_T     = 6.96;   % Temperature '''
 K_E     = 0.0;   % Electrostat '''
 SIGMA_E = 0.0233380;   % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
-KIN_E   = 1;     % 1: kinetic electrons, 2: adiabatic electrons
+KIN_E   = 0;     % 1: kinetic electrons, 2: adiabatic electrons
 %% GRID PARAMETERS
-PMAXE   = 4;     % Hermite basis size of electrons
-JMAXE   = 2;     % Laguerre "
-PMAXI   = 4;     % " ions
-JMAXI   = 2;     % "
-NX      = 32;    % real space x-gridpoints
-NY      = 32;     %     ''     y-gridpoints
+PMAXE   = 8;     % Hermite basis size of electrons
+JMAXE   = 4;     % Laguerre "
+PMAXI   = 8;     % " ions
+JMAXI   = 4;     % "
+NX      = 1;    % real space x-gridpoints
+NY      = 20;     %     ''     y-gridpoints
 LX      = 120;   % Size of the squared frequency domain
 LY      = 120;     % Size of the squared frequency domain
-NZ      = 1;     % number of perpendicular planes (parallel grid)
+NZ      = 16;     % 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
 % GEOMETRY= 's-alpha';
-Q0      = 2.5;    % safety factor
+Q0      = 1.4;    % safety factor
 SHEAR   = 0.0;    % magnetic shear (Not implemented yet)
 EPS     = 0.18;    % inverse aspect ratio
 %% TIME PARMETERS
-TMAX    = 25;  % Maximal time unit
+TMAX    = 50;  % Maximal time unit
 DT      = 2e-2;   % Time step
 SPS0D   = 1;      % Sampling per time unit for 2D arrays
 SPS2D   = 0;      % Sampling per time unit for 2D arrays
@@ -49,14 +49,14 @@ SPSCP   = 0;    % Sampling per time unit for checkpoints
 JOB2LOAD= -1;
 %% OPTIONS
 SIMID   = 'quick_run';  % Name of the simulation
-LINEARITY = 'nonlinear';   % activate non-linearity (is cancelled if KXEQ0 = 1)
+LINEARITY = 'linear';   % activate non-linearity (is cancelled if KXEQ0 = 1)
 % Collision operator
 % (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
-CO      = 'SG';
+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: gyrofluid closure (p+2j<=Pmax))s
+CLOS    = 1;   % Closure model (0: =0 truncation, 1: gyrofluid closure (p+2j<=Pmax))s
 NL_CLOS = 0;   % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS)
 KERN    = 0;   % Kernel model (0 : GK)
 INIT_OPT= 'phi';   % Start simulation with a noisy mom00/phi/allmom
@@ -75,7 +75,7 @@ MU      = 0.0; % Hyperdiffusivity coefficient
 INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
 MU_X    = MU;     %
 MU_Y    = MU;     %
-MU_Z    = 0.0;     %
+MU_Z    = 0.2;     %
 MU_P    = 0.0;     %
 MU_J    = 0.0;     %
 LAMBDAD = 0.0;
@@ -114,12 +114,13 @@ end
 if 0
 %% Ballooning plot
 options.time_2_plot = data.Ts3D(end);
+options.kymodes     = [0.2 0.5 0.7];
 options.normalized  = 1;
 options.field       = 'phi';
 fig = plot_ballooning(data,options);
 end
 
-if 1
+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
-- 
GitLab