From 4edaa07d6309af460c99711b3d73be3f5e9b8d2d Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Tue, 27 Oct 2020 10:47:47 +0100
Subject: [PATCH] Separation between loading results and post processing

---
 wk/analysis_2D.m  | 240 ++++++++++++++++++++++++++--------------------
 wk/load_results.m |  23 +++++
 2 files changed, 159 insertions(+), 104 deletions(-)
 create mode 100644 wk/load_results.m

diff --git a/wk/analysis_2D.m b/wk/analysis_2D.m
index 6d23ad40..2a379eda 100644
--- a/wk/analysis_2D.m
+++ b/wk/analysis_2D.m
@@ -1,34 +1,14 @@
-%% load results %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-JOBNUM = 00;
-FMT = '.fig';
-filename = [BASIC.SIMID,'_','%.2d.h5'];
-filename = sprintf(filename,JOBNUM); disp(['Loading ',filename])
-% Loading from output file
-if strcmp(OUTPUTS.write_moments,'.true.') 
-[Nipj, Pi, Ji, kr, kz, Ts, dt] = load_5D_data(filename, 'moments_i');
-[Nepj, Pe, Je,  ~,  ~,  ~,  ~] = load_5D_data(filename, 'moments_e');
-Ni00    = squeeze(Nipj(1,1,:,:,:));
-Ne00    = squeeze(Nepj(1,1,:,:,:));
-else
-[Ni00, kr, kz, Ts, dt] = load_2D_data(filename, 'Ni00');
- Ne00                  = load_2D_data(filename, 'Ne00');
- Pi = [0]; Ji = Pi; Pe = Pi; Je = Pi;
- Nipj = zeros(1,1,numel(kr),numel(kz),numel(Ts));
- Nepj = Nipj;
- Nipj(1,1,:,:,:) = Ni00; Nepj(1,1,:,:,:) = Ne00;
-end
-PHI                            = load_2D_data(filename, 'phi');
+%% Load results
+load_results
 
-if strcmp(OUTPUTS.write_non_lin,'.true.') 
-    Sipj    = load_5D_data(filename, 'Sipj');
-    Sepj    = load_5D_data(filename, 'Sepj');
-end
 %% Retrieving max polynomial degree and sampling info
 Npe = numel(Pe); Nje = numel(Je); [JE,PE] = meshgrid(Je,Pe);
 Npi = numel(Pi); Nji = numel(Ji); [JI,PI] = meshgrid(Ji,Pi);
-dt_samp = mean(diff(Ts)); Ns      = numel(Ts);
+Ns5D      = numel(Ts5D);
+Ns2D      = numel(Ts2D);
 % renaming and reshaping quantity of interest
-Ts      = Ts';
+Ts5D      = Ts5D';
+Ts2D      = Ts2D';
 if strcmp(OUTPUTS.write_non_lin,'.true.')
     Si00      = squeeze(Sipj(1,1,:,:,:));
     Se00      = squeeze(Sepj(1,1,:,:,:));
@@ -54,36 +34,54 @@ ne00   = zeros(Nr,Nz);
 ni00   = zeros(Nr,Nz);
 si00   = zeros(Nr,Nz);
 phi    = zeros(Nr,Nz);
+drphi  = zeros(Nr,Nz);
+dzphi  = zeros(Nr,Nz);
 
-for it = 1:numel(PHI(1,1,:))
+for it = 1:numel(Ts2D)
     NE_ = Ne00(:,:,it); NI_ = Ni00(:,:,it); PH_ = PHI(:,:,it);
     ne00(:,:,it)  = real(fftshift(ifft2((NE_),Nr,Nz)));
     ni00(:,:,it)  = real(fftshift(ifft2((NI_),Nr,Nz)));
     phi (:,:,it)  = real(fftshift(ifft2((PH_),Nr,Nz)));
-if strcmp(OUTPUTS.write_non_lin,'.true.')
-    SI_ = Si00(:,:,it);
-    si00(:,:,it)  = fftshift(ifft2(half_2_full_cc_2D(SI_),'symmetric'));
-end
+    drphi(:,:,it) = real(fftshift(ifft2(1i*KR.*(PH_),Nr,Nz)));
+    dzphi(:,:,it) = real(fftshift(ifft2(1i*KZ.*(PH_),Nr,Nz)));
 end
 
 % Post processing
 disp('- post processing')
-E_pot    = zeros(1,Ns);      % Potential energy n^2
-E_kin    = zeros(1,Ns);      % Kinetic energy grad(phi)^2
-ExB      = zeros(1,Ns);      % ExB drift intensity \propto |\grad \phi|
-CFL      = zeros(1,Ns);      % CFL
-Ne_norm  = zeros(Npe,Nje,Ns);% Time evol. of the norm of Napj
-Ni_norm  = zeros(Npi,Nji,Ns);% .
+E_pot    = zeros(1,Ns2D);      % Potential energy n^2
+E_kin    = zeros(1,Ns2D);      % Kinetic energy grad(phi)^2
+ExB      = zeros(1,Ns2D);      % ExB drift intensity \propto |\grad \phi|
+Flux_ri  = zeros(1,Ns2D);      % Particle flux Gamma = <ni drphi>
+Flux_zi  = zeros(1,Ns2D);      % Particle flux Gamma = <ni dzphi>
+Flux_re  = zeros(1,Ns2D);      % Particle flux Gamma = <ne drphi>
+Flux_ze  = zeros(1,Ns2D);      % Particle flux Gamma = <ne dzphi>
+Ne_norm  = zeros(Npe,Nje,Ns5D);% Time evol. of the norm of Napj
+Ni_norm  = zeros(Npi,Nji,Ns5D);% .
 if strcmp(OUTPUTS.write_non_lin,'.true.')
-Se_norm  = zeros(Npe,Nje,Ns);% Time evol. of the norm of Sapj
-Si_norm  = zeros(Npi,Nji,Ns);% .
-Sne00_norm = zeros(1,Ns);    % Time evol. of the amp of e nonlin term
-Sni00_norm = zeros(1,Ns);    %
+Se_norm  = zeros(Npe,Nje,Ns5D);% Time evol. of the norm of Sapj
+Si_norm  = zeros(Npi,Nji,Ns5D);% .
+Sne00_norm = zeros(1,Ns2D);    % Time evol. of the amp of e nonlin term
+Sni00_norm = zeros(1,Ns2D);    %
 end
 
 Ddr = 1i*KR; Ddz = 1i*KZ; lapl   = Ddr.^2 + Ddz.^2; 
 
-for it = 1:numel(PHI(1,1,:))
+for it = 1:numel(Ts2D) % Loop over 2D arrays
+    NE_ = Ne00(:,:,it); NI_ = Ni00(:,:,it); PH_ = PHI(:,:,it);
+    E_pot(it)   = pi/Lr/Lz*sum(sum(abs(NI_).^2))/Nkr/Nkr; % integrate through Parseval id
+    E_kin(it)   = pi/Lr/Lz*sum(sum(abs(Ddr.*PH_).^2+abs(Ddz.*PH_).^2))/Nkr/Nkr;
+    ExB(it)     = max(max(max(abs(phi(3:end,:,it)-phi(1:end-2,:,it))/(2*dr))),max(max(abs(phi(:,3:end,it)-phi(:,1:end-2,it))'/(2*dz))));
+    Flux_ri(it)  = sum(sum(ni00(:,:,it).*drphi(:,:,it)))*dr*dz;
+    Flux_zi(it)  = sum(sum(ni00(:,:,it).*dzphi(:,:,it)))*dr*dz;
+    Flux_re(it)  = sum(sum(ne00(:,:,it).*drphi(:,:,it)))*dr*dz;
+    Flux_ze(it)  = sum(sum(ne00(:,:,it).*dzphi(:,:,it)))*dr*dz;
+end
+
+E_kin_KZ = mean(mean(abs(Ddr.*PHI(:,:,it)).^2+abs(Ddz.*PHI(:,:,it)).^2,3),2);
+E_kin_KR = mean(mean(abs(Ddr.*PHI(:,:,it)).^2+abs(Ddz.*PHI(:,:,it)).^2,3),2);
+dEdt     = diff(E_pot+E_kin)./dt2D;
+
+for it = 1:numel(Ts5D) % Loop over 5D arrays
     NE_ = Ne00(:,:,it); NI_ = Ni00(:,:,it); PH_ = PHI(:,:,it);
     Ne_norm(:,:,it)= sum(sum(abs(Nepj(:,:,:,:,it)),3),4);
     Ni_norm(:,:,it)= sum(sum(abs(Nipj(:,:,:,:,it)),3),4);
@@ -93,22 +91,17 @@ if strcmp(OUTPUTS.write_non_lin,'.true.')
     Si_norm(:,:,it)= sum(sum(abs(Sipj(:,:,:,:,it)),3),4);
     Sni00_norm(it) = sum(sum(abs(Si00(:,:,it))));
 end
-    E_pot(it)   = pi/Lr/Lz*sum(sum(abs(NI_).^2))/Nkr/Nkr; % integrate through Parseval id
-    E_kin(it)   = pi/Lr/Lz*sum(sum(abs(Ddr.*PH_).^2+abs(Ddz.*PH_).^2))/Nkr/Nkr;
-    ExB(it)     = max(max(max(abs(phi(3:end,:,it)-phi(1:end-2,:,it))/(2*dr))),max(max(abs(phi(:,3:end,it)-phi(:,1:end-2,it))'/(2*dz))));
 end
-E_kin_KZ = mean(mean(abs(Ddr.*PHI(:,:,it)).^2+abs(Ddz.*PHI(:,:,it)).^2,3),2);
-E_kin_KR = mean(mean(abs(Ddr.*PHI(:,:,it)).^2+abs(Ddz.*PHI(:,:,it)).^2,3),2);
-dEdt     = diff(E_pot+E_kin)./diff(Ts);
+
 
 %% Growth rate
 disp('- growth rate')
-tend   = Ts(end); tstart   = 0.8*tend; 
+tend   = Ts2D(end); tstart   = 0.6*tend; 
 g_          = zeros(Nkr,Nkz);
 [~,ikr0KH] = min(abs(kr-KR0KH));
 for ikr = 1:Nkr
     for ikz = 1:Nkz
-        g_(ikr,ikz) = LinearFit_s(Ts,squeeze(abs(Ni00(ikr,ikz,:))),tstart,tend);
+        g_(ikr,ikz) = LinearFit_s(Ts2D,squeeze(abs(Ni00(ikr,ikz,:))),tstart,tend);
     end
 end
 % gkr0kz_Ni00 = max(real(g_(:,:)),[],1);
@@ -123,7 +116,7 @@ fig = figure; FIGNAME = ['t_evolutions',sprintf('_%.2d',JOBNUM)];
         for ij = 1:Nje
             plt      = @(x) squeeze(x(ip,ij,:));
             plotname = ['$N_e^{',num2str(ip-1),num2str(ij-1),'}$'];
-            semilogy(Ts,plt(Ne_norm),'DisplayName',plotname); hold on;
+            semilogy(Ts5D,plt(Ne_norm),'DisplayName',plotname); hold on;
         end
     end
     grid on; xlabel('$t$'); ylabel('$\sum_{k_r,k_z}|N_e^{pj}|$');
@@ -132,21 +125,23 @@ fig = figure; FIGNAME = ['t_evolutions',sprintf('_%.2d',JOBNUM)];
         for ij = 1:Nji
             plt      = @(x) squeeze(x(ip,ij,:));
             plotname = ['$N_i^{',num2str(ip-1),num2str(ij-1),'}$'];
-            semilogy(Ts,plt(Ni_norm),'DisplayName',plotname); hold on;
+            semilogy(Ts5D,plt(Ni_norm),'DisplayName',plotname); hold on;
         end
     end
     grid on; xlabel('$t$'); ylabel('$\sum_{k_r,k_z}|N_i^{pj}|$');
     subplot(223)
-        semilogy(Ts,E_kin+E_pot,'-','DisplayName','$\sum|ik\tilde\phi_i|^2+\sum|\tilde n_i|^2$')
-        hold on;
-        grid on; xlabel('$t$'); ylabel('$E$'); %legend('show');
+        plot(Ts2D,Flux_ri,'-','DisplayName','$\Gamma_{ri}$'); hold on;
+        plot(Ts2D,Flux_zi,'-','DisplayName','$\Gamma_{zi}$'); hold on;
+        plot(Ts2D,Flux_re,'-','DisplayName','$\Gamma_{re}$')
+        plot(Ts2D,Flux_ze,'-','DisplayName','$\Gamma_{ze}$')
+        grid on; xlabel('$t$'); ylabel('$\Gamma$'); %legend('show');
 if strcmp(OUTPUTS.write_non_lin,'.true.')
     subplot(224)
     for ip = 1:Npi
         for ij = 1:Nji
             plt      = @(x) squeeze(x(ip,ij,:));
             plotname = ['$S_i^{',num2str(ip-1),num2str(ij-1),'}$'];
-            semilogy(Ts,plt(Si_norm),'DisplayName',plotname); hold on;
+            semilogy(Ts5D,plt(Si_norm),'DisplayName',plotname); hold on;
         end
     end
     grid on; xlabel('$t$'); ylabel('$S$'); %legend('show');
@@ -165,53 +160,90 @@ save_figure
 t0    = 0;
 skip_ = 1; 
 DELAY = 0.01*skip_;
-FRAMES = floor(t0/dt_samp)+1:skip_:numel(Ts);
+FRAMES = floor(t0/dt2D)+1:skip_:numel(Ts2D);
 if 0
 %% GIFS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% Density ion
-GIFNAME = ['ni00',sprintf('_%.2d',JOBNUM)]; INTERP = 1;
-FIELD = real(ni00); X = RR; Y = ZZ; T = Ts;
+GIFNAME = ['ni00',sprintf('_%.2d',JOBNUM)]; INTERP = 0;
+FIELD = real(ni00); X = RR; Y = ZZ; T = Ts2D;
 FIELDNAME = '$n_i^{00}$'; XNAME = '$r$'; YNAME = '$z$';
 create_gif
 %% Density electron
 GIFNAME = ['ne00',sprintf('_%.2d',JOBNUM)]; INTERP = 1;
-FIELD = real(ne00); X = RR; Y = ZZ; T = Ts;
+FIELD = real(ne00); X = RR; Y = ZZ; T = Ts2D;
 FIELDNAME = '$n_e^{00}$'; XNAME = '$r$'; YNAME = '$z$';
 create_gif
 %% Phi
 GIFNAME = ['phi',sprintf('_%.2d',JOBNUM)];INTERP = 1;
-FIELD = real(phi); X = RR; Y = ZZ; T = Ts;
+FIELD = real(phi); X = RR; Y = ZZ; T = Ts2D;
 FIELDNAME = '$\phi$'; XNAME = '$r$'; YNAME = '$z$';
 create_gif
 %% Density ion frequency
 GIFNAME = ['Ni00',sprintf('_%.2d',JOBNUM)]; INTERP = 0;
-FIELD =ifftshift((abs(Ni00)),2); X = fftshift(KR,2); Y = fftshift(KZ,2); T = Ts;
+FIELD =ifftshift((abs(Ni00)),2); X = fftshift(KR,2); Y = fftshift(KZ,2); T = Ts2D;
 FIELDNAME = '$N_i^{00}$'; XNAME = '$k_r$'; YNAME = '$k_z$';
 create_gif
-%% PJ ion moment frequency
-p_ = 0+1; j_ = 1+1;
-GIFNAME = ['Ni',num2str(p_),num2str(j_),sprintf('_%.2d',JOBNUM)]; INTERP = 0;
-FIELD =ifftshift(abs(squeeze(Nipj(p_,j_,:,:,:))),2); X = fftshift(KR,2); Y = fftshift(KZ,2); T = Ts;
-FIELDNAME = ['$N_i^{',num2str(p_-1),num2str(j_-1),'}$']; XNAME = '$k_r$'; YNAME = '$k_z$';
-create_gif
-%% non linear term frequ. space
-GIFNAME = ['Si00',sprintf('_%.2d',JOBNUM)]; INTERP = 0;
-FIELD = fftshift((abs(Si00)),2); X = fftshift(KR,2); Y = fftshift(KZ,2); T = Ts;
-FIELDNAME = '$S_i^{00}$'; XNAME = '$k_r$'; YNAME = '$k_z$';
-create_gif
-%% non linear term real space
-GIFNAME = ['si00',sprintf('_%.2d',JOBNUM)]; INTERP = 1;
-FIELD = real(si00); X = RR; Y = ZZ; T = Ts;
-FIELDNAME = '$s_i^{00}$'; XNAME = '$r$'; YNAME = '$z$';
-create_gif
-%% Electron mode growth
-GIFNAME = ['norm_Nepj',sprintf('_%.2d',JOBNUM)]; INTERP = 1;
-FIELD = Ne_norm; X = PE; Y = JE; T = Ts;
-FIELDNAME = '$\max_k{\gamma_e}$'; XNAME = '$p$'; YNAME = '$j$';
-create_gif
+%% Density ion frequency @ kr = 0
+GIFNAME = ['Ni00_kr0',sprintf('_%.2d',JOBNUM)]; INTERP = 0;
+FIELD =(squeeze(abs(Ni00(1,:,:)))); linestyle = 'o-.';
+X = (kz); T = Ts2D; YMIN = -.1; YMAX = 1.1; XMIN = min(kz); XMAX = max(kz);
+FIELDNAME = '$N_i^{00}(kr=0)$'; XNAME = '$k_r$';
+create_gif_1D
+% %% PJ ion moment frequency
+% p_ = 0+1; j_ = 3+1;
+% GIFNAME = ['Ni',num2str(p_),num2str(j_),sprintf('_%.2d',JOBNUM)]; INTERP = 0;
+% FIELD =ifftshift(abs(squeeze(Nipj(p_,j_,:,:,:))),2); X = fftshift(KR,2); Y = fftshift(KZ,2); T = Ts;
+% FIELDNAME = ['$N_i^{',num2str(p_-1),num2str(j_-1),'}$']; XNAME = '$k_r$'; YNAME = '$k_z$';
+% create_gif
+% %% non linear term frequ. space
+% GIFNAME = ['Si00',sprintf('_%.2d',JOBNUM)]; INTERP = 0;
+% FIELD = fftshift((abs(Si00)),2); X = fftshift(KR,2); Y = fftshift(KZ,2); T = Ts;
+% FIELDNAME = '$S_i^{00}$'; XNAME = '$k_r$'; YNAME = '$k_z$';
+% create_gif
+% %% Electron mode growth
+% GIFNAME = ['norm_Nepj',sprintf('_%.2d',JOBNUM)]; INTERP = 1;
+% FIELD = Ne_norm; X = PE; Y = JE; T = Ts;
+% FIELDNAME = '$\max_k{\gamma_e}$'; XNAME = '$p$'; YNAME = '$j$';
+% create_gif
 end
 %%
 
+if 0
+%% Asymmetry error on kr = 0
+up   = 2:Nkz/2;
+down = Nkz:-1:Nkz/2+2;
+eps_asym_r = Ts5D;
+eps_asym_i = Ts5D;
+for it = 1:numel(Ts5D)
+    eps_asym_r(it) = max(squeeze(abs(real(Ni00(1,up,it))-real(Ni00(1,down,it)))));
+    eps_asym_i(it) = max(abs(imag(Ni00(1,up,it))+imag(Ni00(1,down,it))));
+end
+
+figure
+    it = 2;
+    subplot(311)
+        plot(kz(up),real(Ni00(1,up,it))); hold on
+        plot(kz(up),real(Ni00(1,down,it)),'--')
+        ylabel('Real'); title(['$N_i^{00}(t=',num2str(Ts5D(it)),')$'])
+    subplot(312)
+        plot(kz(up),imag(Ni00(1,up,it))); hold on
+        plot(kz(up),-imag(Ni00(1,down,it)),'--')
+        ylabel('Imag'); xlabel('$kz,kr=0$')
+    subplot(313)
+        plot(kz(2:Nkz/2),real(Ni00(1,up,it))-real(Ni00(1,down,it))); hold on
+        plot(kz(2:Nkz/2),imag(Ni00(1,up,it))+imag(Ni00(1,down,it)),'--')
+        legend('Re','Im');
+        ylabel('err'); xlabel('$kz,kr=0$')
+%%       
+figure
+     semilogy(Ts5D,eps_asym_r); hold on
+     semilogy(Ts5D,eps_asym_i,'--');
+     title('$\max_{k_z}|N_i^{00}(0,k_z)-N_i^{00}(0,-k_z)|$');
+     legend('Re','Im'); grid on; xlabel('$t$')
+
+end
+
+
 if 0
 %% Growth rate
 fig = figure; FIGNAME = ['growth_rate',sprintf('_%.2d',JOBNUM)];
@@ -223,7 +255,7 @@ save_figure
 
 end
 %%
-if 1
+if 0
 %% Mode time evolution
 [~,ik00] = min(abs(kz));
 [~,idk]  = min(abs(kz-dkz));
@@ -231,23 +263,23 @@ if 1
 [~,ik75] = min(abs(kz-0.7*KR0KH));
 [~,ik10] = min(abs(kz-1.00*KR0KH));
 plt = @(x) abs(squeeze(x));
-fig = figure; FIGNAME = ['frame',sprintf('_%.2d',JOBNUM)];
-        semilogy(Ts,plt(Ni00(ikr0KH,ik00,:)),'DisplayName', ...
+fig = figure; FIGNAME = ['mode_time_evolution',sprintf('_%.2d',JOBNUM)];
+        semilogy(Ts2D,plt(Ni00(ikr0KH,ik00,:)),'DisplayName', ...
             ['$k_z/k_{r0} = $',num2str(kz(ik00)/KR0KH)]); hold on
-        semilogy(Ts,plt(Ni00(ikr0KH,idk,:)),'DisplayName', ...
+        semilogy(Ts2D,plt(Ni00(ikr0KH,idk,:)),'DisplayName', ...
             ['$k_z/k_{r0} = $',num2str(kz(idk)/KR0KH)]); hold on
-        semilogy(Ts,plt(Ni00(ikr0KH,ik50,:)),'DisplayName', ...
+        semilogy(Ts2D,plt(Ni00(ikr0KH,ik50,:)),'DisplayName', ...
             ['$k_z/k_{r0} = $',num2str(kz(ik50)/KR0KH)]); hold on
-        semilogy(Ts,plt(Ni00(ikr0KH,ik75,:)),'DisplayName', ...
+        semilogy(Ts2D,plt(Ni00(ikr0KH,ik75,:)),'DisplayName', ...
             ['$k_z/k_{r0} = $',num2str(kz(ik75)/KR0KH)]); hold on
-        semilogy(Ts,plt(Ni00(ikr0KH,ik10,:)),'DisplayName', ...
+        semilogy(Ts2D,plt(Ni00(ikr0KH,ik10,:)),'DisplayName', ...
             ['$k_z/k_{r0} = $',num2str(kz(ik10)/KR0KH)]); hold on
         xlabel('$t$'); ylabel('$\hat n_i^{00}$'); legend('show');
-        semilogy(Ts,plt(Ni00(ikr0KH,ik00,end))*exp(gkr0kz_Ni00(ik00)*(Ts-Ts(end))),'k--')
-        semilogy(Ts,plt(Ni00(ikr0KH,idk,end))*exp(gkr0kz_Ni00(idk)*(Ts-Ts(end))),'k--')
-        semilogy(Ts,plt(Ni00(ikr0KH,ik50,end))*exp(gkr0kz_Ni00(ik50)*(Ts-Ts(end))),'k--')
-        semilogy(Ts,plt(Ni00(ikr0KH,ik75,end))*exp(gkr0kz_Ni00(ik75)*(Ts-Ts(end))),'k--')
-        semilogy(Ts,plt(Ni00(ikr0KH,ik10,end))*exp(gkr0kz_Ni00(ik10)*(Ts-Ts(end))),'k--')
+        semilogy(Ts2D,plt(Ni00(ikr0KH,ik00,end))*exp(gkr0kz_Ni00(ik00)*(Ts2D-Ts2D(end))),'k--')
+        semilogy(Ts2D,plt(Ni00(ikr0KH,idk,end))*exp(gkr0kz_Ni00(idk)*(Ts2D-Ts2D(end))),'k--')
+        semilogy(Ts2D,plt(Ni00(ikr0KH,ik50,end))*exp(gkr0kz_Ni00(ik50)*(Ts2D-Ts2D(end))),'k--')
+        semilogy(Ts2D,plt(Ni00(ikr0KH,ik75,end))*exp(gkr0kz_Ni00(ik75)*(Ts2D-Ts2D(end))),'k--')
+        semilogy(Ts2D,plt(Ni00(ikr0KH,ik10,end))*exp(gkr0kz_Ni00(ik10)*(Ts2D-Ts2D(end))),'k--')
         plot(tstart*[1 1],ylim,'k-','LineWidth',0.5);
         plot(tend*[1 1],ylim,'k-','LineWidth',0.5);
 save_figure
@@ -256,21 +288,21 @@ end
 %%
 if 0
 %% Show frame in real space
-tf = 20; [~,it] = min(abs(Ts-tf));
+tf = 20; [~,it] = min(abs(Ts5D-tf));
 fig = figure; FIGNAME = ['rz_frame',sprintf('_%.2d',JOBNUM)];
     subplot(221); plt = @(x) (((x)));
         pclr = pcolor((RR),(ZZ),plt(ne00(:,:,it))); set(pclr, 'edgecolor','none'); colorbar;
-        xlabel('$r$'); ylabel('$z$'); title(sprintf('t=%.3d',Ts(it))); legend('$|\hat n_e^{00}|$');
+        xlabel('$r$'); ylabel('$z$'); title(sprintf('t=%.3d',Ts5D(it))); legend('$|\hat n_e^{00}|$');
     subplot(222); plt = @(x) ((x));
         pclr = pcolor((RR),(ZZ),plt(ni00(:,:,it))); set(pclr, 'edgecolor','none'); colorbar;
-        xlabel('$r$'); ylabel('$z$'); title(sprintf('t=%.3d',Ts(it))); legend('$|\hat n_i^{00}|$');
+        xlabel('$r$'); ylabel('$z$'); title(sprintf('t=%.3d',Ts5D(it))); legend('$|\hat n_i^{00}|$');
     subplot(223); plt = @(x) ((x));
         pclr = pcolor((RR),(ZZ),plt(phi(:,:,it))); set(pclr, 'edgecolor','none'); colorbar;
-        xlabel('$r$'); ylabel('$z$'); title(sprintf('t=%.3d',Ts(it))); legend('$|\hat\phi|$');
+        xlabel('$r$'); ylabel('$z$'); title(sprintf('t=%.3d',Ts5D(it))); legend('$|\hat\phi|$');
 if strcmp(OUTPUTS.write_non_lin,'.true.')
     subplot(224); plt = @(x) fftshift((abs(x)),2);
         pclr = pcolor((RR),(ZZ),plt(si00(:,:,it))); set(pclr, 'edgecolor','none'); colorbar;
-        xlabel('$r$'); ylabel('$z$'); title(sprintf('t=%.3d',Ts(it))); legend('$|S_i^{00}|$');
+        xlabel('$r$'); ylabel('$z$'); title(sprintf('t=%.3d',Ts5D(it))); legend('$|S_i^{00}|$');
 end
 save_figure
 end
@@ -278,17 +310,17 @@ end
 %%
 if 0
 %% Show frame in kspace
-tf = 20; [~,it] = min(abs(Ts-tf));
+tf = 20; [~,it] = min(abs(Ts5D-tf));
 fig = figure; FIGNAME = ['krkz_frame',sprintf('_%.2d',JOBNUM)];
     subplot(221); plt = @(x) fftshift((abs(x)),2);
         pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(PHI(:,:,it))); set(pclr, 'edgecolor','none'); colorbar;
-        xlabel('$k_r$'); ylabel('$k_z$'); title(sprintf('t=%.3d',Ts(it))); legend('$|\hat\phi|$');
+        xlabel('$k_r$'); ylabel('$k_z$'); title(sprintf('t=%.3d',Ts5D(it))); legend('$|\hat\phi|$');
     subplot(222); plt = @(x) fftshift(abs(x),2);
         pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(Ni00(:,:,it))); set(pclr, 'edgecolor','none'); colorbar;
-        xlabel('$k_r$'); ylabel('$k_z$'); title(sprintf('t=%.3d',Ts(it))); legend('$|\hat n_i^{00}|$');
+        xlabel('$k_r$'); ylabel('$k_z$'); title(sprintf('t=%.3d',Ts5D(it))); legend('$|\hat n_i^{00}|$');
     subplot(223); plt = @(x) fftshift(abs(x),2);
         pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(Ne00(:,:,it))); set(pclr, 'edgecolor','none'); colorbar;
-        xlabel('$k_r$'); ylabel('$k_z$'); title(sprintf('t=%.3d',Ts(it))); legend('$|\hat n_e^{00}|$');
+        xlabel('$k_r$'); ylabel('$k_z$'); title(sprintf('t=%.3d',Ts5D(it))); legend('$|\hat n_e^{00}|$');
 if strcmp(OUTPUTS.write_non_lin,'.true.')
     subplot(224); plt = @(x) fftshift((abs(x)),2);
         pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(Si00(:,:,it))); set(pclr, 'edgecolor','none'); colorbar;
diff --git a/wk/load_results.m b/wk/load_results.m
new file mode 100644
index 00000000..a8c70c17
--- /dev/null
+++ b/wk/load_results.m
@@ -0,0 +1,23 @@
+%% load results %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+JOBNUM = 00;
+FMT = '.fig';
+filename = [BASIC.SIMID,'_','%.2d.h5'];
+filename = sprintf(filename,JOBNUM); disp(['Loading ',filename])
+% Loading from output file
+% CPUTIME   = h5readatt(filename,'/data/input','cpu_time');
+if strcmp(OUTPUTS.write_moments,'.true.') 
+[Nipj, Pi, Ji, kr, kz, Ts5D, dt5D] = load_5D_data(filename, 'moments_i');
+[Nepj, Pe, Je                    ] = load_5D_data(filename, 'moments_e');
+end
+[Ni00, kr, kz, Ts2D, dt2D] = load_2D_data(filename, 'Ni00');
+ Ne00                      = load_2D_data(filename, 'Ne00');
+%Pi = [0]; Ji = Pi; Pe = Pi; Je = Pi;
+% Nipj = zeros(1,1,numel(kr),numel(kz),numel(Ts5D));
+% Nepj = Nipj;
+% Nipj(1,1,:,:,:) = Ni00; Nepj(1,1,:,:,:) = Ne00;
+PHI                            = load_2D_data(filename, 'phi');
+
+if strcmp(OUTPUTS.write_non_lin,'.true.') 
+    Sipj    = load_5D_data(filename, 'Sipj');
+    Sepj    = load_5D_data(filename, 'Sepj');
+end
\ No newline at end of file
-- 
GitLab