From ca8720fdb8da53f19ea048f98c48974d2a7e6eb7 Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Sat, 1 Apr 2023 12:37:06 +0200
Subject: [PATCH] faster and more versatile matlab scripts
---
matlab/load/compile_results_3Da.m | 37 ++++
matlab/load/compile_results_low_mem.m | 95 ++++----
matlab/load/load_3D_data.m | 22 +-
matlab/load/load_3Da_data.m | 43 ++++
matlab/load/load_grid_data.m | 8 +-
matlab/load/load_params.m | 148 +++++++------
.../plot_radial_transport_and_spacetime.m | 75 +------
matlab/plot/show_moments_spectrum.m | 208 ++++--------------
matlab/process_field_v2.m | 29 ++-
9 files changed, 286 insertions(+), 379 deletions(-)
create mode 100644 matlab/load/compile_results_3Da.m
create mode 100644 matlab/load/load_3Da_data.m
diff --git a/matlab/load/compile_results_3Da.m b/matlab/load/compile_results_3Da.m
new file mode 100644
index 00000000..3fffa329
--- /dev/null
+++ b/matlab/load/compile_results_3Da.m
@@ -0,0 +1,37 @@
+function [field, Ts3D] = compile_results_3Da(DIRECTORY,JOBNUMMIN,JOBNUMMAX,fieldname)
+CONTINUE = 1;
+JOBNUM = JOBNUMMIN; JOBFOUND = 0;
+% field
+field = [];
+Ts3D = [];
+ii = 1;
+while(CONTINUE)
+ filename = sprintf([DIRECTORY,'outputs_%.2d.h5'],JOBNUM);
+ % Check presence and jobnummax
+ if (exist(filename, 'file') == 2 && JOBNUM <= JOBNUMMAX)
+ %test if it is corrupted or currently running
+ try
+ openable = ~isempty(h5read(filename,'/data/var3d/time'));
+ catch
+ openable = 0;
+ end
+ if openable
+ % load field %%
+ [ F, T, ~] = load_3Da_data(filename, fieldname);
+ field = cat(5,field,F);
+ Ts3D = cat(1,Ts3D,T);
+ ii = ii + 1;
+ JOBFOUND = JOBFOUND + 1;
+ end
+ elseif (JOBNUM > JOBNUMMAX)
+ CONTINUE = 0;
+ end
+ JOBNUM = JOBNUM + 1;
+end
+
+if(JOBFOUND == 0)
+ disp('no results found, please verify the paths');
+ return;
+end
+
+end
\ No newline at end of file
diff --git a/matlab/load/compile_results_low_mem.m b/matlab/load/compile_results_low_mem.m
index bfe83682..ea3275d8 100644
--- a/matlab/load/compile_results_low_mem.m
+++ b/matlab/load/compile_results_low_mem.m
@@ -1,21 +1,21 @@
function [DATA] = compile_results_low_mem(DATA,DIRECTORY,JOBNUMMIN,JOBNUMMAX)
CONTINUE = 1;
JOBNUM = JOBNUMMIN; JOBFOUND = 0;
-DATA.TJOB_SE = []; % Start and end times of jobs
-DATA.NU_EVOL = []; % evolution of parameter nu between jobs
-DATA.CO_EVOL = []; % evolution of CO
-DATA.MUx_EVOL = []; % evolution of parameter mu between jobs
-DATA.MUy_EVOL = []; % evolution of parameter mu between jobs
-DATA.MUz_EVOL = []; % evolution of parameter mu between jobs
-DATA.K_N_EVOL = []; %
-DATA.K_T_EVOL = []; %
-DATA.L_EVOL = []; %
-DATA.DT_EVOL = []; %
+DATA.params_evol.TJOB_SE = []; % Start and end times of jobs
+DATA.params_evol.NU_EVOL = []; % evolution of parameter nu between jobs
+DATA.params_evol.CO_EVOL = []; % evolution of CO
+DATA.params_evol.MUx_EVOL = []; % evolution of parameter mu between jobs
+DATA.params_evol.MUy_EVOL = []; % evolution of parameter mu between jobs
+DATA.params_evol.MUz_EVOL = []; % evolution of parameter mu between jobs
+DATA.params_evol.K_N_EVOL = []; %
+DATA.params_evol.K_T_EVOL = []; %
+DATA.params_evol.L_EVOL = []; %
+DATA.params_evol.DT_EVOL = []; %
% Low memoery cost data
-HFLUXI_ = [];
+HFLUX_ = [];
HFLUXE_ = [];
-GGAMMAI_ = [];
-PGAMMAI_ = [];
+GGAMMA_ = [];
+PGAMMA_ = [];
Ts0D_ = [];
DATA.outfilenames = {};
ii = 1;
@@ -36,38 +36,38 @@ while(CONTINUE)
% Loading from output file
CPUTIME = h5readatt(filename,'/data/input','cpu_time');
DT_SIM = h5readatt(filename,'/data/input/basic','dt');
- [Pe, Je, Pi, Ji, kx, ky, z] = load_grid_data(filename);
+ [Parray, Jarray, kx, ky, z] = load_grid_data(filename);
W_GAMMA = strcmp(h5readatt(filename,'/data/input/diag_par','write_gamma'),'y');
W_HF = strcmp(h5readatt(filename,'/data/input/diag_par','write_hf' ),'y');
KIN_E = strcmp(h5readatt(filename,'/data/input/model', 'ADIAB_E' ),'n');
BETA = h5readatt(filename,'/data/input/model','beta');
if W_GAMMA
- [ GGAMMA_RI, Ts0D, ~] = load_0D_data(filename, 'gflux_xi');
- PGAMMA_RI = load_0D_data(filename, 'pflux_xi');
- GGAMMAI_ = cat(1,GGAMMAI_,GGAMMA_RI); clear GGAMMA_RI
- PGAMMAI_ = cat(1,PGAMMAI_,PGAMMA_RI); clear PGAMMA_RI
+ [ GGAMMA_R, Ts0D, ~] = load_0D_data(filename, 'gflux_x');
+ PGAMMA_R = load_0D_data(filename, 'pflux_x');
+ GGAMMA_ = cat(1,GGAMMA_,GGAMMA_R); clear GGAMMA_R
+ PGAMMA_ = cat(1,PGAMMA_,PGAMMA_R); clear PGAMMA_R
end
if W_HF
- [ HFLUX_XI, Ts0D, ~] = load_0D_data(filename, 'hflux_xi');
- HFLUXI_ = cat(1,HFLUXI_,HFLUX_XI); clear HFLUX_XI
+ [ HFLUX_X, Ts0D, ~] = load_0D_data(filename, 'hflux_x');
+ HFLUX_ = cat(1,HFLUX_,HFLUX_X); clear HFLUX_X
end
Ts0D_ = cat(1,Ts0D_,Ts0D);
% Evolution of simulation parameters
- load_params
- DATA.TJOB_SE = [DATA.TJOB_SE Ts0D(1) Ts0D(end)];
- DATA.NU_EVOL = [DATA.NU_EVOL DATA.NU DATA.NU];
- DATA.CO_EVOL = [DATA.CO_EVOL DATA.CO DATA.CO];
- DATA.MUx_EVOL = [DATA.MUx_EVOL DATA.MUx DATA.MUx];
- DATA.MUy_EVOL = [DATA.MUy_EVOL DATA.MUy DATA.MUy];
- DATA.MUz_EVOL = [DATA.MUz_EVOL DATA.MUz DATA.MUz];
- DATA.K_N_EVOL = [DATA.K_N_EVOL DATA.K_N DATA.K_N];
- DATA.K_T_EVOL = [DATA.K_T_EVOL DATA.K_T DATA.K_T];
- DATA.L_EVOL = [DATA.L_EVOL DATA.L DATA.L];
- DATA.DT_EVOL = [DATA.DT_EVOL DATA.DT_SIM DATA.DT_SIM];
+ DATA = load_params(DATA,filename);
+ DATA.params_evol.TJOB_SE = [DATA.params_evol.TJOB_SE Ts0D(1) Ts0D(end)];
+ DATA.params_evol.NU_EVOL = [DATA.params_evol.NU_EVOL DATA.inputs.NU DATA.inputs.NU];
+ DATA.params_evol.CO_EVOL = [DATA.params_evol.CO_EVOL DATA.inputs.CO DATA.inputs.CO];
+ DATA.params_evol.MUx_EVOL = [DATA.params_evol.MUx_EVOL DATA.inputs.MUx DATA.inputs.MUx];
+ DATA.params_evol.MUy_EVOL = [DATA.params_evol.MUy_EVOL DATA.inputs.MUy DATA.inputs.MUy];
+ DATA.params_evol.MUz_EVOL = [DATA.params_evol.MUz_EVOL DATA.inputs.MUz DATA.inputs.MUz];
+ DATA.params_evol.K_N_EVOL = [DATA.params_evol.K_N_EVOL DATA.inputs.K_N DATA.inputs.K_N];
+ DATA.params_evol.K_T_EVOL = [DATA.params_evol.K_T_EVOL DATA.inputs.K_T DATA.inputs.K_T];
+ DATA.params_evol.L_EVOL = [DATA.params_evol.L_EVOL DATA.inputs.L DATA.inputs.L];
+ DATA.params_evol.DT_EVOL = [DATA.params_evol.DT_EVOL DATA.inputs.DT_SIM DATA.inputs.DT_SIM];
ii = ii + 1;
JOBFOUND = JOBFOUND + 1;
@@ -118,30 +118,31 @@ else
% scaling
DATA.scale = 1;%(1/Nx/Ny)^2;
% Fields
- DATA.GGAMMA_RI = GGAMMAI_; DATA.PGAMMA_RI = PGAMMAI_; DATA.HFLUX_X = HFLUXI_;
+ DATA.GGAMMA_RI = GGAMMA_; DATA.PGAMMA_RI = PGAMMA_; DATA.HFLUX_X = HFLUX_;
if(KIN_E)
DATA.HFLUX_XE = HFLUXE_;
end
DATA.Ts0D = Ts0D_;
DATA.KIN_E=KIN_E;
% grids
- DATA.Pe = Pe; DATA.Pi = Pi;
- DATA.Je = Je; DATA.Ji = Ji;
- DATA.kx = kx; DATA.ky = ky; DATA.z = z; DATA.Npol = -z(1)/pi;
- DATA.x = x; DATA.y = y;
- DATA.ikx0 = ikx0; DATA.iky0 = iky0;
- DATA.Nx = Nx; DATA.Ny = Ny; DATA.Nz = Nz; DATA.Nkx = Nkx; DATA.Nky = Nky;
- DATA.Pmaxe = numel(Pe); DATA.Pmaxi = numel(Pi); DATA.Jmaxe = numel(Je); DATA.Jmaxi = numel(Ji);
+ DATA.grids.Parray = Parray;
+ DATA.grids.Jarray = Jarray;
+ DATA.grids.kx = kx; DATA.grids.ky = ky; DATA.grids.z = z; DATA.grids.Npol = -z(1)/pi;
+ DATA.grids.x = x; DATA.grids.y = y;
+ DATA.grids.ikx0 = ikx0; DATA.grids.iky0 = iky0;
+ DATA.grids.Nx = Nx; DATA.grids.Ny = Ny; DATA.grids.Nz = Nz; DATA.grids.Nkx = Nkx; DATA.grids.Nky = Nky;
+ DATA.grids.Np = numel(Parray);DATA.grids.Nj = numel(Jarray);
+ DATA.grids.deltap = Parray(2)-Parray(1);
DATA.dir = DIRECTORY;
DATA.localdir = DIRECTORY;
- DATA.param_title=['$\nu_{',DATA.CONAME,'}=$', num2str(DATA.NU), ...
- ', $\kappa_{Ni}=$',num2str(DATA.K_N),', $\kappa_{Ti}=$',num2str(DATA.K_T),...
- ', $L=',num2str(DATA.L),'$, $N=',...
- num2str(DATA.Nx),'$, $(P,J)=(',num2str(DATA.Pi(end)),',',...
- num2str(DATA.Ji(end)),')$,',' $\mu_{hd}=$(',num2str(DATA.MUx),...
- ',',num2str(DATA.MUy),')'];
- DATA.paramshort = [num2str(DATA.Pmaxi),'x',num2str(DATA.Jmaxi),'x',...
- num2str(DATA.Nkx),'x',num2str(DATA.Nky),'x',num2str(DATA.Nz)];
+ DATA.param_title=['$\nu_{',DATA.inputs.CONAME,'}=$', num2str(DATA.inputs.NU), ...
+ ', $\kappa_{Ni}=$',num2str(DATA.inputs.K_N),', $\kappa_{Ti}=$',num2str(DATA.inputs.K_T),...
+ ', $L=',num2str(DATA.inputs.L),'$, $N=',...
+ num2str(DATA.grids.Nx),'$, $(P,J)=(',num2str(DATA.grids.Parray(end)),',',...
+ num2str(DATA.grids.Jarray(end)),')$,',' $\mu_{hd}=$(',num2str(DATA.inputs.MUx),...
+ ',',num2str(DATA.inputs.MUy),')'];
+ DATA.paramshort = [num2str(DATA.grids.Np),'x',num2str(DATA.grids.Nj),'x',...
+ num2str(DATA.grids.Nkx),'x',num2str(DATA.grids.Nky),'x',num2str(DATA.grids.Nz)];
JOBNUM = LASTJOB;
filename = sprintf([DIRECTORY,'outputs_%.2d.h5'],JOBNUM);
diff --git a/matlab/load/load_3D_data.m b/matlab/load/load_3D_data.m
index 47bc4bf0..491b1921 100644
--- a/matlab/load/load_3D_data.m
+++ b/matlab/load/load_3D_data.m
@@ -13,6 +13,10 @@ function [ data, time, dt ] = load_3D_data( filename, variablename )
sz = size(tmp);
cmpx = 0;
end
+ % add a z dimension even if 2D
+ if(numel(sz) == 2)
+ sz = [sz 1];
+ end
% add time dimension
sz = [sz numel(time)];
data = zeros(sz);
@@ -20,16 +24,22 @@ function [ data, time, dt ] = load_3D_data( filename, variablename )
for it = 1:numel(time)
tmp = h5read(filename,['/data/var3d/',variablename,'/', num2str(cstart+it,'%06d')]);
if cmpx
- if(numel(sz) == 3)
- data(:,:,it) = tmp.real + 1i * tmp.imaginary;
- else
+ switch numel(sz)
+ case(3)
+ data(:,:,1,it) = tmp.real + 1i * tmp.imaginary;
+ case(4)
data(:,:,:,it) = tmp.real + 1i * tmp.imaginary;
+ case(5)
+ data(:,:,:,:,it) = tmp.real + 1i * tmp.imaginary;
end
else
- if(numel(sz) == 3)
- data(:,:,it) = tmp;
- else
+ switch numel(sz)
+ case(3)
+ data(:,:,1,it) = tmp;
+ case(4)
data(:,:,:,it) = tmp;
+ case(5)
+ data(:,:,:,:,it) = tmp;
end
end
end
diff --git a/matlab/load/load_3Da_data.m b/matlab/load/load_3Da_data.m
new file mode 100644
index 00000000..d3dc8bec
--- /dev/null
+++ b/matlab/load/load_3Da_data.m
@@ -0,0 +1,43 @@
+function [ data, time, dt ] = load_3Da_data( filename, variablename )
+%LOAD_3D_DATA load a 3D variable stored in a hdf5 result file from HeLaZ
+ time = h5read(filename,'/data/var3d/time');
+ dt = h5readatt(filename,'/data/input/basic','dt');
+ cstart= h5readatt(filename,'/data/input/basic','start_iframe3d');
+ Na = h5readatt(filename,'/data/input/model','Na');
+ Np = h5readatt(filename,'/data/input/grid', 'Np');
+ Nj = h5readatt(filename,'/data/input/grid', 'Nj');
+ Nky = h5readatt(filename,'/data/input/grid', 'Nky');
+ Nkx = h5readatt(filename,'/data/input/grid', 'Nkx');
+ Nz = h5readatt(filename,'/data/input/grid', 'Nz');
+
+
+ % Find array size by loading the first output
+ tmp = h5read(filename,['/data/var3d/',variablename,'/', num2str(cstart+1,'%06d')]);
+ try % check if it is complex or real
+ sz = size(tmp.real);
+ cmpx = 1;
+ catch
+ sz = size(tmp);
+ cmpx = 0;
+ end
+ % add a dimension if nz=1 or na=1
+% if Na == 1
+% sz = [1 sz];
+% end
+ if Nz == 1
+ sz = [sz 1];
+ end
+ % add time dimension
+ data = zeros([sz numel(time)]);
+
+ for it = 1:numel(time)
+ tmp = h5read(filename,['/data/var3d/',variablename,'/', num2str(cstart+it,'%06d')]);
+ if cmpx
+ data(:,:,:,:,it) = reshape(tmp.real,sz) + 1i * reshape(tmp.imaginary,sz);
+ else
+ data(:,:,:,:,it) = reshape(tmp,sz);
+ end
+ end
+
+end
+
diff --git a/matlab/load/load_grid_data.m b/matlab/load/load_grid_data.m
index 2a1dd8c9..6df1d511 100644
--- a/matlab/load/load_grid_data.m
+++ b/matlab/load/load_grid_data.m
@@ -1,9 +1,7 @@
-function [ pe, je, pi, ji, kx, ky, z ] = load_grid_data( filename )
+function [ p, j, kx, ky, z ] = load_grid_data( filename )
%LOAD_GRID_DATA stored in a hdf5 result file from HeLaZ
- pe = h5read(filename,'/data/grid/coordp');
- je = h5read(filename,'/data/grid/coordj');
- pi = h5read(filename,'/data/grid/coordp');
- ji = h5read(filename,'/data/grid/coordj');
+ p = h5read(filename,'/data/grid/coordp');
+ j = h5read(filename,'/data/grid/coordj');
kx = h5read(filename,'/data/grid/coordkx');
ky = h5read(filename,'/data/grid/coordky');
z = h5read(filename,'/data/grid/coordz');
diff --git a/matlab/load/load_params.m b/matlab/load/load_params.m
index 56e12861..73f72a92 100644
--- a/matlab/load/load_params.m
+++ b/matlab/load/load_params.m
@@ -1,82 +1,84 @@
-DATA.CO = h5readatt(filename,'/data/input/coll','CO');
-DATA.K_N = h5readatt(filename,'/data/input/ions','k_N');
-DATA.K_T = h5readatt(filename,'/data/input/ions','k_T');
-DATA.Q0 = h5readatt(filename,'/data/input/geometry','q0');
-DATA.EPS = h5readatt(filename,'/data/input/geometry','eps');
-DATA.SHEAR = h5readatt(filename,'/data/input/geometry','shear');
-DATA.GEOM = h5readatt(filename,'/data/input/geometry','geometry');
-% DATA.KAPPA = h5readatt(filename,'/data/input/geometry','kappa');
-% DATA.DELTA = h5readatt(filename,'/data/input/geometry','delta');
+function [DATA] = load_params(DATA,filename)
+ DATA.inputs.CO = h5readatt(filename,'/data/input/coll','CO');
+ DATA.inputs.K_N = h5readatt(filename,'/data/input/ions','k_N');
+ DATA.inputs.K_T = h5readatt(filename,'/data/input/ions','k_T');
+ DATA.inputs.Q0 = h5readatt(filename,'/data/input/geometry','q0');
+ DATA.inputs.EPS = h5readatt(filename,'/data/input/geometry','eps');
+ DATA.inputs.SHEAR = h5readatt(filename,'/data/input/geometry','shear');
+ DATA.inputs.GEOM = h5readatt(filename,'/data/input/geometry','geometry');
+ % DATA.KAPPA = h5readatt(filename,'/data/input/geometry','kappa');
+ % DATA.DELTA = h5readatt(filename,'/data/input/geometry','delta');
-DATA.DT_SIM = h5readatt(filename,'/data/input/basic','dt');
-DATA.PMAX = h5readatt(filename,'/data/input/grid','pmax');
-DATA.JMAX = h5readatt(filename,'/data/input/grid','jmax');
-DATA.Nx = h5readatt(filename,'/data/input/grid','Nx');
-DATA.Ny = h5readatt(filename,'/data/input/grid','Ny');
-DATA.L = h5readatt(filename,'/data/input/grid','Lx');
-try
-DATA.CLOS = h5readatt(filename,'/data/input/model','CLOS');
-DATA.NL_CLOS = h5readatt(filename,'/data/input/model','NL_CLOS');
-catch
+ DATA.inputs.DT_SIM = h5readatt(filename,'/data/input/basic','dt');
+ DATA.inputs.PMAX = h5readatt(filename,'/data/input/grid','pmax');
+ DATA.inputs.JMAX = h5readatt(filename,'/data/input/grid','jmax');
+ DATA.inputs.Nx = h5readatt(filename,'/data/input/grid','Nx');
+ DATA.inputs.Ny = h5readatt(filename,'/data/input/grid','Ny');
+ DATA.inputs.L = h5readatt(filename,'/data/input/grid','Lx');
try
- DATA.ha_cl = h5readatt(filename,'/data/input/closure','hierarchy_closure');
- DATA.CLOS = h5readatt(filename,'/data/input/closure','dmax');
- DATA.nl_cl = h5readatt(filename,'/data/input/closure','nonlinear_closure');
- DATA.NL_CLOS = h5readatt(filename,'/data/input/closure','nmax');
+ DATA.inputs.CLOS = h5readatt(filename,'/data/input/model','CLOS');
+ DATA.inputs.NL_CLOS = h5readatt(filename,'/data/input/model','NL_CLOS');
catch
- DATA.CLOS = 99;
- DATA.NL_CLOS = 99;
+ try
+ DATA.inputs.ha_cl = h5readatt(filename,'/data/input/closure','hierarchy_closure');
+ DATA.inputs.CLOS = h5readatt(filename,'/data/input/closure','dmax');
+ DATA.inputs.nl_cl = h5readatt(filename,'/data/input/closure','nonlinear_closure');
+ DATA.inputs.NL_CLOS = h5readatt(filename,'/data/input/closure','nmax');
+ catch
+ DATA.inputs.CLOS = 99;
+ DATA.inputs.NL_CLOS = 99;
+ end
end
-end
-DATA.Na = h5readatt(filename,'/data/input/model','Na');
-DATA.NU = h5readatt(filename,'/data/input/model','nu');
-DATA.MUp = h5readatt(filename,'/data/input/model','mu_p');
-DATA.MUj = h5readatt(filename,'/data/input/model','mu_j');
-DATA.MUx = h5readatt(filename,'/data/input/model','mu_x');
-DATA.MUy = h5readatt(filename,'/data/input/model','mu_y');
-DATA.MUz = h5readatt(filename,'/data/input/model','mu_z');
-DATA.LINEARITY = h5readatt(filename,'/data/input/model','LINEARITY');
-DATA.BETA = h5readatt(filename,'/data/input/model','beta');
-DATA.TAU_E = h5readatt(filename,'/data/input/model','tau_e');
-DATA.HYP_V = h5readatt(filename,'/data/input/model','HYP_V');
-DATA.K_cB = h5readatt(filename,'/data/input/model','k_cB');
-DATA.K_gB = h5readatt(filename,'/data/input/model','k_gB');
+ DATA.inputs.Na = h5readatt(filename,'/data/input/model','Na');
+ DATA.inputs.NU = h5readatt(filename,'/data/input/model','nu');
+ DATA.inputs.MUp = h5readatt(filename,'/data/input/model','mu_p');
+ DATA.inputs.MUj = h5readatt(filename,'/data/input/model','mu_j');
+ DATA.inputs.MUx = h5readatt(filename,'/data/input/model','mu_x');
+ DATA.inputs.MUy = h5readatt(filename,'/data/input/model','mu_y');
+ DATA.inputs.MUz = h5readatt(filename,'/data/input/model','mu_z');
+ DATA.inputs.LINEARITY = h5readatt(filename,'/data/input/model','LINEARITY');
+ DATA.inputs.BETA = h5readatt(filename,'/data/input/model','beta');
+ DATA.inputs.TAU_E = h5readatt(filename,'/data/input/model','tau_e');
+ DATA.inputs.HYP_V = h5readatt(filename,'/data/input/model','HYP_V');
+ DATA.inputs.K_cB = h5readatt(filename,'/data/input/model','k_cB');
+ DATA.inputs.K_gB = h5readatt(filename,'/data/input/model','k_gB');
-DATA.W_GAMMA = h5readatt(filename,'/data/input/diag_par','write_gamma') == 'y';
-DATA.W_PHI = h5readatt(filename,'/data/input/diag_par','write_phi') == 'y';
-DATA.W_NA00 = h5readatt(filename,'/data/input/diag_par','write_Na00') == 'y';
-DATA.W_NAPJ = h5readatt(filename,'/data/input/diag_par','write_Napj') == 'y';
-DATA.W_SAPJ = h5readatt(filename,'/data/input/diag_par','write_Sapj') == 'y';
+ DATA.inputs.W_GAMMA = h5readatt(filename,'/data/input/diag_par','write_gamma') == 'y';
+ DATA.inputs.W_PHI = h5readatt(filename,'/data/input/diag_par','write_phi') == 'y';
+ DATA.inputs.W_NA00 = h5readatt(filename,'/data/input/diag_par','write_Na00') == 'y';
+ DATA.inputs.W_NAPJ = h5readatt(filename,'/data/input/diag_par','write_Napj') == 'y';
+ DATA.inputs.W_SAPJ = h5readatt(filename,'/data/input/diag_par','write_Sapj') == 'y';
-% Species dependent parameters
-DATA.sigma = zeros(1,DATA.Na);
-DATA.tau = zeros(1,DATA.Na);
-DATA.q = zeros(1,DATA.Na);
-DATA.K_N = zeros(1,DATA.Na);
-DATA.K_T = zeros(1,DATA.Na);
-spnames = {'ions','electrons'};
-for ia=1:DATA.Na
- spdata = ['/data/input/',spnames{ia}];
- DATA.sigma(ia) = h5readatt(filename,spdata,'sigma');
- DATA.tau(ia) = h5readatt(filename,spdata,'tau');
- DATA.q(ia) = h5readatt(filename,spdata,'q');
- DATA.K_N(ia) = h5readatt(filename,spdata,'k_N');
- DATA.K_T(ia) = h5readatt(filename,spdata,'k_T');
-end
-DATA.spnames = spnames{1:DATA.Na};
-DATA.CONAME = DATA.CO;
+ % Species dependent parameters
+ DATA.inputs.sigma = zeros(1,DATA.inputs.Na);
+ DATA.inputs.tau = zeros(1,DATA.inputs.Na);
+ DATA.inputs.q = zeros(1,DATA.inputs.Na);
+ DATA.inputs.K_N = zeros(1,DATA.inputs.Na);
+ DATA.inputs.K_T = zeros(1,DATA.inputs.Na);
+ spnames = {'ions','electrons'};
+ for ia=1:DATA.inputs.Na
+ spdata = ['/data/input/',spnames{ia}];
+ DATA.inputs.sigma(ia) = h5readatt(filename,spdata,'sigma');
+ DATA.inputs.tau(ia) = h5readatt(filename,spdata,'tau');
+ DATA.inputs.q(ia) = h5readatt(filename,spdata,'q');
+ DATA.inputs.K_N(ia) = h5readatt(filename,spdata,'k_N');
+ DATA.inputs.K_T(ia) = h5readatt(filename,spdata,'k_T');
+ end
+ DATA.inputs.spnames = spnames{1:DATA.inputs.Na};
+ DATA.inputs.CONAME = DATA.inputs.CO;
-if (DATA.CLOS == 0); DATA.CLOSNAME = 'Trunc.';
-elseif(DATA.CLOS == 1); DATA.CLOSNAME = 'Clos. 1';
-elseif(DATA.CLOS == 2); DATA.CLOSNAME = 'Clos. 2';
-end
+ if (DATA.inputs.CLOS == 0); DATA.CLOSNAME = 'Trunc.';
+ elseif(DATA.inputs.CLOS == 1); DATA.CLOSNAME = 'Clos. 1';
+ elseif(DATA.inputs.CLOS == 2); DATA.CLOSNAME = 'Clos. 2';
+ end
-degngrad = ['P_',num2str(DATA.PMAX),'_J_',num2str(DATA.JMAX)];
+ degngrad = ['P_',num2str(DATA.inputs.PMAX),'_J_',num2str(DATA.inputs.JMAX)];
-degngrad = [degngrad,'_Kni_%1.1f_nu_%0.0e_',...
- DATA.CONAME,'_CLOS_',num2str(DATA.CLOS),'_mu_%0.0e'];
-degngrad = sprintf(degngrad,[DATA.K_N,DATA.NU,DATA.MUx]);
-% if ~DATA.LINEARITY; degngrad = ['lin_',degngrad]; end
-resolution = [num2str(DATA.Nx),'x',num2str(DATA.Ny),'_'];
-gridname = ['L_',num2str(DATA.L),'_'];
-DATA.PARAMS = [resolution,gridname,degngrad];
\ No newline at end of file
+ degngrad = [degngrad,'_Kni_%1.1f_nu_%0.0e_',...
+ DATA.inputs.CONAME,'_CLOS_',num2str(DATA.inputs.CLOS),'_mu_%0.0e'];
+ degngrad = sprintf(degngrad,[DATA.inputs.K_N,DATA.inputs.NU,DATA.inputs.MUx]);
+ % if ~DATA.LINEARITY; degngrad = ['lin_',degngrad]; end
+ resolution = [num2str(DATA.inputs.Nx),'x',num2str(DATA.inputs.Ny),'_'];
+ gridname = ['L_',num2str(DATA.inputs.L),'_'];
+ DATA.params_string = [resolution,gridname,degngrad];
+end
\ No newline at end of file
diff --git a/matlab/plot/plot_radial_transport_and_spacetime.m b/matlab/plot/plot_radial_transport_and_spacetime.m
index 3ee858a7..0234d379 100644
--- a/matlab/plot/plot_radial_transport_and_spacetime.m
+++ b/matlab/plot/plot_radial_transport_and_spacetime.m
@@ -10,12 +10,6 @@ function [FIGURE] = plot_radial_transport_and_spacetime(DATA, OPTIONS)
Gx_infty_std = std (DATA.PGAMMA_RI(its0D:ite0D))*SCALE;
Qx_infty_avg = mean(DATA.HFLUX_X(its0D:ite0D))*SCALE;
Qx_infty_std = std (DATA.HFLUX_X(its0D:ite0D))*SCALE;
-% disp(['Q_x=',sprintf('%2.2e',Qx_infty_avg),'+-',sprintf('%2.2e',Qx_infty_std)]);
-% f_avg_z = squeeze(mean(DATA.PHI(:,:,:,:),3));
-% [~,ikzf] = max(squeeze(mean(abs(f_avg_z(1,:,its3D:ite3D)),3)));
-% ikzf = min([ikzf,DATA.Nky]);
-% Ns3D = numel(DATA.Ts3D);
-% [KX, KY] = meshgrid(DATA.kx, DATA.ky);
%% error estimation
DT_ = (tend-tstart)/OPTIONS.NCUT;
Qx_ee = zeros(1,OPTIONS.NCUT);
@@ -27,41 +21,16 @@ function [FIGURE] = plot_radial_transport_and_spacetime(DATA, OPTIONS)
Qx_avg = mean(Qx_ee);
Qx_err = std(Qx_ee);
disp(['Q_avg=',sprintf('%2.2e',Qx_avg),'+-',sprintf('%2.2e',Qx_err)]);
-% %% computations
-%
-% % Compute zonal and non zonal energies
-% E_Zmode_SK = zeros(1,Ns3D);
-% E_NZmode_SK = zeros(1,Ns3D);
-% for it = 1:numel(DATA.Ts3D)
-% E_Zmode_SK(it) = squeeze(DATA.ky(ikzf).^2.*abs(squeeze(f_avg_z(ikzf,1,it))).^2);
-% E_NZmode_SK(it) = squeeze(sum(sum(((1+KX.^2+KY.^2).*abs(squeeze(f_avg_z(:,:,it))).^2.*(KY~=0)))));
-% end
-% % Compute thermodynamic entropy Eq.(5) Navarro et al. 2012 PoP
-% % 1/2 sum_p sum_j Napj^2(k=0) (avg z)
-% Nipjz = sum(sum(sum(sum(conj(DATA.Nipj).*DATA.Nipj))));
-% ff = trapz(DATA.z,Nipjz,5);
-% E_TE = 0.5*squeeze(ff);
-% % Compute electrostatic energy
-% E_ES = zeros(size(DATA.Ts5D));
-% bi = sqrt(KX.^2+KY.^2)*DATA.sigma(1)*sqrt(2*DATA.tau(1)); %argument of the kernel
-% for it5D = 1:numel(DATA.Ts5D)
-% [~,it3D] = min(abs(DATA.Ts3D-DATA.Ts5D(it5D)));
-% for in = 1:DATA.Jmaxi
-% Knphi = kernel(in-1,bi).*squeeze(trapz(DATA.z,DATA.PHI(:,:,:,it3D),3));
-% Ni0n_z= squeeze(trapz(DATA.z,DATA.Nipj(1,in,:,:,:,it5D),5));
-% E_ES(it5D) = 0.5*sum(sum(abs(conj(Knphi).*Ni0n_z)));
-% end
-% end
%% Figure
clr_ = lines(20);
mvm = @(x) movmean(x,OPTIONS.NMVA);
- FIGURE.fig = figure; FIGURE.FIGNAME = ['ZF_transport_drphi','_',DATA.PARAMS]; %set(gcf, 'Position', [500, 1000, 1000, 600])
+ FIGURE.fig = figure; FIGURE.FIGNAME = ['ZF_transport_drphi','_',DATA.params_string]; %set(gcf, 'Position', [500, 1000, 1000, 600])
FIGURE.ax1 = subplot(2,1,1,'parent',FIGURE.fig);
plot(mvm(DATA.Ts0D),mvm(DATA.PGAMMA_RI*SCALE),'--',...
- 'color',clr_((DATA.Pmaxi-1)/2,:),...
+ 'color',clr_((DATA.grids.Np-1)/2,:),...
'DisplayName',['$\Gamma_x$ ',DATA.paramshort]); hold on;
plot(mvm(DATA.Ts0D),mvm(DATA.HFLUX_X*SCALE),'-',...
- 'color',clr_((DATA.Pmaxi-1)/2,:),...
+ 'color',clr_((DATA.grids.Np-1)/2,:),...
'DisplayName',['$Q_x$ ',DATA.paramshort]); hold on;
ylabel('Transport')
if(~isnan(Qx_infty_avg))
@@ -79,24 +48,13 @@ mvm = @(x) movmean(x,OPTIONS.NMVA);
grid on; set(gca,'xticklabel',[]);
title({DATA.param_title,...
['$\Gamma^{\infty} = $',num2str(Gx_infty_avg),'$, Q^{\infty} = $',num2str(Qx_infty_avg)]});
-
-% %% Free energy
-% FIGURE.ax2 = subplot(3,1,2,'parent',FIGURE.fig);
-% yyaxis left
-% plot(DATA.Ts5D,E_TE,'DisplayName','$\epsilon_f$'); hold on;
-% ylabel('Entropy');%('$\epsilon_f$')
-% yyaxis right
-% plot(DATA.Ts5D,E_ES,'DisplayName','$\epsilon_\phi$');
-% ylabel('ES energy');%('$\epsilon_\phi$')
-% xlim([DATA.Ts5D(1), DATA.Ts5D(end)]);
-% xlabel('$t c_s/R$'); grid on; set(gca,'xticklabel',[]);% xlim([0 500]);
%% radial shear radial profile
% computation
Ns3D = numel(DATA.Ts3D);
- [KX, KY] = meshgrid(DATA.kx, DATA.ky);
+ [KX, KY] = meshgrid(DATA.grids.kx, DATA.grids.ky);
plt = @(x) mean(x(:,:,:),1);
- kycut = max(DATA.ky);
- kxcut = max(DATA.kx);
+ kycut = max(DATA.grids.ky);
+ kxcut = max(DATA.grids.kx);
LP = (abs(KY)<kycut).*(abs(KX)<kxcut); %Low pass filter
OPTIONS.NAME = OPTIONS.ST_FIELD;
@@ -104,12 +62,12 @@ mvm = @(x) movmean(x,OPTIONS.NMVA);
OPTIONS.COMP = 'avg';
OPTIONS.TIME = DATA.Ts3D;
OPTIONS.POLARPLOT = 0;
- toplot = process_field(DATA,OPTIONS);
+ toplot = process_field_v2(DATA,OPTIONS);
f2plot = toplot.FIELD;
dframe = ite3D - its3D;
clim = max(max(max(abs(plt(f2plot(:,:,:))))));
FIGURE.ax3 = subplot(2,1,2,'parent',FIGURE.fig);
- [TY,TX] = meshgrid(DATA.x,DATA.Ts3D(toplot.FRAMES));
+ [TY,TX] = meshgrid(DATA.grids.x,DATA.Ts3D(toplot.FRAMES));
pclr = pcolor(TX,TY,squeeze(plt(f2plot))');
set(pclr, 'edgecolor','none');
legend(['$\langle ',OPTIONS.ST_FIELD,' \rangle_{y,z}$']) %colorbar;
@@ -124,21 +82,4 @@ mvm = @(x) movmean(x,OPTIONS.NMVA);
subplot(311)
plot(DATA.Ts3D,squeeze(mean(plt(f2plot),1)));
end
- %% Zonal vs NZonal energies
-% subplot(312)
-% it0 = 1; itend = Ns3D;
-% trange = toplot.FRAMES;
-% plt1 = @(x) x;%-x(1);
-% plt2 = @(x) x./max((x(:)));
-% toplot = sum(squeeze(plt(f2plot)).^2,1); %ST from before
-% % plty = @(x) x(500:end)./max(squeeze(x(500:end)));
-% yyaxis left
-% % plot(plt1(DATA.Ts3D(trange)),plt2(E_Zmode_SK(trange)),'DisplayName','$k_{zf}^2|\phi_{kzf}|^2$');
-% plot(plt1(DATA.Ts3D(trange)),plt2(toplot(:)),'DisplayName','Sum $A^2$');
-% ylim([-0.1, 1.5]); ylabel('$E_{Z}$')
-% yyaxis right
-% plot(plt1(DATA.Ts3D(trange)),plt2(E_NZmode_SK(trange)),'DisplayName','$(1+k^2)|\phi^k|^2$');
-% xlim([DATA.Ts3D(it0), DATA.Ts3D(itend)]);
-% ylim([-0.1, 1.5]); ylabel('$E_{NZ}$')
-% xlabel('$t c_s/R$'); grid on; set(gca,'xticklabel',[]);% xlim([0 500]);
end
\ No newline at end of file
diff --git a/matlab/plot/show_moments_spectrum.m b/matlab/plot/show_moments_spectrum.m
index d3839297..1f93b35d 100644
--- a/matlab/plot/show_moments_spectrum.m
+++ b/matlab/plot/show_moments_spectrum.m
@@ -1,145 +1,42 @@
function [ FIGURE ] = show_moments_spectrum( DATA, OPTIONS )
+species_name = {'i','e'}; % hard coded because this list is not output yet
-Pi = DATA.Pi;
-Ji = DATA.Ji;
-if ~(isempty(DATA.Nipjz))
- Time_ = DATA.Ts3D;
- Nipj = sum(abs(DATA.Nipjz),3);
-else
- Time_ = DATA.Ts5D;
-% Nipj = sum(sum(sum(abs(DATA.Nipj).^2,3),4),5);
- Nipj = sum(sum(sum(conj(DATA.Nipj).*DATA.Nipj,3),4),5);
-end
-Nipj = squeeze(Nipj);
-
-if DATA.KIN_E
-Pe = DATA.Pe;
-Je = DATA.Je;
- if ~(isempty(DATA.Nepjz))
- Nepj = sum(abs(DATA.Nepjz),3);
- else
- Nepj = sum(sum(sum(abs(DATA.Nepj),3),4),5);
- end
- Nepj = squeeze(Nepj);
-end
-
-FIGURE.fig = figure; FIGURE.FIGNAME = ['mom_spectrum_',DATA.PARAMS];
+Pa = DATA.grids.Parray;
+Ja = DATA.grids.Jarray;
+Time_ = DATA.Ts3D;
+FIGURE.fig = figure; FIGURE.FIGNAME = ['mom_spectrum_',DATA.params_string];
set(gcf, 'Position', [100 50 1000 400])
-if ~OPTIONS.ST
- t0 = OPTIONS.TIME(1);
- t1 = OPTIONS.TIME(end);
- [~,it0] = min(abs(t0-Time_));
- [~,it1] = min(abs(t1-Time_));
- Nipj = mean(Nipj(:,:,it0:it1),3);
- if DATA.KIN_E
- Nepj = mean(Nepj(:,:,it0:it1),3);
- end
- if numel(OPTIONS.TIME) == 1
- TITLE=['$t=',num2str(OPTIONS.TIME),'$'];
- else
- TITLE=['$t\in[',num2str(t0),',',num2str(t1),']$'];
- end
- Nipj = squeeze(Nipj);
-
- ymini = min(Nipj); ymaxi = max(Nipj);
-
- if DATA.KIN_E
- Nepj = squeeze(Nepj);
- ymine = min(Nepj); ymaxe = max(Nepj);
- ymax = max([ymaxi ymaxe]);
- ymin = min([ymini ymine]);
- end
- if DATA.KIN_E
- subplot(121)
- end
- if ~OPTIONS.P2J
- for ij = 1:numel(Ji)
- name = ['$j=',num2str(Ji(ij)),'$'];
- semilogy(Pi,Nipj(:,ij),'o-','DisplayName',name); hold on;
- end
- xlabel('$p$');
- else
- for ij = 1:numel(Ji)
- name = ['$j=',num2str(Ji(ij)),'$'];
- semilogy(Pi+2*Ji(ij),Nipj(:,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.KIN_E
- subplot(122)
- if ~OPTIONS.P2J
- for ij = 1:numel(Je)
- name = ['$j=',num2str(Je(ij)),'$'];
- semilogy(Pe,Nepj(:,ij),'o-','DisplayName',name); hold on;
- end
- xlabel('p');
- else
- for ij = 1:numel(Je)
- name = ['$j=',num2str(Je(ij)),'$'];
- semilogy(Pe+2*Je(ij),Nepj(:,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
-else
+for ia = 1:DATA.inputs.Na
+ Napjz = sum(abs(squeeze(DATA.Napjz(ia,:,:,:,:))),3);
+ subplot(double(DATA.inputs.Na),1,double(ia))
if OPTIONS.P2J
- plotname = '$\langle\sum_k |N_i^{pj}|\rangle_{p+2j=const}$';
- [JJ,PP] = meshgrid(Ji,Pi);
- P2Ji = PP + 2*JJ;
+ plotname = ['$\langle\sum_k |N_',species_name{ia},'^{pj}|\rangle_{p+2j=const}$'];
+ [JJ,PP] = meshgrid(Ja,Pa);
+ P2J = PP + 2*JJ;
% form an axis of velocity ordered moments
- y_ = unique(P2Ji); yname = '$p+2j$';
+ p2j = unique(P2J); yname = '$p+2j$';
% weights to average
- weights = zeros(size(y_));
+ weights = zeros(size(p2j));
% space time of moments amplitude wrt to p+2j degree
- Ni_ST = zeros(numel(y_),numel(Time_));
+ Na_ST = zeros(numel(p2j),numel(Time_));
% fill the st diagramm by averaging every p+2j=n moments together
- for ip = 1:numel(Pi)
- for ij = 1:numel(Ji)
- [~,ip2j] = min(abs(y_-(Pi(ip)+2*Ji(ij))));
- Ni_ST(ip2j,:) = Ni_ST(ip2j,:) + transpose(squeeze(Nipj(ip,ij,:)));
+ for ip = 1:numel(Pa)
+ for ij = 1:numel(Ja)
+ [~,ip2j] = min(abs(p2j-(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(y_)
- Ni_ST(ip2j,:) = Ni_ST(ip2j,:)/weights(ip2j);
- end
- if DATA.KIN_E
- % same for electrons!!
- [JJ,PP] = meshgrid(Je,Pe);
- P2Je = PP + 2*JJ;
- % form an axis of velocity ordered moments
- p2je = unique(P2Je);
- % weights to average
- weights = zeros(size(p2je));
- % space time of moments amplitude wrt to p+2j degree
- Ne_ST = zeros(numel(p2je),numel(Time_));
- % fill the st diagramm by averaging every p+2j=n moments together
- for ip = 1:numel(Pe)
- for ij = 1:numel(Je)
- [~,ip2j] = min(abs(y_-(Pe(ip)+2*Je(ij))));
- Ne_ST(ip2j,:) = Ne_ST(ip2j,:) + transpose(squeeze(Nepj(ip,ij,:)));
- weights(ip2j) = weights(ip2j) + 1;
- end
- end
- % doing the average
- for ip2j = 1:numel(y_)
- Ne_ST(ip2j,:) = Ne_ST(ip2j,:)/weights(ip2j);
- end
+ for ip2j = 1:numel(p2j)
+ Na_ST(ip2j,:) = Na_ST(ip2j,:)/weights(ip2j);
end
ticks_labels = {};
else % We just order our moments w.r.t. to the convention ..
% (0,0) (1,0) (2,0) (0,1) (3,0) (1,1) (4,0) (2,1) (0,2) etc.
plotname = '$\langle\sum_k |N_i^{pj}|^2\rangle_z$';
- Nmoments = numel(Nipj(:,:,1)); % total number of moments
+ Nmoments = numel(Napjz(:,:,1)); % total number of moments
HL_deg = zeros(Nmoments,2); % list of degrees, first column Hermite, second Laguerre
im = 2;
deg = 1; % the zero degree is always here first place
@@ -149,42 +46,34 @@ else
FOUND = 1;
while(FOUND) % As we find a pair matching the degree we retry
FOUND = 0;
- for ij = 1:DATA.Jmaxi
- for ip = 1:DATA.Pmaxi
- if((ip-1)+2*(ij-1) == deg)
- % Check if the pair is already added
- check_ = HL_deg == [DATA.Pi(ip) DATA.Pi(ij)];
- check_ = sum(check_(:,1) .* check_(:,2));
- if ~check_ % if not add it
- HL_deg(im,1) = DATA.Pi(ip);
- HL_deg(im,2) = DATA.Pi(ij);
- ticks_labels{im} = ['(',num2str(DATA.Pi(ip)),',',num2str(DATA.Ji(ij)),')'];
- im = im + 1; FOUND = 1;
+ for ij = 1:DATA.grids.Nj
+ for ip = 1:DATA.grids.Np
+ if((ip-1)+2*(ij-1) == deg)
+ % Check if the pair is already added
+ check_ = HL_deg == [DATA.grids.Parray(ip) DATA.grids.Jarray(ij)];
+ check_ = sum(check_(:,1) .* check_(:,2));
+ if ~check_ % if not add it
+ HL_deg(im,1) = DATA.grids.Parray(ip);
+ HL_deg(im,2) = DATA.grids.Jarray(ij);
+ ticks_labels{im} = ['(',num2str(DATA.grids.Parray(ip)),',',num2str(DATA.grids.Jarray(ij)),')'];
+ im = im + 1; FOUND = 1;
+ end
+ end
end
- end
end
end
- end
% No pair found anymore, increase the degree
deg = deg + 1;
end
-
+
% form an axis of velocity ordered moments
- y_ = 1:Nmoments; yname = '$(P,J)$';
- % space time of moments amplitude wrt to p+2j degree
- Ni_ST = zeros(Nmoments,numel(Time_));
- for i_ = 1:Nmoments
- Ni_ST(i_,:) = Nipj(HL_deg(i_,1)+1,HL_deg(i_,2)+1,:);
- end
- if DATA.KIN_E
+ p2j = 1:Nmoments; yname = '$(P,J)$';
% space time of moments amplitude wrt to p+2j degree
- Ne_ST = zeros(Nmoments,numel(Time_));
+ Na_ST = zeros(Nmoments,numel(Time_));
for i_ = 1:Nmoments
- Ne_ST = Nepj(HL_deg(i_,1)+1,HL_deg(i_,2)+1,:);
+ Na_ST(i_,:) = Napjz((HL_deg(i_,1)/DATA.grids.deltap)+1,HL_deg(i_,2)+1,:);
end
- end
-
- end
+ end
% plots
% scaling
if OPTIONS.NORMALIZED
@@ -192,29 +81,16 @@ else
else
plt = @(x,i) x;
end
- if DATA.KIN_E
- subplot(2,1,1)
- end
-
- imagesc(Time_,y_,plt(Ni_ST,1:numel(y_)));
+
+ imagesc(Time_,p2j,plt(Na_ST,1:numel(p2j)));
set(gca,'YDir','normal')
xlabel('$t$'); ylabel(yname);
if ~isempty(ticks_labels)
- yticks(y_);
+ yticks(p2j);
yticklabels(ticks_labels)
end
title(plotname)
-
- if DATA.KIN_E
- subplot(2,1,2)
- imagesc(Time_,p2je,plt(Ne_ST,1:numel(y_)));
- set(gca,'YDir','normal')
- xlabel('$t$'); ylabel(yname)
- title(plotname)
- suptitle(DATA.param_title);
- end
end
-
end
diff --git a/matlab/process_field_v2.m b/matlab/process_field_v2.m
index f744d7e4..5b2d103e 100644
--- a/matlab/process_field_v2.m
+++ b/matlab/process_field_v2.m
@@ -9,35 +9,35 @@ for i = 1:numel(OPTIONS.TIME)
end
FRAMES = unique(FRAMES);
%% Setup the plot geometry
-[KX, KY] = meshgrid(DATA.kx, DATA.ky);
+[KX, KY] = meshgrid(DATA.grids.kx, DATA.grids.ky);
directions = {'y','x','z'};
-Nx = DATA.Nx; Ny = DATA.Ny; Nz = DATA.Nz; Nt = numel(FRAMES);
+Nx = DATA.grids.Nx; Ny = DATA.grids.Ny; Nz = DATA.grids.Nz; Nt = numel(FRAMES);
POLARPLOT = OPTIONS.POLARPLOT;
LTXNAME = OPTIONS.NAME;
switch OPTIONS.PLAN
case 'xy'
XNAME = '$x$'; YNAME = '$y$';
- [X,Y] = meshgrid(DATA.x,DATA.y);
+ [X,Y] = meshgrid(DATA.grids.x,DATA.grids.y);
REALP = 1; COMPDIM = 3; POLARPLOT = 0; SCALE = 1;
case 'xz'
XNAME = '$x$'; YNAME = '$z$';
- [Y,X] = meshgrid(DATA.z,DATA.x);
+ [Y,X] = meshgrid(DATA.grids.z,DATA.grids.x);
REALP = 1; COMPDIM = 1; SCALE = 0;
case 'yz'
XNAME = '$y$'; YNAME = '$z$';
- [Y,X] = meshgrid(DATA.z,DATA.y);
+ [Y,X] = meshgrid(DATA.grids.z,DATA.grids.y);
REALP = 1; COMPDIM = 2; SCALE = 0;
case 'kxky'
XNAME = '$k_x$'; YNAME = '$k_y$';
- [X,Y] = meshgrid(DATA.kx,DATA.ky);
+ [X,Y] = meshgrid(DATA.grids.kx,DATA.grids.ky);
REALP = 0; COMPDIM = 3; POLARPLOT = 0; SCALE = 1;
case 'kxz'
XNAME = '$k_x$'; YNAME = '$z$';
- [Y,X] = meshgrid(DATA.z,DATA.kx);
+ [Y,X] = meshgrid(DATA.grids.z,DATA.grids.kx);
REALP = 0; COMPDIM = 1; POLARPLOT = 0; SCALE = 0;
case 'kyz'
XNAME = '$k_y$'; YNAME = '$z$';
- [Y,X] = meshgrid(DATA.z,DATA.ky);
+ [Y,X] = meshgrid(DATA.grids.z,DATA.grids.ky);
REALP = 0; COMPDIM = 2; POLARPLOT = 0; SCALE = 0;
case 'sx'
XNAME = '$v_\parallel$'; YNAME = '$\mu$';
@@ -112,24 +112,24 @@ switch OPTIONS.COMP
i = OPTIONS.COMP;
compr = @(x) x(i,:,:);
if REALP
- COMPNAME = sprintf(['y=','%2.1f'],DATA.x(i));
+ COMPNAME = sprintf(['y=','%2.1f'],DATA.grids.x(i));
else
- COMPNAME = sprintf(['k_y=','%2.1f'],DATA.kx(i));
+ COMPNAME = sprintf(['k_y=','%2.1f'],DATA.grids.kx(i));
end
FIELDNAME = [LTXNAME,'(',COMPNAME,')'];
case 2
i = OPTIONS.COMP;
compr = @(x) x(:,i,:);
if REALP
- COMPNAME = sprintf(['x=','%2.1f'],DATA.y(i));
+ COMPNAME = sprintf(['x=','%2.1f'],DATA.grids.y(i));
else
- COMPNAME = sprintf(['k_x=','%2.1f'],DATA.ky(i));
+ COMPNAME = sprintf(['k_x=','%2.1f'],DATA.grids.ky(i));
end
FIELDNAME = [LTXNAME,'(',COMPNAME,')'];
case 3
i = OPTIONS.COMP;
compr = @(x) x(:,:,i);
- COMPNAME = sprintf(['z=','%2.1f'],DATA.z(i));
+ COMPNAME = sprintf(['z=','%2.1f'],DATA.grids.z(i));
FIELDNAME = [LTXNAME,'(',COMPNAME,')'];
end
end
@@ -161,7 +161,6 @@ end
switch OPTIONS.NAME
case '\phi' %ES pot
NAME = 'phi';
-% FLD_ = DATA.PHI(:,:,:,FRAMES); % data to plot
FLD_ = DATA.PHI(:,:,:,FRAMES); % data to plot
OPE_ = 1; % Operation on data
case '\psi' %EM pot
@@ -296,7 +295,7 @@ else
tmp = zeros(DATA.Nky,DATA.Nkx,Nz);
end
for it = 1:numel(FRAMES)
- for iz = 1:numel(DATA.z)
+ for iz = 1:numel(DATA.grids.z)
tmp(:,:,iz) = squeeze(process(OPE_.*FLD_(:,:,iz,it)));
end
FIELD(:,:,it) = squeeze(compr(tmp));
--
GitLab