diff --git a/matlab/compute/compute_fa_2D.m b/matlab/compute/compute_fa_2D.m
index 8d9758e0b27a106faad7bce9a9bce4058ee67f57..f2e9552e97653b21476b4888f0fd77b7aaa31f4f 100644
--- a/matlab/compute/compute_fa_2D.m
+++ b/matlab/compute/compute_fa_2D.m
@@ -1,4 +1,4 @@
-function [SS,XX,FF] = compute_fa_2D(data, options)
+function [SS,XX,FF] = compute_fa_2D(data, species, s, x, T)
%% Compute the dispersion function from the moment hierarchi decomp.
% Normalized Hermite
Hp = @(p,s) polyval(HermitePoly(p),s)./sqrt(2.^p.*factorial(p));
@@ -9,99 +9,90 @@ Lj = @(j,x) polyval(LaguerrePoly(j),x);
FaM = @(s,x) exp(-s.^2-x);
%meshgrid for efficient evaluation
-[SS, XX] = meshgrid(options.SPAR, options.XPERP);
+[SS, XX] = meshgrid(s, x);
-[~, ikx0] = min(abs(data.kx));
-[~, iky0] = min(abs(data.ky));
-kx_ = data.kx; ky_ = data.ky;
-
-switch options.SPECIES
+switch species
case 'e'
- Napj_ = data.Nepj;
- parray = double(data.Pe);
- jarray = double(data.Je);
- case 'i'
- Napj_ = data.Nipj;
- parray = double(data.Pi);
- jarray = double(data.Ji);
-end
+ Napj_ = data.Napjz(2,:,:,:,:);
-switch options.iz
- case 'avg'
- Napj_ = mean(Napj_,5);
- phi_ = mean(data.PHI,3);
- otherwise
- iz = options.iz;
- Napj_ = Napj_(:,:,:,:,iz,:);
- phi_ = data.PHI(:,:,iz);
+ case 'i'
+ Napj_ = data.Napjz(1,:,:,:,:);
end
+parray = double(data.grids.Parray);
+jarray = double(data.grids.Jarray);
+% switch options.iz
+ % case 'avg'
+ options.SHOW_FLUXSURF = 0;
+ options.SHOW_METRICS = 0;
+ options.SHOW_CURVOP = 0;
+ [~, geo_arrays] = plot_metric(data,options);
+ J = geo_arrays.Jacobian;
+ Nz = data.grids.Nz;
+ tmp_ = 0;
+ for iz = 1:Nz
+ tmp_ = tmp_ + J(iz)*Napj_(:,:,:,iz,:);
+ end
+ Napj_ = tmp_/sum(J(iz));
+ % Napj_ = mean(Napj_,4);
+ % Napj_ = Napj_(:,:,:,Nz/2+1,:);
+ % phi_ = mean(data.PHI,3);
+ % otherwise
+ % iz = options.iz;
+ % Napj_ = Napj_(:,:,:,:,iz,:);
+ % phi_ = data.PHI(:,:,iz);
+% end
% Napj_ = squeeze(Napj_);
-frames = options.T;
-for it = 1:numel(options.T)
- [~,frames(it)] = min(abs(options.T(it)-data.Ts5D));
+frames = T;
+for it = 1:numel(T)
+ [~,frames(it)] = min(abs(T(it)-data.Ts3D));
end
frames = unique(frames);
-Napj_ = mean(Napj_(:,:,:,:,frames),5);
+Napj_ = mean(Napj_(:,:,:,:,frames),5);
-% Napj_ = squeeze(Napj_);
+Napj_ = squeeze(Napj_);
Np = numel(parray); Nj = numel(jarray);
-if options.non_adiab
- for ij_ = 1:Nj
- for ikx = 1:data.Nkx
- for iky = 1:data.Nky
- kp_ = sqrt(kx_(ikx)^2 + ky_(iky)^2);
- Napj_(1,ij_,iky,ikx) = Napj_(1,ij_,iky,ikx) + kernel(ij_,kp_)*phi_(iky,ikx);
- end
- end
- end
-end
-
-if options.RMS
- FF = zeros(data.Nky,data.Nkx,numel(options.XPERP),numel(options.SPAR));
+% if options.RMS
+ FF = zeros(numel(x),numel(s));
FAM = FaM(SS,XX);
for ip_ = 1:Np
p_ = parray(ip_);
HH = Hp(p_,SS);
for ij_ = 1:Nj
- j_ = jarray(ij_);
- LL = Lj(j_,XX);
+ j_ = jarray(ij_);
+ LL = Lj(j_,XX);
HLF = HH.*LL.*FAM;
- for ikx = 1:data.Nkx
- for iky = 1:data.Nky
- FF(iky,ikx,:,:) = squeeze(FF(iky,ikx,:,:)) + Napj_(ip_,ij_,iky,ikx)*HLF;
- end
- end
+ FF = FF + Napj_(ip_,ij_)*HLF;
end
end
-else
- FF = zeros(numel(options.XPERP),numel(options.SPAR));
- FAM = FaM(SS,XX);
- for ip_ = 1:Np
- p_ = parray(ip_);
- HH = Hp(p_,SS);
- for ij_ = 1:Nj
- j_ = jarray(ij_);
- LL = Lj(j_,XX);
- FF = FF + squeeze(Napj_(ip_,ij_,ikx0,iky0))*HH.*LL.*FAM;
- end
- end
-end
+% else
+% FF = zeros(numel(options.XPERP),numel(options.SPAR));
+% FAM = FaM(SS,XX);
+% for ip_ = 1:Np
+% p_ = parray(ip_);
+% HH = Hp(p_,SS);
+% for ij_ = 1:Nj
+% j_ = jarray(ij_);
+% LL = Lj(j_,XX);
+% FF = FF + squeeze(Napj_(ip_,ij_,ikx0,iky0))*HH.*LL.*FAM;
+% end
+% end
+% end
FF = (FF.*conj(FF)); %|f_a|^2
% FF = abs(FF); %|f_a|
-if options.RMS
+% if options.RMS
% FF = squeeze(mean(mean(sqrt(FF),1),2)); %sqrt(<|f_a|^2>kx,ky)
- FF = sqrt(squeeze(mean(mean(FF,1),2))); %<|f_a|>kx,ky
-else
- FF = sqrt(squeeze(FF)); %sqrt(<|f_a|>x,y)
-end
+ FF = sqrt(FF); %<|f_a|>kx,ky
+% else
+% FF = sqrt(squeeze(FF)); %sqrt(<|f_a|>x,y)
+% end
-FF = FF./max(max(FF));
-FF = FF';
+% FF = FF./max(max(FF));
+% FF = FF';
% FF = sqrt(FF);
% FF = FF';
end
\ No newline at end of file
diff --git a/matlab/compute/fourier_GENE.m b/matlab/compute/fourier_GENE.m
new file mode 100644
index 0000000000000000000000000000000000000000..2039f29ca9ccad0f446dfbe218642b6a97687fac
--- /dev/null
+++ b/matlab/compute/fourier_GENE.m
@@ -0,0 +1,9 @@
+function [ field_c ] = fourier_GENE( field_r )
+
+%fft, symmetric as we are using real data
+spectrumKxYZ=nx*fft(fieldxyz,[],1);
+field_c=ny*fft(spectrumKxYZ,[],2);
+clear spectrumKxKyZ
+
+end
+
diff --git a/matlab/compute/ifourier_GENE.m b/matlab/compute/ifourier_GENE.m
index 1538c36dfb1cde9ae4803a82e3ebf32b7ecddaea..e05abaa4644297f91490770467a25bb20bd35166 100644
--- a/matlab/compute/ifourier_GENE.m
+++ b/matlab/compute/ifourier_GENE.m
@@ -26,22 +26,6 @@ end
spectrumXKyZ=nx*ifft(spectrumKyKxZ,[],1);
field_r=ny*ifft(spectrumXKyZ,[],2,'symmetric');
clear spectrumKxKyZ
-
-%% Adapted for HeLaZ old representation
-% [nkx,ny,nz]=size(field_c);
-% %extend to whole ky by imposing reality condition
-% nx=2*nkx-1;
-%
-% %note, we need one extra point which we set to zero for the ifft
-% spectrumKxKyZ=zeros(nx,ny,nz);
-% spectrumKxKyZ(1:nkx,:,:)=field_c(:,:,:);
-% spectrumKxKyZ((nkx+1):(nx),1,:)=conj(field_c(nkx:-1:2,1,:));
-% spectrumKxKyZ((nkx+1):(nx),2:ny,:)=conj(field_c(nkx:-1:2,ny:-1:2,:));
-%
-% %inverse fft, symmetric as we are using real data
-% spectrumXKyZ=nx*ifft(spectrumKxKyZ,[],1);
-% field_r=ny*ifft(spectrumXKyZ,[],2,'symmetric');
-% clear spectrumKxKyZ
end
diff --git a/matlab/compute/mode_growth_meter.m b/matlab/compute/mode_growth_meter.m
index ab5736b4f47592164e2158576a322d4c088c1d71..c8998c0a8c22905ed37cf65810a219fa8ccd0fb3 100644
--- a/matlab/compute/mode_growth_meter.m
+++ b/matlab/compute/mode_growth_meter.m
@@ -6,10 +6,16 @@ d = OPTIONS.fftz.flag; % To add spectral evolution of z (useful for 3d zpin
switch OPTIONS.FIELD
case 'Ni00'
FIELD = DATA.Ni00;
- fname = 'N^{00}';
+ fname = 'Ni^{00}';
+ case 'Ne00'
+ FIELD = DATA.Ne00;
+ fname = 'Ne^{00}';
case 'phi'
FIELD = DATA.PHI;
fname = '\phi';
+ case 'T_i'
+ FIELD = DATA.TEMP_I;
+ fname = 'T_i';
end
if numel(size(FIELD)) == 3
field = squeeze(FIELD);
@@ -121,8 +127,13 @@ for i = 1:2
%plot
% subplot(2+d,3,1+3*(i-1))
FIGURE.axes(1+3*(i-1)) = subplot(2+d,3,1+3*(i-1),'parent',FIGURE.fig);
+ if MODES_SELECTOR == 2
[YY,XX] = meshgrid(t,fftshift(k));
pclr = pcolor(XX,YY,abs(plt(fftshift(field,MODES_SELECTOR),1:numel(k))));set(pclr, 'edgecolor','none'); hold on;
+ else
+ [YY,XX] = meshgrid(t,(k));
+ pclr = pcolor(XX,YY,abs(plt((field),1:numel(k))));set(pclr, 'edgecolor','none'); hold on;
+ end
set(gca,'YDir','normal')
% xlim([t(1) t(end)]); %ylim([1e-5 1])
xlabel(['$',kstr,'\rho_s$']); ylabel('$t c_s /R_0$');
@@ -138,10 +149,10 @@ for i = 1:2
semilogy(t,plt(field,ikzf),'--k');
end
%plot the time window where the gr are measured
- plot(t(it1)*[1 1],[1e-10 1],'--k')
- plot(t(it2)*[1 1],[1e-10 1],'--k')
- xlim([t(1) t(end)]); %ylim([1e-5 1])
+ xline(t(it1),'--k')
+ xline(t(it2),'--k')
xlabel('$t c_s /R_0$'); ylabel(['$|',fname,'_{',kstr,'}|$']);
+ axis tight
title('Measure time window')
% subplot(2+d,3,3+3*(i-1))
diff --git a/matlab/load/compile_results_2Da.m b/matlab/load/compile_results_2Da.m
new file mode 100644
index 0000000000000000000000000000000000000000..f89f413bd464f076b155de47693d5eabb401480b
--- /dev/null
+++ b/matlab/load/compile_results_2Da.m
@@ -0,0 +1,37 @@
+function [field, Ts2D] = compile_results_2Da(DIRECTORY,JOBNUMMIN,JOBNUMMAX,fieldname)
+CONTINUE = 1;
+JOBNUM = JOBNUMMIN; JOBFOUND = 0;
+% field
+field = [];
+Ts2D = [];
+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/var2d/time'));
+ catch
+ openable = 0;
+ end
+ if openable
+ % load field %%
+ [ F, T, ~] = load_2Da_data(filename, fieldname);
+ field = cat(4,field,F);
+ Ts2D = cat(1,Ts2D,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 33fd76b5eb6086530d02faf9553cbc2b7dfddc2f..8ef6048ccd5d717c9c2f64699ddbeb84b77c267f 100644
--- a/matlab/load/compile_results_low_mem.m
+++ b/matlab/load/compile_results_low_mem.m
@@ -29,7 +29,7 @@ while(CONTINUE)
DATA.outfilenames{ii} = filename;
%test if it is corrupted or currently running
try
- openable = ~isempty(h5read(filename,'/data/var3d/time'));
+ openable = ~isempty(h5read(filename,'/data/var0d/time'));
catch
openable = 0;
end
@@ -38,13 +38,25 @@ while(CONTINUE)
fprintf('Loading ID %.2d (%s)\n',JOBNUM,filename);
% loading input parameters
DATA.(sprintf('fort_%.2d',JOBNUM)) = struct();
- DATA.(sprintf('fort_%.2d',JOBNUM)) = read_namelist(fortname);
+ % cp the STDIN.XX from the output file to the directory to load
+ tmpfort = [fortname,'.tmp'];
+ fid = fopen(tmpfort,'wt');
+ input_string = h5read(filename,sprintf('/files/STDIN.%.2d',JOBNUM));
+ fprintf(fid, input_string);
+ fclose(fid);
+ DATA.(sprintf('fort_%.2d',JOBNUM)) = read_namelist(tmpfort);
+ system(['rm ',tmpfort]); % clean it
% Loading from output file
CPUTIME = h5readatt(filename,'/data/input','cpu_time');
DT_SIM = h5readatt(filename,'/data/input/basic','dt');
[Parray, Jarray, kx, ky, z] = load_grid_data(filename);
- W_GAMMA = strcmp(h5readatt(filename,'/data/input/diagnostics','write_gamma'),'y');
- W_HF = strcmp(h5readatt(filename,'/data/input/diagnostics','write_hf' ),'y');
+ try
+ W_GAMMA = strcmp(h5readatt(filename,'/data/input/diagnostics','write_gamma'),'y');
+ W_HF = strcmp(h5readatt(filename,'/data/input/diagnostics','write_hf' ),'y');
+ catch
+ W_GAMMA = strcmp(h5readatt(filename,'/data/input/diag_par','write_gamma'),'y');
+ W_HF = strcmp(h5readatt(filename,'/data/input/diag_par','write_hf' ),'y');
+ end
KIN_E = strcmp(h5readatt(filename,'/data/input/model', 'ADIAB_E' ),'n');
BETA = h5readatt(filename,'/data/input/model','beta');
diff --git a/matlab/load/load_2D_data.m b/matlab/load/load_2D_data.m
index 67da60c9dd1931e73cde73c7520ecb6370a0c153..6ca82e8430e399ec31c1d9feb74d31e3d2203a6e 100644
--- a/matlab/load/load_2D_data.m
+++ b/matlab/load/load_2D_data.m
@@ -1,5 +1,5 @@
function [ data, time, dt ] = load_2D_data( filename, variablename )
-%LOAD_3D_DATA load a 3D variable stored in a hdf5 result file from HeLaZ
+%LOAD_2D_DATA load a 2D variable stored in a hdf5 result file
time = h5read(filename,'/data/var2d/time');
dt = h5readatt(filename,'/data/input/basic','dt');
cstart= h5readatt(filename,'/data/input/basic','start_iframe3d');
diff --git a/matlab/load/load_2Da_data.m b/matlab/load/load_2Da_data.m
new file mode 100644
index 0000000000000000000000000000000000000000..3cdd94bfb0f0db936ee9745c89e8755843deea76
--- /dev/null
+++ b/matlab/load/load_2Da_data.m
@@ -0,0 +1,43 @@
+function [ data, time, dt ] = load_2Da_data( filename, variablename )
+%LOAD_2Da_DATA load a 2D variable stored in a hdf5 result file
+ time = h5read(filename,'/data/var2d/time');
+ dt = h5readatt(filename,'/data/input/basic','dt');
+ cstart= h5readatt(filename,'/data/input/basic','start_iframe2d');
+ 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/var2d/',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/var2d/',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_multiple_outputs_marconi.m b/matlab/load/load_multiple_outputs_marconi.m
deleted file mode 100644
index 1fdd3b8608236ae4a056034abf7852fc50a814a0..0000000000000000000000000000000000000000
--- a/matlab/load/load_multiple_outputs_marconi.m
+++ /dev/null
@@ -1,6 +0,0 @@
-addpath(genpath('../matlab')) % ... add
-load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/300x150_L_120_P_8_J_4_eta_0.6_nu_1e-01_SGGK_mu_0e+00/')
-load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_B/300x150_L_120_P_8_J_4_eta_0.6_nu_5e-01_SGGK_mu_0e+00/')
-load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/300x150_L_120_P_8_J_4_eta_0.6_nu_1e-01_PAGK_mu_0e+00/')
-load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/500x500_L_120_P_4_J_2_eta_0.6_nu_1e-01_DGGK_mu_0e+00/')
-load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/500x500_L_120_P_2_J_1_eta_0.6_nu_1e-01_DGGK_mu_0e+00/')
diff --git a/matlab/load/load_params.m b/matlab/load/load_params.m
index 967baf3c4587ab4217c2763ad1c0260becddb13d..85d7c00d0ba3f29c41297abb58c7a6e107cf8b06 100644
--- a/matlab/load/load_params.m
+++ b/matlab/load/load_params.m
@@ -51,12 +51,19 @@ function [DATA] = load_params(DATA,filename)
catch
DATA.inputs.ADIAB_I = 0;
end
- DATA.inputs.W_GAMMA = h5readatt(filename,'/data/input/diagnostics','write_gamma') == 'y';
- DATA.inputs.W_PHI = h5readatt(filename,'/data/input/diagnostics','write_phi') == 'y';
- DATA.inputs.W_NA00 = h5readatt(filename,'/data/input/diagnostics','write_Na00') == 'y';
- DATA.inputs.W_NAPJ = h5readatt(filename,'/data/input/diagnostics','write_Napj') == 'y';
- DATA.inputs.W_SAPJ = h5readatt(filename,'/data/input/diagnostics','write_Sapj') == 'y';
-
+ try
+ DATA.inputs.W_GAMMA = h5readatt(filename,'/data/input/diagnostics','write_gamma') == 'y';
+ DATA.inputs.W_PHI = h5readatt(filename,'/data/input/diagnostics','write_phi') == 'y';
+ DATA.inputs.W_NA00 = h5readatt(filename,'/data/input/diagnostics','write_Na00') == 'y';
+ DATA.inputs.W_NAPJ = h5readatt(filename,'/data/input/diagnostics','write_Napj') == 'y';
+ DATA.inputs.W_SAPJ = h5readatt(filename,'/data/input/diagnostics','write_Sapj') == 'y';
+ catch
+ 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';
+ end
% Species dependent parameters
DATA.inputs.sigma = zeros(1,DATA.inputs.Na);
DATA.inputs.tau = zeros(1,DATA.inputs.Na);
diff --git a/matlab/load/readFortranNamelist.m b/matlab/load/readFortranNamelist.m
deleted file mode 100644
index 5dbbda9337ba82b2ff2674d9d30573dc51ba30e5..0000000000000000000000000000000000000000
--- a/matlab/load/readFortranNamelist.m
+++ /dev/null
@@ -1,213 +0,0 @@
-function S = read_namelist(filename)
-% S = READ_NAMELIST(FILENAME) returns the struct S containg namelists and
-% variables in the file FILENAME organised in hierachical way:
-%
-% |--VAR1
-% |--VAR2
-% |-- NMLST_A--|...
-% | |--VARNa
-% |
-% | |--VAR1
-% |-- NMLST_B--|--VAR2
-% | |...
-% S --| ... |--VARNb
-% |
-% | |--VAR1
-% |-- NMLST_M--|--VAR2
-% |...
-% |--VARNm
-%
-% Note: The function can read multidimensioal variables as well. The
-% function assumes that there is no more than one namelist section per
-% line. At this time there is no syntax checking functionality so the
-% function will crash in case of errors.
-%
-% Example:
-% NMLST = read_namelist('OPTIONS.nam');
-% NMLST.NAM_FRAC.XUNIF_NATURE = 0.1;
-% write_namelist(NMlST, 'MOD_OPTIONS.nam');
-%
-% Written by: Darien Pardinas Diaz (darien.pardinas-diaz@monash.edu)
-% Version: 1.0
-% Date: 16 Dec 2011
-S = struct();
-% Open and read the text file containing the namelists
-fid = fopen(filename,'r');
-c = 0;
-lines = cell(1);
-% Read all the text lines in namelist file
-while ~feof(fid)
- line = fgetl(fid);
- % Remove comments if any on the line
- idx = find(line == '!');
- if ~isempty(idx),
- line = line(1:idx(1)-1);
- end
- if ~isempty(line),
- c = c + 1;
- lines{c} = line;
- end
-end
-fclose(fid);
-i = 0;
-while i < c;
- % Find a record
- i = i + 1;
- line = lines{i};
- idx = find(line == '&');
- if ~isempty(idx), % i.e. a namelist start
- line = line(idx(1) + 1:end);
- % find next space
- idx = find(line == ' ');
- if ~isempty(idx),
- namelst = line(1:idx(1) - 1);
- line = line(idx(1) + 1:end);
- else
- namelst = line;
- line = [];
- end
- nmlst_bdy = [];
- idx = strfind(line,'/');
- % Get the variable specification section
- while isempty(idx) && i < c,
- nmlst_bdy = [nmlst_bdy ' ' line];
- i = i + 1;
- line = lines{i};
- idx = strfind(line,'/');
- end
- if ~isempty(idx) && idx(1) > 1,
- nmlst_bdy = [nmlst_bdy ' ' line];
- end
- % Parse current namelist (set of variables)
- S.(namelst) = parse_namelist(nmlst_bdy);
- end
-end
-function S = parse_namelist(strng)
-% Internal function to parse the body text of a namelist section.
-% Limitations: the following patterns are prohibited inside the literal
-% strings: '.t.' '.f.' '.true.' '.false.' '(:)'
-% Get all .true., .t. and .false., .f. to T and F
-strng = regexprep(strng,'\.true\.' ,'T','ignorecase');
-strng = regexprep(strng,'\.false\.','F','ignorecase');
-strng = regexprep(strng,'\.t\.','T','ignorecase');
-strng = regexprep(strng,'\.f\.','F','ignorecase');
-% Make evaluable the (:) expression in MATLAB if any
-strng = regexprep(strng, '\(:\)', '(1,:)');
-[strng, islit] = parse_literal_strings([strng ' ']);
-% Find the position of all the '='
-eq_idx = find(strng == '=');
-nvars = length(eq_idx);
-arg_start = eq_idx + 1;
-arg_end = zeros(size(eq_idx));
-vars = cell(nvars,1);
-S = struct;
-% Loop through every variable
-for k = 1:nvars,
- i = eq_idx(k) - 1;
- % Move to the left and discard blank spaces
- while strng(i) == ' ', i = i - 1; end
- % Now we are over the variable name or closing parentesis
- j = i;
- if strng(i) == ')',
- while strng(i) ~= '(', i = i - 1; end
- i = i - 1;
- % Move to the left and discard any possible blank spaces
- while strng(i) == ' ', i = i - 1; end
- end
-
- % Now we are over the last character of the variable name
- while strng(i) ~= ' ', i = i - 1; end
-
- if k > 1, arg_end(k - 1) = i; end
- vars{k} = ['S.' strng(i + 1: j)];
-end
-arg_end(end) = length(strng);
-% This variables are used in the eval function to evaluate True/False,
-% so don't remove it!
-T = '.true.';
-F = '.false.';
-% Loop through every variable guess varible type
-for k = 1:nvars,
- arg = strng(arg_start(k):arg_end(k));
- arglit = islit(arg_start(k):arg_end(k))';
-
- % Remove commas in non literal string...
- commas = ~arglit & arg == ',';
- if any(commas)
- arg(commas) = ' ';
- end
-
- if any(arglit),
- % We are parsing a variable that is literal string
- arg = ['{' arg '};'];
- elseif ~isempty(find( arg == 'T' | arg == 'F', 1)),
- % We are parsing a boolean variable
- arg = ['{' arg '};'];
- else
- % We are parsing a numerical array
- arg = ['[' arg '];'];
- end
- % Eval the modified syntax in Matlab
- eval([vars{k} ' = ' arg]);
-end
-function [strng, is_lit] = parse_literal_strings(strng)
-% Parse the literal declarations of strings and change to Matlab syntax
-len = length(strng);
-add_squote = []; % Positions to add a scape single quote on syntax
-rem_dquote = []; % Positions to remove a double quote scape on syntax
-i = 1;
-while i < len,
- if strng(i) == '''', % Opening string with single quote...
- i = i + 1;
- while i < len && strng(i) ~= '''' || strcmp(strng(i:i+1),'''''') ,
- i = i + 1;
- if strcmp(strng(i-1:i),''''''),
- i = i + 1;
- end
- end
- end
- if strng(i) == '"', % Opening string with double quote...
- strng(i) = ''''; % Change to single quote
- i = i + 1;
- while strng(i) ~= '"' || strcmp(strng(i:i+1),'""') && i < len,
- % Check for a possible sigle quote here
- if strng(i) == '''',
- add_squote = [add_squote i];
- end
- i = i + 1;
- if strcmp(strng(i-1:i),'""'),
- rem_dquote = [rem_dquote i-1];
- i = i + 1;
- end
- end
- strng(i) = ''''; % Change to single quote
- end
- i = i + 1;
-end
-for i = 1:length(add_squote);
- strng = [strng(1:add_squote(i)) strng(add_squote(i):end)];
-end
-for i = 1:length(rem_dquote);
- strng = [strng(1:rem_dquote(i)-1) strng(rem_squote(i)+1:end)];
-end
-% Now everything should be in Matlab string syntax
-% Classify syntax as literal or regular expression
-i = 1;
-len = length(strng);
-is_lit = zeros(len,1);
-while i < len,
- if strng(i) == '''', % Opening string with single quote...
- is_lit(i) = 1;
- i = i + 1;
- while i < len && strng(i) ~= '''' || strcmp(strng(i:i+1),''''''),
- is_lit(i) = 1;
- i = i + 1;
- if strcmp(strng(i-1:i),''''''),
- is_lit(i) = 1;
- i = i + 1;
- end
- end
- is_lit(i) = 1;
- end
- i = i + 1;
-end
diff --git a/matlab/load/read_namelist.m b/matlab/load/read_namelist.m
index 3d49e6d0406d2532eaa3fc51ee3f7b2fe84d7631..954dbc9af151a177268b600963ef51968cade8b8 100644
--- a/matlab/load/read_namelist.m
+++ b/matlab/load/read_namelist.m
@@ -40,10 +40,10 @@ while ~feof(fid)
line = fgetl(fid);
% Remove comments if any on the line
idx = find(line == '!');
- if ~isempty(idx),
+ if ~isempty(idx)
line = line(1:idx(1)-1);
end
- if ~isempty(line),
+ if ~isempty(line)
c = c + 1;
lines{c} = line;
end
diff --git a/matlab/plot/Ursula_Meier.m b/matlab/plot/Ursula_Meier.m
new file mode 100644
index 0000000000000000000000000000000000000000..affabdea29979e8a4f585ad3433beb13caf1df26
--- /dev/null
+++ b/matlab/plot/Ursula_Meier.m
@@ -0,0 +1,106 @@
+function [filename] = Ursula_Meier(MOVIE)
+ FPS = 16; BWR = 0; ncolor_level = 128; format = '.gif'; SAT = 0.75;
+ INTERP = 1;
+ if isfield(MOVIE,'FPS')
+ FPS = MOVIE.FPS;
+ end
+ if isfield(MOVIE,'BWR')
+ BWR = MOVIE.BWR;
+ end
+ if isfield(MOVIE,'format')
+ format = MOVIE.format;
+ end
+ if isfield(MOVIE,'SAT')
+ SAT = MOVIE.SAT;
+ end
+ if isfield(MOVIE,'INTERP')
+ INTERP = MOVIE.INTERP;
+ end
+ DELAY = 1/MOVIE.FPS;
+ FILENAME = MOVIE.FILENAME; FIELDNAME = MOVIE.FIELDNAME;
+ X = MOVIE.X; Y = MOVIE.Y; T = MOVIE.T; F = MOVIE.F;
+ XNAME = MOVIE.XNAME; YNAME = MOVIE.YNAME;
+ switch format
+ case '.avi'
+ vidfile = VideoWriter(FILENAME,'Uncompressed AVI');
+ vidfile.FrameRate = FPS;
+ open(vidfile);
+ end
+ % Set colormap boundaries
+ hmax = max(max(max(F(:,:,:))));
+ hmin = min(min(min(F(:,:,:))));
+ if hmax == hmin
+ disp('Warning : h = hmin = hmax = const')
+ else
+ % Setup figure frame
+ fig = figure('Color','white');%,'Position', toplot.DIMENSIONS.*[0.5 0.5 1.0 1]);
+ pcolor(X,Y,F(:,:,1)); % to set up
+ if BWR
+ colormap(bluewhitered)
+ else
+ colormap(gray)
+ end
+ clim(SAT*[-1 1])
+ axis equal
+ axis tight % this ensures that getframe() returns a consistent size
+ if INTERP
+ shading interp;
+ end
+ in = 1;
+ nbytes = fprintf(2,'frame %d/%d',in,numel(F(1,1,:)));
+ for n = 1:numel(T) % loop over selected frames
+ scale = max(max(abs(F(:,:,n)))); % Scaling to normalize
+ pclr = pcolor(X,Y,F(:,:,n)/scale); % frame plot\
+ clim([-1,1]);
+ if INTERP
+ shading interp;
+ end
+ set(pclr, 'edgecolor','none'); %pbaspect(toplot.ASPECT);
+ if BWR
+ colormap(bluewhitered)
+ else
+ colormap(gray)
+ end
+ clim(SAT*[-1 1])
+ title([FIELDNAME,', $t \approx$', sprintf('%.3d',ceil(T(n)))...
+ ,', scaling = ',sprintf('%.1e',scale)]);
+
+ xlabel(XNAME); ylabel(YNAME); %colorbar;
+ axis tight equal
+ drawnow
+ % Capture the plot as an image
+ frame = getframe(fig);
+ switch format
+ case '.gif'
+ im = frame2im(frame);
+ [imind,cm] = rgb2ind(im,ncolor_level);
+ % Write to the GIF File
+ if in == 1
+ imwrite(imind,cm,FILENAME,'gif', 'Loopcount',inf);
+ else
+ imwrite(imind,cm,FILENAME,'gif','WriteMode','append', 'DelayTime',DELAY);
+ end
+ case '.avi'
+ writeVideo(vidfile,frame);
+ otherwise
+ disp('Unknown format');
+ break
+ end
+ % terminal info
+ while nbytes > 0
+ fprintf('\b')
+ nbytes = nbytes - 1;
+ end
+ nbytes = fprintf(2,'frame %d/%d',n,numel(F(1,1,:)));
+ in = in + 1;
+ end
+ disp(' ')
+ switch format
+ case '.gif'
+ disp(['Gif saved @ : ',FILENAME])
+ case '.avi'
+ disp(['Video saved @ : ',FILENAME])
+ close(vidfile);
+ end
+ end
+end
\ No newline at end of file
diff --git a/matlab/plot/cinecita.m b/matlab/plot/cinecita.m
new file mode 100644
index 0000000000000000000000000000000000000000..fcfcdd03609bd2a2667102006d974a1b2d372e22
--- /dev/null
+++ b/matlab/plot/cinecita.m
@@ -0,0 +1,69 @@
+function [ ] = cinecita(DATADIR,J0,J1,fieldname,plan,save)
+
+data = {};
+[data] = compile_results_low_mem(data,DATADIR,J0,J1);
+%% Select the field
+SKIP_COMP = 0; % Turned on only for 2D plots
+OPE_ = 1; % Operation on data
+switch fieldname
+ case 'phi' %ES pot
+ [FIELD,TIME] = compile_results_3D(DATADIR,J0,J1,'phi');
+ ltxname = '\phi';
+ case 'phi_obmp' %ES pot
+ [FIELD,TIME] = compile_results_2D(DATADIR,J0,J1,'phi_obmp');
+ ltxname = '\phi_{z=0}';
+ SKIP_COMP = 1;
+ case '\psi' %EM pot
+ [FIELD,TIME] = compile_results_3D(DATADIR,J0,J1,'psi');
+ ltxname = '\psi';
+ case 'phi_nz' % non-zonal ES pot
+ [FIELD,TIME] = compile_results_3D(DATADIR,J0,J1,'phi');
+ ltxname = '\phi^{NZ}';
+ OPE_ = (KY~=0);
+ case 'vE_y' % y-comp of ExB velocity
+ [FIELD,TIME] = compile_results_3D(DATADIR,J0,J1,'phi');
+ ltxname = 'v_{Ey}';
+ OPE_ = -1i*KX;
+ case 'vE_x' % x-comp of ExB velocity
+ [FIELD,TIME] = compile_results_3D(DATADIR,J0,J1,'phi');
+ ltxname = 'v_{Ex}';
+ OPE_ = -1i*KY;
+ case 'sE_y' % y-comp of ExB shear
+ [FIELD,TIME] = compile_results_3D(DATADIR,J0,J1,'phi');
+ ltxname = 's_{Ey}';
+ OPE_ = KX.^2;
+ case 'sE_x' % x-comp of ExB shear
+ [FIELD,TIME] = compile_results_3D(DATADIR,J0,J1,'phi');
+ ltxname = 's_{Ex}';
+ OPE_ = KY.^2;
+ case 'wz' % ES pot vorticity
+ [FIELD,TIME] = compile_results_3D(DATADIR,J0,J1,'phi');
+ ltxname = '\omega_z';
+ OPE_ = -(KX.^2+KY.^2);
+ otherwise
+ disp('Fieldname not recognized');
+ return
+end
+%% Setup grids
+% [KX, KY] = meshgrid(DATA.grids.kx, DATA.grids.ky);
+Nx = data.grids.Nx; Ny = data.grids.Ny; Nz = data.grids.Nz; Nt = numel(TIME);
+% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Options
+options.INTERP = 0;
+options.POLARPLOT = 0;
+options.BWR = 1; % bluewhitered plot or gray
+options.CLIMAUTO = 1; % adjust the colormap auto
+options.NAME = ltxname;
+options.PLAN = plan;
+options.COMP = 1;
+options.TIME = data.Time(1:1:end);
+data.EPS = 0.1;
+data.a = data.EPS * 2000;
+options.RESOLUTION = 256;
+create_film(data,options,'.gif')
+
+if ~save
+ system(['rm ',FILENAME]);
+end
+
+end
\ No newline at end of file
diff --git a/matlab/plot/cinecita_2D.m b/matlab/plot/cinecita_2D.m
new file mode 100644
index 0000000000000000000000000000000000000000..5a6f0f4f5498ac4f3cf3f5a5748ddd1e4857e5e4
--- /dev/null
+++ b/matlab/plot/cinecita_2D.m
@@ -0,0 +1,136 @@
+function [MOVIE] = cinecita_2D(DATADIR,fieldname,varargin)
+DATADIR = [DATADIR,'/'];
+J0 = 0; J1 = 20; nskip = 1; FOURIER = 0; SAT = 0.75;
+PLAY = 1;
+if numel(varargin) > 0
+ if isfield(varargin{1},'J0')
+ J0 = varargin{1}.J0;
+ end
+ if isfield(varargin{1},'J1')
+ J1 = varargin{1}.J1;
+ end
+ if isfield(varargin{1},'nskip')
+ nskip = varargin{1}.nskip;
+ end
+ if isfield(varargin{1},'FOURIER')
+ FOURIER = varargin{1}.FOURIER;
+ end
+ if isfield(varargin{1},'SAT')
+ SAT = varargin{1}.SAT;
+ end
+ if isfield(varargin{1},'PLAY')
+ PLAY = varargin{1}.PLAY;
+ end
+end
+%
+data = {};
+data = compile_results_low_mem(data,DATADIR,J0,J1);
+data = load_params(data,data.outfilenames{1});
+Nx = data.grids.Nx; Ny = data.grids.Ny; Na = data.inputs.Na;
+[KX,KY] = meshgrid(data.grids.kx,data.grids.ky);
+format = '.gif';
+%% Select the field
+OPE_ = 1; % Operation on data
+Nspec = 0; % Load a specie
+switch fieldname
+case '\phi' %ES pot
+ toload = 'phi';
+ ltxname = '\phi';
+ vname = 'phi';
+case '\phi_{zf}' %ES pot
+ toload = 'phi';
+ OPE_ = KY==0; % Operation on data
+ ltxname = '\phi_{zf}';
+ vname = 'phizf';
+case '\psi' %EM pot
+ toload = 'psi';
+ ltxname = '\psi';
+ vname = 'psi';
+case 'n_i' % ion density
+ toload = 'dens';
+ Nspec = 1;
+ ltxname = 'n_i';
+ vname = 'dens_i';
+case 'n_e' % ion density
+ toload = 'dens';
+ Nspec = 2;
+ ltxname = 'n_e';
+ vname = 'dens_e';
+case 'v_{Ex}' %ES pot
+ toload = 'phi';
+ OPE_ = -1i*KY; % Operation on data
+ ltxname= 'v_{Ex}';
+ vname = 'vEx';
+case 'v_{Ey}' %ES pot
+ toload = 'phi';
+ OPE_ = 1i*KX; % Operation on data
+ ltxname= 'v_{Ey}';
+ vname = 'vEy';
+case '\delta B_x'
+ toload = 'psi';
+ OPE_ = -1i*KY; % Operation on data
+ ltxname = '\delta B_x';
+ vname = 'dBx';
+case '\delta B_y'
+ toload = 'psi';
+ OPE_ = 1i*KX; % Operation on data
+ ltxname = '\delta B_y';
+ vname = 'dBy';
+otherwise
+ toload = fieldname;
+ OPE_ = 1;
+end
+
+if Nspec
+ [FIELD,TIME] = compile_results_2Da(DATADIR,J0,J1,[toload,'_obmp']);
+ FIELD = squeeze(FIELD(Nspec,:,:,:));
+else
+ [FIELD,TIME] = compile_results_2D(DATADIR,J0,J1,[toload,'_obmp']);
+end
+
+TIME = TIME(1:nskip:end); FIELD = FIELD(:,:,1:nskip:end);
+
+switch ~FOURIER
+ case 1 % Real space plot
+ XNAME = '$x$'; YNAME = '$y$'; plane ='xy';
+ [X,Y] = meshgrid(data.grids.x,data.grids.y);
+ INTERP = 1;
+ process = @(x) real(fftshift(ifourier_GENE(x)));
+ shift_x = @(x) x;
+ shift_y = @(x) x;
+ case 0 % Frequencies plot
+ XNAME = '$k_x$'; YNAME = '$k_y$'; plane ='kxky';
+ [X,Y] = meshgrid(data.grids.kx,data.grids.ky);
+ INTERP = 0;
+ process = @(x) abs(fftshift(x,2));
+ shift_x = @(x) fftshift(x,2);
+ shift_y = @(x) fftshift(x,2);
+end
+
+for it = 1:numel(TIME)
+ tmp = squeeze(OPE_.*FIELD(:,:,it));
+ F(:,:,it) = squeeze(process(tmp));
+end
+
+X = shift_x(X);
+Y = shift_y(Y);
+FILENAME = [vname,'_',plane,'_obmp'];
+FILENAME = [DATADIR,FILENAME,format];
+FIELDNAME = ['$',ltxname,'$'];
+MOVIE.F = F;
+MOVIE.X = X;
+MOVIE.Y = Y;
+MOVIE.T = TIME;
+MOVIE.XNAME = XNAME;
+MOVIE.YNAME = YNAME;
+MOVIE.FPS = 16;
+MOVIE.BWR = 0;
+MOVIE.SAT = 0.75;
+MOVIE.FIELDNAME = FIELDNAME;
+MOVIE.FILENAME = FILENAME;
+% Shoot the movie
+if PLAY
+ Ursula_Meier(MOVIE)
+end
+
+end
\ No newline at end of file
diff --git a/matlab/plot/create_film.m b/matlab/plot/create_film.m
index 6becccec8c45ab7c06a95c19e033bbf3d28363c9..ca56570d89752bf52d946ac16a07bc43a0a40f99 100644
--- a/matlab/plot/create_film.m
+++ b/matlab/plot/create_film.m
@@ -1,16 +1,22 @@
function create_film(DATA,OPTIONS,format)
%% Plot options
-FPS = 16; DELAY = 1/FPS;
+FPS = OPTIONS.FPS; DELAY = 1/FPS;
BWR = OPTIONS.BWR; NORMALIZED = 1;
-if ~strcmp(OPTIONS.PLAN,'sx')
- T = DATA.Ts3D;
-else
- T = DATA.Ts5D;
-end
%% Processing
switch OPTIONS.PLAN
case 'RZ'
toplot = poloidal_plot(DATA,OPTIONS);
+ case '3D'
+ OPTIONS.PLAN = 'xy';
+ [~,OPTIONS.COMP] = min(abs(OPTIONS.XYZ(3) - DATA.grids.z));
+ toplot = process_field(DATA,OPTIONS);
+ OPTIONS.PLAN = 'xz';
+ [~,OPTIONS.COMP] = min(abs(OPTIONS.XYZ(2) - DATA.grids.y));
+ toplot_xz = process_field(DATA,OPTIONS);
+ OPTIONS.PLAN = 'yz';
+ [~,OPTIONS.COMP] = min(abs(OPTIONS.XYZ(1) - DATA.grids.x));
+ toplot_yz = process_field(DATA,OPTIONS);
+ OPTIONS.PLAN = '3D';
otherwise
toplot = process_field(DATA,OPTIONS);
end
@@ -20,7 +26,7 @@ XNAME = ['$',toplot.XNAME,'$'];
YNAME = ['$',toplot.YNAME,'$'];
FIELDNAME = ['$',toplot.FIELDNAME,'$'];
FIELD = toplot.FIELD; X = toplot.X; Y = toplot.Y;
-FRAMES = toplot.FRAMES;
+TIME = toplot.TIME;
switch format
case '.avi'
vidfile = VideoWriter(FILENAME,'Uncompressed AVI');
@@ -36,9 +42,22 @@ if hmax == hmin
disp('Warning : h = hmin = hmax = const')
else
% Setup figure frame
-fig = figure('Color','white');%,'Position', toplot.DIMENSIONS.*[0.5 0.5 1.0 1]);
if ~strcmp(OPTIONS.PLAN,'sx')
- pcolor(X,Y,FIELD(:,:,1)); % to set up
+ if ~strcmp(OPTIONS.PLAN,'3D')
+ fig = figure('Color','white');
+ pclr = pcolor(X,Y,FIELD(:,:,1)/scale); % frame plot\
+ set(pclr, 'edgecolor','none'); %pbaspect(toplot.ASPECT);
+ else
+ fig = figure('Color','white','Position', 1024*[0.5 0.5 0.5 1]);
+ s = surface(toplot.X,toplot.Y,OPTIONS.XYZ(3)+0*toplot.X,FIELD(:,:,1)/scale); hold on;
+ s.EdgeColor = 'none';
+ s = surface(toplot_xz.X,OPTIONS.XYZ(2)+0*toplot_xz.X,toplot_xz.Y,toplot_xz.FIELD(:,:,1)./scale);
+ s.EdgeColor = 'none';
+ s = surface(OPTIONS.XYZ(1)+0*toplot_yz.X,toplot_yz.X,toplot_yz.Y,toplot_yz.FIELD(:,:,1)./scale);
+ s.EdgeColor = 'none';
+ zlabel('z');
+ view([1 -1 0.25])
+ end
if BWR
colormap(bluewhitered)
else
@@ -52,30 +71,31 @@ fig = figure('Color','white');%,'Position', toplot.DIMENSIONS.*[0.5 0.5 1.0 1])
shading interp;
end
else
+ fig = figure('Color','white');
contour(toplot.X,toplot.Y,FIELD(:,:,1),128);
end
in = 1;
nbytes = fprintf(2,'frame %d/%d',in,numel(FIELD(1,1,:)));
- for n = 1:numel(FRAMES) % loop over selected frames
+ for n = 1:numel(TIME) % loop over selected frames
scale = max(max(abs(FIELD(:,:,n)))); % Scaling to normalize
if ~strcmp(OPTIONS.PLAN,'sx')
- if NORMALIZED
+ if ~strcmp(OPTIONS.PLAN,'3D')
pclr = pcolor(X,Y,FIELD(:,:,n)/scale); % frame plot\
- caxis([-1,1]);
+ set(pclr, 'edgecolor','none'); %pbaspect(toplot.ASPECT);
else
- pclr = pcolor(X,Y,FIELD(:,:,n)); % frame plot\
- if CONST_CMAP
- caxis([-1,1]*max(abs([hmin hmax]))); % adaptive color map
- else
- caxis([-1,1]*scale); % adaptive color map
- end
- title([FIELDNAME,', $t \approx$', sprintf('%.3d',ceil(T(n)))...
- ,', scale = ',sprintf('%.1e',scale)]);
+ s = surface(toplot.X,toplot.Y,OPTIONS.XYZ(3)+0*toplot.X,FIELD(:,:,n)/scale); hold on;
+ s.EdgeColor = 'none';
+ s = surface(toplot_xz.X,OPTIONS.XYZ(2)+0*toplot_xz.X,toplot_xz.Y,toplot_xz.FIELD(:,:,n)./scale);
+ s.EdgeColor = 'none';
+ s = surface(OPTIONS.XYZ(1)+0*toplot_yz.X,toplot_yz.X,toplot_yz.Y,toplot_yz.FIELD(:,:,n)./scale);
+ s.EdgeColor = 'none';
+ zlabel('z');
+ view([1 -1 0.25])
end
+ clim([-1,1]);
if toplot.INTERP
shading interp;
end
- set(pclr, 'edgecolor','none'); %pbaspect(toplot.ASPECT);
if BWR
colormap(bluewhitered)
else
@@ -87,10 +107,14 @@ fig = figure('Color','white');%,'Position', toplot.DIMENSIONS.*[0.5 0.5 1.0 1])
else % show velocity distr.
contour(toplot.X,toplot.Y,FIELD(:,:,n)/scale,128);
end
- title([FIELDNAME,', $t \approx$', sprintf('%.3d',ceil(T(FRAMES(n))))...
- ,', scaling = ',sprintf('%.1e',scale)]);
-
- xlabel(XNAME); ylabel(YNAME); %colorbar;
+ if ~OPTIONS.RMAXIS
+ title([FIELDNAME,', $t \approx$', sprintf('%.3d',ceil(TIME(n)))...
+ ,', scaling = ',sprintf('%.1e',scale)]);
+ xlabel(XNAME); ylabel(YNAME); %colorbar;
+ else
+ axis off
+ axis tight
+ end
drawnow
% Capture the plot as an image
frame = getframe(fig);
diff --git a/matlab/plot/photomaton.m b/matlab/plot/photomaton.m
index 67b02747b9ca4d0f12ce50d12ff2a476018ecc3e..d563030d9b1ae66ad1ab1cb54ac69dec98bd3b4d 100644
--- a/matlab/plot/photomaton.m
+++ b/matlab/plot/photomaton.m
@@ -4,11 +4,26 @@ function [ FIGURE ] = photomaton( DATA,OPTIONS )
switch OPTIONS.PLAN
case 'RZ'
toplot = poloidal_plot(DATA,OPTIONS);
+ case '3D'
+ OPTIONS.PLAN = 'xy';
+ [~,OPTIONS.COMP] = min(abs(OPTIONS.XYZ(3) - DATA.grids.z));
+ toplot = process_field(DATA,OPTIONS);
+ OPTIONS.PLAN = 'xz';
+ [~,OPTIONS.COMP] = min(abs(OPTIONS.XYZ(2) - DATA.grids.y));
+ toplot_xz = process_field(DATA,OPTIONS);
+ OPTIONS.PLAN = 'yz';
+ [~,OPTIONS.COMP] = min(abs(OPTIONS.XYZ(1) - DATA.grids.x));
+ toplot_yz = process_field(DATA,OPTIONS);
+ OPTIONS.PLAN = '3D';
otherwise
toplot = process_field(DATA,OPTIONS);
end
FNAME = toplot.FILENAME;
FRAMES = toplot.FRAMES;
+if OPTIONS.TAVG
+ toplot.FIELD = mean(toplot.FIELD,3);
+ FRAMES = FRAMES(1);
+end
Nframes= numel(FRAMES);
Nrows = ceil(Nframes/4);
Ncols = ceil(Nframes/Nrows);
@@ -29,17 +44,31 @@ FIGURE.fig = figure; %set(gcf, 'Position', toplot.DIMENSIONS.*[1 1 Ncols Nrows]
end
if ~strcmp(OPTIONS.PLAN,'sx')
tshot = DATA.Ts3D(FRAMES(i_));
- pclr = pcolor(toplot.X,toplot.Y,toplot.FIELD(:,:,i_)./scale); set(pclr, 'edgecolor','none');
- if OPTIONS.AXISEQUAL
- pbaspect(toplot.ASPECT)
- end
- if ~strcmp(OPTIONS.PLAN,'kxky')
- clim([-1,1]*frame_max/scale);
- colormap(bluewhitered);
- if OPTIONS.INTERP
- shading interp;
- end
- end
+ if ~strcmp(OPTIONS.PLAN,'3D')
+ pclr = pcolor(toplot.X,toplot.Y,toplot.FIELD(:,:,i_)./scale); set(pclr, 'edgecolor','none');
+ else
+ % xy plane
+ s = surface(toplot.X,toplot.Y,OPTIONS.XYZ(3)+0*toplot.Y,toplot.FIELD(:,:,i_)./scale); hold on;
+ s.EdgeColor = 'none';
+ % s.FaceAlpha = 0.5;
+ % xz plane
+ s = surface(toplot_xz.X,OPTIONS.XYZ(2)+0*toplot_xz.X,toplot_xz.Y,toplot_xz.FIELD(:,:,i_)./scale);
+ s.EdgeColor = 'none';
+ % s.FaceAlpha = 0.7;
+ % yz plane
+ s = surface(OPTIONS.XYZ(1)+0*toplot_yz.X,toplot_yz.X,toplot_yz.Y,toplot_yz.FIELD(:,:,i_)./scale);
+ s.EdgeColor = 'none';
+ % s.FaceAlpha = 0.7;
+ xlabel('x');
+ ylabel('y');
+ zlabel('z');
+ h = gca;
+ plot3(h.XLim,OPTIONS.XYZ(2)*[1 1],OPTIONS.XYZ(3)*[1 1],'--k','LineWidth',1.)
+ plot3(OPTIONS.XYZ(1)*[1 1],h.YLim,OPTIONS.XYZ(3)*[1 1],'--k','LineWidth',1.)
+ plot3(OPTIONS.XYZ(1)*[1 1],OPTIONS.XYZ(2)*[1 1],h.ZLim,'--k','LineWidth',1.)
+ % zline(OPTIONS.XYZ(3),'-k','LineWidth',1.5)
+ view([1 -1 0.25])
+ end
else
contour(toplot.X,toplot.Y,toplot.FIELD(:,:,i_)./scale,128);
% pclr = pcolor(toplot.X,toplot.Y,toplot.FIELD(:,:,FRAMES(i_))./scale); set(pclr, 'edgecolor','none'); shading interp
@@ -47,12 +76,23 @@ FIGURE.fig = figure; %set(gcf, 'Position', toplot.DIMENSIONS.*[1 1 Ncols Nrows]
end
xlabel(toplot.XNAME); ylabel(toplot.YNAME);
% if i_ > 1; set(gca,'ytick',[]); end;
- title([sprintf('$t c_s/R=%.0f$',tshot),', max = ',sprintf('%.1e',frame_max)]);
+ title([sprintf('$t c_s/R=%5.2f$',tshot),', max = ',sprintf('%.1e',frame_max)]);
if OPTIONS.CLIMAUTO
clim('auto')
end
end
+ if OPTIONS.AXISEQUAL
+ pbaspect(toplot.ASPECT)
+ end
+ if ~strcmp(OPTIONS.PLAN,'kxky')
+ clim([-1,1]*frame_max/scale);
+ colormap(bluewhitered);
+ if OPTIONS.INTERP
+ shading interp;
+ end
+ end
legend(['$',toplot.FIELDNAME,'$']);
+ legend('Location','northeast');
FIGURE.FIGNAME = [FNAME,'_snaps',TNAME];
end
diff --git a/matlab/plot/plot_ballooning.m b/matlab/plot/plot_ballooning.m
index f43866e82dd9cc6cf5cb02dad988bbd859b9f083..aa6565c71b034dbd4be27b87d4be6e42c20b534d 100644
--- a/matlab/plot/plot_ballooning.m
+++ b/matlab/plot/plot_ballooning.m
@@ -6,12 +6,12 @@ function [FIG, b_angle, phib, psib, ky_] = plot_ballooning(DATA,OPTIONS)
[~,ikyarray] = min(abs(DATA.grids.ky - OPTIONS.kymodes));
ikyarray = unique(ikyarray);
dz = DATA.grids.z(2) - DATA.grids.z(1);
- phi_real=real(DATA.PHI(:,:,:,it1));
- phi_imag=imag(DATA.PHI(:,:,:,it1));
+ phi_real=mean(real(DATA.PHI(:,:,:,it0:it1)),4);
+ phi_imag=mean(imag(DATA.PHI(:,:,:,it0:it1)),4);
ncol = 1;
if DATA.inputs.BETA > 0
- psi_real=real(DATA.PSI(:,:,:,it1));
- psi_imag=imag(DATA.PSI(:,:,:,it1));
+ psi_real=mean(real(DATA.PSI(:,:,:,it0:it1)),4);
+ psi_imag=mean(imag(DATA.PSI(:,:,:,it0:it1)),4);
ncol = 2;
end
if OPTIONS.PLOT_KP
@@ -50,12 +50,12 @@ function [FIG, b_angle, phib, psib, ky_] = plot_ballooning(DATA,OPTIONS)
title(['$k_y=',sprintf('%2.2f',DATA.grids.ky(iky)),...
',t_{avg}\in [',sprintf('%1.1f',DATA.Ts3D(it0)),',',sprintf('%1.1f',DATA.Ts3D(it1)),']$']);
% z domain start end point
- for i_ = 2*[1:DATA.grids.Nkx-1]-(DATA.grids.Nkx+1)
- xline(DATA.grids.Npol*(i_),'HandleVisibility','off'); hold on;
- end
- for i_ = 2*[2:DATA.grids.Nkx]-(DATA.grids.Nkx+1)
- xline(DATA.grids.Npol*(i_)-dz/pi,'HandleVisibility','off'); hold on;
- end
+ % for i_ = 2*[1:DATA.grids.Nkx-1]-(DATA.grids.Nkx+1)
+ % xline(DATA.grids.Npol*(i_),'HandleVisibility','off'); hold on;
+ % end
+ % for i_ = 2*[2:DATA.grids.Nkx]-(DATA.grids.Nkx+1)
+ % xline(DATA.grids.Npol*(i_)-dz/pi,'HandleVisibility','off'); hold on;
+ % end
if DATA.inputs.BETA > 0
[psib_real,b_angle] = ballooning_representation(psi_real,DATA.inputs.SHEAR,DATA.grids.Npol,DATA.grids.kx,iky,DATA.grids.ky,DATA.grids.z);
[psib_imag,~] = ballooning_representation(psi_imag,DATA.inputs.SHEAR,DATA.grids.Npol,DATA.grids.kx,iky,DATA.grids.ky,DATA.grids.z);
@@ -71,12 +71,12 @@ function [FIG, b_angle, phib, psib, ky_] = plot_ballooning(DATA,OPTIONS)
legend('real','imag','norm')
% z domain start end point
- for i_ = 2*[1:DATA.grids.Nkx-1]-(DATA.grids.Nkx+1)
- xline(DATA.grids.Npol*(i_),'HandleVisibility','off'); hold on;
- end
- for i_ = 2*[2:DATA.grids.Nkx]-(DATA.grids.Nkx+1)
- xline(DATA.grids.Npol*(i_)-dz/pi,'HandleVisibility','off'); hold on;
- end
+ % for i_ = 2*[1:DATA.grids.Nkx-1]-(DATA.grids.Nkx+1)
+ % xline(DATA.grids.Npol*(i_),'HandleVisibility','off'); hold on;
+ % end
+ % for i_ = 2*[2:DATA.grids.Nkx]-(DATA.grids.Nkx+1)
+ % xline(DATA.grids.Npol*(i_)-dz/pi,'HandleVisibility','off'); hold on;
+ % end
xlabel('$\chi / \pi$')
ylabel('$\psi_B (\chi)$');
title(['$k_y=',sprintf('%2.2f',DATA.grids.ky(iky)),...
@@ -89,12 +89,12 @@ function [FIG, b_angle, phib, psib, ky_] = plot_ballooning(DATA,OPTIONS)
subplot(numel(ikyarray),ncol,ncol*(iplot-1)+3)
plot(b_angle/pi, kperpb,'-k'); hold on;
% z domain start end point
- for i_ = 2*[1:DATA.grids.Nkx-1]-(DATA.grids.Nkx+1)
- xline(DATA.grids.Npol*(i_),'HandleVisibility','off'); hold on;
- end
- for i_ = 2*[2:DATA.grids.Nkx]-(DATA.grids.Nkx+1)
- xline(DATA.grids.Npol*(i_)-dz/pi,'HandleVisibility','off'); hold on;
- end
+ % for i_ = 2*[1:DATA.grids.Nkx-1]-(DATA.grids.Nkx+1)
+ % xline(DATA.grids.Npol*(i_),'HandleVisibility','off'); hold on;
+ % end
+ % for i_ = 2*[2:DATA.grids.Nkx]-(DATA.grids.Nkx+1)
+ % xline(DATA.grids.Npol*(i_)-dz/pi,'HandleVisibility','off'); hold on;
+ % end
xlabel('$\chi / \pi$')
ylabel('$k_\perp (\chi)$');
title(['$k_y=',sprintf('%2.2f',DATA.grids.ky(iky)),...
diff --git a/matlab/plot/plot_cosol_mat.m b/matlab/plot/plot_cosol_mat.m
index ff31e4d26519c7a4f4f2870216fdd478d01a2698..3954ced9b260fe2f729e172b68e32cddb4c0a6c4 100644
--- a/matlab/plot/plot_cosol_mat.m
+++ b/matlab/plot/plot_cosol_mat.m
@@ -1,97 +1,103 @@
-% MAT = results.iCa;
-% figure
-% suptitle('FCGK,P=6,J=3');
-% subplot(221)
-% imagesc(log(abs(MAT)));
-% title('log abs')
-% subplot(222)
-% imagesc(imag(MAT)); colormap(bluewhitered)
-% title('imag'); colorbar
-% subplot(223)
-% imagesc(imag(MAT)>0);
-% title('imag$>$0');
-% subplot(224)
-% imagesc(imag(MAT)<0);
-% title('imag$<$0');
-%
-%
- %% SGGK
-P_ = 20; J_ = 10;
-mat_file_name = '/home/ahoffman/HeLaZ/iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5';
-SGGK_self = h5read(mat_file_name,'/00000/Caapj/Ciipj');
-SGGK_ei = h5read(mat_file_name,'/00000/Ceipj/CeipjF')+h5read(mat_file_name,'/00000/Ceipj/CeipjT');
-SGGK_ie = h5read(mat_file_name,'/00000/Ciepj/CiepjF')+h5read(mat_file_name,'/00000/Ciepj/CiepjT');
+%% Coeff trajectory
+% matfilename = 'gk_landau_P10_J5_dk_5e-2_km_2.0_NFLR_12.h5'; JMAX = 5;
+% matfilename = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'; JMAX = 10;
+% matfilename = 'gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5'; JMAX = 10;
+matfilename = 'gk_landau_P16_J9_dk_5e-2_km_2.0_NFLR_8.h5'; JMAX = 10;
+P_ = 10; J_ = P_/2;
+n = 1:((P_+1)*(J_+1));
+for p_ = 0:P_
+ n((0:J_)+p_*(J_+1)+1) = (0:J_)+p_*(JMAX+1);
+end
+mat_file_name = ['/home/ahoffman/gyacomo/iCa/',matfilename];
+kp_a = h5read(mat_file_name,'/coordkperp');
+coeff = kp_a*0;
+gmax = kp_a*0;
+p1 = 0; j1 = 0;
+p2 = 0; j2 = 0;
+for ik = 1:numel(kp_a)
+ matidx = sprintf('%5.5i',ik-1);
+ M_self = h5read(mat_file_name,['/',matidx,'/Caapj/Ciipj']);
+ MAT = M_self(n+1,n+1);
+ ipj1 = n(j1+p1*(J_+1)+1);
+ ipj2 = n(j2+p2*(J_+1)+1);
+ coeff(ik) = MAT(ipj1+1,ipj2+1);
+ gmax(ik) = -max((real(eig(MAT))));
+
+end
+figure
+plot(kp_a,gmax);
+ %% SGGK
+P_ = 10; J_ = 5;
+n = 1:((P_+1)*(J_+1));
+for p_ = 0:P_
+ n((0:J_)+p_*(J_+1)+1) = (0:J_)+p_*(10+1);
+end
+kp = 5.0;
+kp_a = h5read(mat_file_name,'/coordkperp');
+[~,matidx] = min(abs(kp_a-kp));
+matidx = sprintf('%5.5i',matidx);
+mat_file_name = '/home/ahoffman/gyacomo/iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5';
+M_self = h5read(mat_file_name,['/',matidx,'/Caapj/Ciipj']);
figure
-MAT = 1i*SGGK_self; suptitle('SGGK ii,P=20,J=10, k=0');
-% MAT = 1i*SGGK_ei; suptitle('SGGK ei,P=20,J=10');
-% MAT = 1i*SGGK_ie; suptitle('SGGK ie,P=20,J=10');
-subplot(221)
- imagesc(abs(MAT));
- title('log abs')
-subplot(222)
- imagesc(imag(MAT)); colormap(bluewhitered)
- title('imag'); colorbar
-subplot(223)
- imagesc(imag(MAT)>0);
- title('imag$>$0');
-subplot(224)
- imagesc(imag(MAT)<0);
- title('imag$<$0');
+MAT = M_self; %suptitle('SGGK ii,P=20,J=10, k=0');
+imagesc((MAT(n+1,n+1)));
+title('Sugama')
+xlim('tight')
+ylim('tight')
+axis equal
+% clim(1*[-1 1])
+clim auto
+colormap(bluewhitered)
%% PAGK
-P_ = 20; J_ = 10;
-mat_file_name = '/home/ahoffman/HeLaZ/iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5';
-PAGK_self = h5read(mat_file_name,'/00000/Caapj/Ceepj');
-% PAGK_ei = h5read(mat_file_name,'/00000/Ceipj/CeipjF')+h5read(mat_file_name,'/00000/Ceipj/CeipjT');
-% PAGK_ie = h5read(mat_file_name,'/00000/Ciepj/CiepjF')+h5read(mat_file_name,'/00000/Ciepj/CiepjT');
-
+P_ = 10; J_ = 5;
+n = 1:((P_+1)*(J_+1));
+for p_ = 0:P_
+ n((0:J_)+p_*(J_+1)+1) = (0:J_)+p_*(10+1);
+end
+kp = 1.0;
+kp_a = h5read(mat_file_name,'/coordkperp');
+[~,matidx] = min(abs(kp_a-kp));
+matidx = sprintf('%5.5i',matidx);
+mat_file_name = '/home/ahoffman/gyacomo/iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5';
+PAGK_self = h5read(mat_file_name,['/',matidx,'/Caapj/Ciipj']);
figure
-MAT = 1i*PAGK_self; suptitle('PAGK ii,P=20,J=10, k=0');
-
-subplot(221)
- imagesc(abs(MAT));
- title('log abs')
-subplot(222)
- imagesc(imag(MAT)); colormap(bluewhitered)
- title('imag'); colorbar
-subplot(223)
- imagesc(imag(MAT)>0);
- title('imag$>$0');
-subplot(224)
- imagesc(imag(MAT)<0);
- title('imag$<$0');
+PA_MAT = PAGK_self;
+ imagesc((PA_MAT(n+1,n+1)));
+ title('Lorentz')
+ xlim('tight')
+ ylim('tight')
+ axis equal
+ % clim(1*[-1 1])
+ clim auto
+ colormap(bluewhitered)
%% FCGK
-P_ = 4; J_ = 2;
-mat_file_name = '/home/ahoffman/gyacomo/iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5';
-% mat_file_name = '/home/ahoffman/gyacomo/iCa/LDGK_P6_J3_dk_5e-2_km_2.5_NFLR_20.h5';
-% mat_file_name = '/home/ahoffman/gyacomo/iCa/LDGK_P10_J5_dk_5e-2_km_5_NFLR_12.h5';
-
-kp = 0.0;
+P_ = 10; J_ = 5;
+n = 1:((P_+1)*(J_+1));
+for p_ = 0:P_
+ n((0:J_)+p_*(J_+1)+1) = (0:J_)+p_*(5+1);
+end
+mat_file_name = '/home/ahoffman/gyacomo/iCa/gk_landau_P10_J5_dk_5e-2_km_2.0_NFLR_12.h5';
+kp = 1.0;
kp_a = h5read(mat_file_name,'/coordkperp');
[~,matidx] = min(abs(kp_a-kp));
matidx = sprintf('%5.5i',matidx);
FCGK_self = h5read(mat_file_name,['/',matidx,'/Caapj/Ciipj']);
-FCGK_ei = h5read(mat_file_name,['/',matidx,'/Ceipj/CeipjF'])+h5read(mat_file_name,'/00000/Ceipj/CeipjT');
-FCGK_ie = h5read(mat_file_name,['/',matidx,'/Ciepj/CiepjF'])+h5read(mat_file_name,'/00000/Ciepj/CiepjT');
figure
-MAT = 1i*FCGK_self; suptitle(['FCGK ii,P=6,J=3, k=',num2str(kp),' (',matidx,')']);
+FC_MAT = FCGK_self;
% MAT = 1i*FCGK_ei; suptitle('FCGK ei,P=20,J=10');
% MAT = 1i*FCGK_ie; suptitle('FCGK ie,P=20,J=10');
-subplot(221)
- imagesc(abs(MAT));
- title('log abs')
-subplot(222)
- imagesc(imag(MAT)); colormap(bluewhitered)
- title('imag'); colorbar
-subplot(223)
- imagesc(imag(MAT)>0);
- title('imag$>$0');
-subplot(224)
- imagesc(imag(MAT)<0);
- title('imag$<$0');
+ imagesc((FC_MAT(n+1,n+1)));
+ title('Landau')
+ xlim('tight')
+ ylim('tight')
+ axis equal
+ % clim(1*[-1 1])
+ clim auto
+ colormap(bluewhitered)
%% Eigenvalue spectrum analysis
if 0
@@ -129,16 +135,19 @@ for j_ = 1:numel(mfns)
% kp_a = kp_a(kp_a<=3);
gmax = zeros(size(kp_a));
wmax = zeros(size(kp_a));
+ c44 = zeros(size(kp_a));
for idx_ = 0:numel(kp_a)-1
matidx = sprintf('%5.5i',idx_);
- MAT = h5read(mat_file_name,['/',matidx,'/Caapj/Ciipj']);
- gmax(idx_+1) = max((real(eig(MAT))));
- wmax(idx_+1) = max((imag(eig(MAT))));
+ FC_MAT = h5read(mat_file_name,['/',matidx,'/Caapj/Ciipj']);
+ gmax(idx_+1) = max((real(eig(FC_MAT))));
+ wmax(idx_+1) = max((imag(eig(FC_MAT))));
+ c44 (idx_+1) = FC_MAT(4,4);
end
subplot(121)
plot(kp_a,gmax,'DisplayName',CONAME_A{j_}); hold on;
subplot(122)
- plot(kp_a,wmax,'DisplayName',CONAME_A{j_}); hold on;
+ % plot(kp_a,wmax,'DisplayName',CONAME_A{j_}); hold on;
+ plot(kp_a,c44,'DisplayName',CONAME_A{j_}); hold on;
end
subplot(121)
legend('show'); grid on;
@@ -195,9 +204,9 @@ for j_ = 1:numel(mfns)
kp_a = h5read(mat_file_name,'/coordkperp');
[~,idx_] = min(abs(kp_a-kperp));
matidx = sprintf('%5.5i',idx_);
- MAT = h5read(mat_file_name,['/',matidx,'/Caapj/Ciipj']);
- grow{j_} = real(eig(MAT))*TAU_A(j_);
- puls{j_} = imag(eig(MAT))*TAU_A(j_);
+ FC_MAT = h5read(mat_file_name,['/',matidx,'/Caapj/Ciipj']);
+ grow{j_} = real(eig(FC_MAT))*TAU_A(j_);
+ puls{j_} = imag(eig(FC_MAT))*TAU_A(j_);
end
figure
diff --git a/matlab/plot/plot_isosurf_wip.m b/matlab/plot/plot_isosurf_wip.m
new file mode 100644
index 0000000000000000000000000000000000000000..886c29fa643311a250986760a1197105ab3b7de7
--- /dev/null
+++ b/matlab/plot/plot_isosurf_wip.m
@@ -0,0 +1,20 @@
+options.NAME = ['N_i^{00}'];
+options.PLAN = 'xyz';
+options.TIME = [30];
+TOPLOT = process_field( data, options );
+
+[X,Y,Z] = meshgrid(data.grids.x,data.grids.y,data.grids.z);
+X = smooth3(X,'gaussian',5);
+Y = smooth3(Y,'gaussian',5);
+Z = smooth3(Z,'gaussian',5);
+TOPLOT.FIELD = smooth3(TOPLOT.FIELD,'gaussian',5);
+s = isosurface(X,Y,Z,TOPLOT.FIELD,0);
+s = reducepatch(s,0.1);
+%% figure
+
+figure
+
+p=patch(s);
+set(p,'EdgeColor','none');
+set(p,'FaceColor','r');
+set(p,'FaceAlpha',0.4);
diff --git a/matlab/plot/plot_kernels.m b/matlab/plot/plot_kernels.m
index a02b82cae6367eb22d4012c19859286eec4407cb..2e7f700e9878c13d7766329bfc8f7bf4f915d86b 100644
--- a/matlab/plot/plot_kernels.m
+++ b/matlab/plot/plot_kernels.m
@@ -2,15 +2,15 @@ if 0
%% Kernels
kmax=5;
nmax=6;
-kx = linspace(0,kmax,100);
+kp = linspace(0,kmax,100);
figure
for n_ = 0:nmax
- plot(kx,kernel(n_,kx),'DisplayName',['$\mathcal{K}_{',num2str(n_),'}$']);hold on;
+ plot(kp,kernel(n_,kp),'DisplayName',['$\mathcal{K}_{',num2str(n_),'}$']);hold on;
end
ylim_ = ylim;
-plot(kx(end)*[2/3 2/3],ylim_,'--k','DisplayName','AA');
+plot(kp(end)*[2/3 2/3],ylim_,'--k','DisplayName','AA');
% J0 = besselj(0,kx);
% plot(kx,J0,'-r','DisplayName','$J_0$');
legend('show'); xlabel('k'); title('Kernel functions $\mathcal{K}_n$');
@@ -19,32 +19,32 @@ end
if 0
%% Bessels and approx
-vperp = linspace(0,1.5,4);
+wperp = linspace(0,1.5,4);
nmax1=5;
nmax2=10;
kmax=7;
figure
-for i = 1:4
-subplot(2,2,i)
- v_ = vperp(i);
- kx = linspace(0,kmax,100);
+for id2 = 1:4
+subplot(2,2,id2)
+ v_ = wperp(id2);
+ kp = linspace(0,kmax,100);
- J0 = besselj(0,kx*v_);
- A1 = 1 - kx.^2*v_^2/4;
- K1 = zeros(size(kx));
- K2 = zeros(size(kx));
+ J0 = besselj(0,kp*v_);
+ A1 = 1 - kp.^2*v_^2/4;
+ K1 = zeros(size(kp));
+ K2 = zeros(size(kp));
for n_ = 0:nmax1
- K1 = K1 + kernel(n_,kx).*polyval(LaguerrePoly(n_),v_^2);
+ K1 = K1 + kernel(n_,kp).*polyval(LaguerrePoly(n_),v_^2);
end
for n_ = 0:nmax2
- K2 = K2 + kernel(n_,kx).*polyval(LaguerrePoly(n_),v_^2);
+ K2 = K2 + kernel(n_,kp).*polyval(LaguerrePoly(n_),v_^2);
end
- plot(kx,J0,'-k','DisplayName','$J_0(kv)$'); hold on;
- plot(kx,A1,'-r','DisplayName','$1 - k^2 v^2/4$');
- plot(kx,K1,'--b','DisplayName',['$\sum_{n=0}^{',num2str(nmax1),'}\mathcal K_n(k) L^n(v)$']);
- plot(kx,K2,'-b','DisplayName',['$\sum_{n=0}^{',num2str(nmax2),'}\mathcal K_n(k) L^n(v)$']);
+ plot(kp,J0,'-k','DisplayName','$J_0(kv)$'); hold on;
+ plot(kp,A1,'-r','DisplayName','$1 - k^2 v^2/4$');
+ plot(kp,K1,'--b','DisplayName',['$\sum_{n=0}^{',num2str(nmax1),'}\mathcal K_n(k) L^n(v)$']);
+ plot(kp,K2,'-b','DisplayName',['$\sum_{n=0}^{',num2str(nmax2),'}\mathcal K_n(k) L^n(v)$']);
ylim_ = [-0.5, 1.0];
- plot(kx(end)*[2/3 2/3],ylim_,'--k','DisplayName','AA');
+ plot(kp(end)*[2/3 2/3],ylim_,'--k','DisplayName','AA');
ylim(ylim_); xlabel('$k$')
legend('show'); grid on; title(['$v = ',num2str(v_),'$'])
end
@@ -53,165 +53,223 @@ end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if 0
%% from Sum_n Kn x Ln x Lj to Sum_n Kn Sum_s dnjs x Ls
-vperp = [2];
-kx = linspace(0,10,200);
-Jmax = 15;
-j = 10;
-fig = figure; set(gcf, 'Position', [100, 100, numel(vperp)*300, 250])
-% suptitle(['$J_{max}=',num2str(Jmax),', j=',num2str(j),'$'])
+wperp = [0 0.2 0.4];
+kp = linspace(0,3,50);
+Jmax = 8;
+js = 4;
+fig = figure; set(gcf, 'Position', [100, 100, numel(wperp)*300, 250])
+% suptitle(['$J_{ma}=',num2str(Jmax),', j=',num2str(j),'$'])
Kn = @(x__,y__) kernel(x__,y__);
-for i = 1:numel(vperp)
-subplot(1,numel(vperp),i)
- v_ = vperp(i);
+% dnjs_ = @(n,j,s) dnjs(n,j,s);
+dnjs_ = @(n,j,s) dnjs_GYAC(n+1,j+1,s+1);
+for id2 = 1:numel(wperp)
+subplot(1,numel(wperp),id2)
+ v_ = wperp(id2);
% J0 x Lj
- J0Lj = besselj(0,kx*v_)*polyval(LaguerrePoly(j),v_^2);
+ J0Lj_ = besselj(0,kp*v_)*polyval(LaguerrePoly(js),v_^2);
% Taylor exp of J0
- A1 = (1 - kx.^2*v_^2/4)*polyval(LaguerrePoly(j),v_^2);
+ A1 = (1 - kp.^2*v_^2/4)*polyval(LaguerrePoly(js),v_^2);
% Sum_n^Nmax Kn x Ln x Lj
- KLL = zeros(size(kx));
+ KLL = zeros(size(kp));
for n_ = 0:Jmax
- KLL = KLL + Kn(n_,kx).*polyval(LaguerrePoly(n_),v_^2);
+ KLL = KLL + Kn(n_,kp).*polyval(LaguerrePoly(n_),v_^2);
end
- KLL = KLL.*polyval(LaguerrePoly(j),v_^2);
+ KLL = KLL.*polyval(LaguerrePoly(js),v_^2);
% Sum_n^Nmax Kn Sum_s^Smax dnjs x Ls
- KdL1 = zeros(size(kx));
- for n_ = 0:Jmax-j
+ KdL1_ = zeros(size(kp));
+ for n_ = 0:Jmax-js
tmp_ = 0;
- for s_ = 0:n_+j
- tmp_ = tmp_ + dnjs(n_,j,s_)*polyval(LaguerrePoly(s_),v_^2);
+ for s_ = 0:n_+js
+ tmp_ = tmp_ + dnjs_(n_,js,s_)*polyval(LaguerrePoly(s_),v_^2);
end
- KdL1 = KdL1 + Kn(n_,kx).*tmp_;
+ KdL1_ = KdL1_ + Kn(n_,kp).*tmp_;
end
% Sum_n^Nmax Kn Sum_s^Smax dnjs x Ls
- KdL2 = zeros(size(kx));
+ KdL2_ = zeros(size(kp));
for n_ = 0:Jmax
tmp_ = 0;
- for s_ = 0:min(Jmax,n_+j)
- tmp_ = tmp_ + dnjs(n_,j,s_)*polyval(LaguerrePoly(s_),v_^2);
+ for s_ = 0:min(Jmax,n_+js)
+ tmp_ = tmp_ + dnjs_(n_,js,s_)*polyval(LaguerrePoly(s_),v_^2);
end
- KdL2 = KdL2 + Kn(n_,kx).*tmp_;
+ KdL2_ = KdL2_ + Kn(n_,kp).*tmp_;
end
% plots
- plot(kx,J0Lj,'-k','DisplayName','$J_0 L_j$'); hold on;
- plot(kx,KLL,'-','DisplayName',['$\sum_{n=0}^{J}\mathcal K_n L^n L_j$']);
- plot(kx,KdL1,'-.','DisplayName',['$\sum_{n=0}^{J-j}\mathcal K_n \sum_{s=0}^{n+j} d_{njs} L^s$']);
+ plot(kp,J0Lj_,'-k','DisplayName','$J_0 L_j$');
ylim_ = ylim;
- plot(kx,A1,'-','DisplayName','$(1 - k^2 v^2/4) L_j$'); hold on;
- plot(kx,KdL2,'--','DisplayName',['$\sum_{n=0}^{J}\mathcal K_n \sum_{s=0}^{\min(J,n+j)} d_{njs} L^s$']);
+ plot(kp,KLL,'-','DisplayName',['$\sum_{n=0}^{J}\mathcal K_n L_n L_j$']); hold on;
+ plot(kp,KdL2_,':o','DisplayName',['$\sum_{n=0}^{J}\mathcal K_n \sum_{s=0}^{\min(J,n+j)} d_{njs} L_s$']);
+ plot(kp,J0Lj_,'-k','DisplayName','$J_0 L_j$'); hold on;
+ ylim_ = ylim;
+ plot(kp,A1,':s','DisplayName','$(1 - l_\perp) L_j$'); hold on;
+ plot(kp,KdL1_,':x','MarkerSize',10,...
+ 'DisplayName',['$\sum_{n=0}^{J-j}\mathcal K_n \sum_{s=0}^{n+j} d_{njs} L_s$']);
+ plot(kp,J0Lj_,'-k','DisplayName','$J_0 L_j$'); hold on;
% plot(kx(end)*[2/3 2/3],ylim_,'--k','DisplayName','AA');
ylim(ylim_);
- xlabel('$k_{\perp}$')
- legend('show');
- grid on; title(['$v = ',num2str(v_),'$',' $J_{max}=',num2str(Jmax),', j=',num2str(j),'$'])
+ xlabel('$k_{\perp}\rho_s$')
+ % legend('show');
+ grid on; title(['$w_\perp = ',num2str(v_),'$',' $J_{max}=',num2str(Jmax),', j=',num2str(js),'$'])
end
-FIGNAME = ['/home/ahoffman/Dropbox/Applications/Overleaf/Hermite-Laguerre Z-pinch (HeLaZ Doc.)/figures/J0Lj_approx_J','_',num2str(Jmax),'_j_',num2str(j),'.eps'];
+FIGNAME = ['/home/ahoffman/Dropbox/Applications/Overleaf/Hermite-Laguerre Z-pinch (HeLaZ Doc.)/figures/J0Lj_approx_J','_',num2str(Jmax),'_j_',num2str(js),'.eps'];
% saveas(fig,FIGNAME,'epsc');
end
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-if 1
-%% Convergence analysis
-kx = linspace(0,10,200);
-v_A = linspace(0,1,10);
-J_A = 1:15;
-[YY,XX] = meshgrid(v_A,J_A);
-klimLL = zeros(size(XX));
-klimdL1= zeros(size(XX));
-Kn = @(x__,y__) kernel(x__,y__);
-fig = figure; set(gcf, 'Position', [100, 100, 500, 400]);
-for j = 5
-for ij_ = 1:numel(J_A)
- iv_ = 1;
- for v_ = v_A
- eps = abs(0.01*polyval(LaguerrePoly(j),v_^2));
- Jmax = J_A(ij_);
- % J0 x Lj
- J0Lj = besselj(0,kx*v_)*polyval(LaguerrePoly(j),v_^2);
- % Taylor exp of J0
- A1 = (1 - kx.^2*v_^2/4)*polyval(LaguerrePoly(j),v_^2);
- % Sum_n^Nmax Kn x Ln x Lj
- KLL = zeros(size(kx));
- for n_ = 0:Jmax
- KLL = KLL + Kn(n_,kx).*polyval(LaguerrePoly(n_),v_^2);
- end
- KLL = KLL.*polyval(LaguerrePoly(j),v_^2);
- % Sum_n^Nmax Kn Sum_s^Smax dnjs x Ls
- KdL1 = zeros(size(kx));
- for n_ = 0:Jmax-j
- tmp_ = 0;
- for s_ = 0:n_+j
- tmp_ = tmp_ + dnjs(n_,j,s_)*polyval(LaguerrePoly(s_),v_^2);
- end
- KdL1 = KdL1 + Kn(n_,kx).*tmp_;
- end
- % Sum_n^Nmax Kn Sum_s^Smax dnjs x Ls
- KdL2 = zeros(size(kx));
- for n_ = 0:Jmax
- tmp_ = 0;
- for s_ = 0:min(Jmax,n_+j)
- tmp_ = tmp_ + dnjs(n_,j,s_)*polyval(LaguerrePoly(s_),v_^2);
- end
- KdL2 = KdL2 + Kn(n_,kx).*tmp_;
- end
- % errors
- err_ = abs(J0Lj-KLL);
- idx_ = 1;
- while(err_(idx_)< eps && idx_ < numel(err_))
- idx_ = idx_ + 1;
- end
- klimLL(ij_,iv_) = kx(idx_);
- err_ = abs(J0Lj-KdL1);
- idx_ = 1;
- while(err_(idx_)< eps && idx_ < numel(err_))
- idx_ = idx_ + 1;
- end
- klimdL1(ij_,iv_) = kx(idx_);
- iv_ = iv_ + 1;
- end
-end
-% plot
-% plot(JmaxA,klimLL,'o-'); hold on
-% plot(J_A,klimdL1,'o-'); hold on
-end
-surf(XX,YY,klimdL1)
-xlabel('$J_{max}$'); ylabel('$v_{\perp}$');
-title(['$k$ s.t. $\epsilon=1\%$, $j=',num2str(j),'$'])
-ylim([0,1]); grid on;
-end
if 0
%% Test Lj Ln = sum_s dnjs Ls
-j = 10;
-n = 10;
-smax = n+j;
-vperp = linspace(0,5,20);
-LjLn = zeros(size(vperp));
-dnjsLs = zeros(size(vperp));
-dnjsLs_bin = zeros(size(vperp));
-dnjsLs_stir = zeros(size(vperp));
-
-i=1;
-for x_ = vperp
- LjLn(i) = polyval(LaguerrePoly(j),x_)*polyval(LaguerrePoly(n),x_);
- dnjsLs(i) = 0;
- dnjsLs_bin(i) = 0;
- dnjsLs_stir(i) = 0;
+js = 4;
+n = 4;
+GYAC = 1;
+smax = n+js;
+wperp = linspace(0,5,50);
+LjLn = zeros(size(wperp));
+dnjsLs = zeros(size(wperp));
+dnjsLs_bin = zeros(size(wperp));
+dnjsLs_stir = zeros(size(wperp));
+if GYAC
+ dnjsLs_GYAC = zeros(size(wperp));
+ % Specify the filename
+ filename = '/home/ahoffman/gyacomo/results/dbg/output.txt';
+ % Use the load function to read data from the text file
+ data = load(filename);
+ % Extract columns
+ indices = data(:, 1:3);
+ values = data(:, 4);
+ is_ = unique(indices(:,1));
+ N_ = numel(is_);
+ dnjs_GYAC = zeros(N_,N_,N_);
+ for id2 = 1:numel(values)
+ in_ = indices(id2,1);
+ ij_ = indices(id2,2);
+ is_ = indices(id2,3);
+ dnjs_GYAC(in_,ij_,is_) = values(id2);
+ end
+end
+
+id2=1;
+for x_ = wperp
+ LjLn(id2) = polyval(LaguerrePoly(js),x_)*polyval(LaguerrePoly(n),x_);
+ dnjsLs(id2) = 0;
+ dnjsLs_bin(id2) = 0;
+ dnjsLs_stir(id2) = 0;
for s_ = 0:smax
- dnjsLs(i) = dnjsLs(i) + dnjs(n,j,s_)*polyval(LaguerrePoly(s_),x_);
- dnjsLs_bin(i) = dnjsLs_bin(i) + dnjs_fact(n,j,s_)*polyval(LaguerrePoly(s_),x_);
- dnjsLs_stir(i) = dnjsLs_stir(i) + dnjs_stir(n,j,s_)*polyval(LaguerrePoly(s_),x_);
+ dnjsLs(id2) = dnjsLs(id2) + dnjs(n,js,s_)*polyval(LaguerrePoly(s_),x_);
+ dnjsLs_bin(id2) = dnjsLs_bin(id2) + dnjs_fact(n,js,s_)*polyval(LaguerrePoly(s_),x_);
+ dnjsLs_stir(id2) = dnjsLs_stir(id2) + dnjs_stir(n,js,s_)*polyval(LaguerrePoly(s_),x_);
+ if GYAC
+ dnjsLs_GYAC(id2) = dnjsLs_GYAC(id2) + dnjs_GYAC(n+1,js+1,s_+1)*polyval(LaguerrePoly(s_),x_);
+ end
end
- i = i+1;
+ id2 = id2+1;
end
figure
-plot(vperp,LjLn); hold on;
-plot(vperp,dnjsLs,'d')
-plot(vperp,dnjsLs_bin,'o')
-plot(vperp,dnjsLs_stir/dnjsLs_stir(1),'x')
+plot(wperp,LjLn,'DisplayName',['$L_',num2str(js),' L_',num2str(n),'$']); hold on;
+plot(wperp,dnjsLs,'d','DisplayName',...
+ ['$\sum_{s=0}^{',num2str(smax),'} d^{coeff}_{',num2str(n),',',num2str(js),',s}L_s$'])
+plot(wperp,dnjsLs_bin,'o','DisplayName','$\sum_s d^{fact}_{njs}L_s$')
+plot(wperp,dnjsLs_stir/dnjsLs_stir(1),'x','DisplayName','$\sum_s d^{stir}_{njs}L_s$')
+plot(wperp,dnjsLs_GYAC,'x','DisplayName','$GYAC$')
+xlabel('$w_\perp$');
% We see that starting arround j = 18, n = 0, the stirling formula is the only
% approximation that does not diverge from the target function. However it
% performs badly for non 0 n...
end
-
+if 0
+%%
+ wperp = linspace(0,3,128);
+ kp = linspace(0,3,128);
+ Jmax = 8;
+ js = Jmax/2;[0 2:Jmax];
+ err_KLL = 1:numel(js);
+ err_KdL1 = 1:numel(js);
+ err_KdL2 = 1:numel(js);
+ err_T1 = 1:numel(js);
+ J0Lj = zeros(numel(kp),numel(wperp));
+ KdL1 = zeros(numel(kp),numel(wperp));
+ KdL2 = zeros(numel(kp),numel(wperp));
+ T1 = zeros(numel(kp),numel(wperp));
+ Kn = @(x__,y__) kernel(x__,y__);
+ % dnjs_ = @(n,j,s) dnjs(n,j,s);
+ dnjs_ = @(n,j,s) dnjs_GYAC(n+1,j+1,s+1);
+ for id1 = 1:numel(js)
+ for id2 = 1:numel(wperp)
+ j_ = js(id1);
+ v_ = wperp(id2);
+ % J0 x Lj
+ J0Lj_ = besselj(0,kp*v_)*polyval(LaguerrePoly(j_),v_^2);
+ % Taylor exp of J0
+ T1_ = (1 - kp.^2*v_^2/4)*polyval(LaguerrePoly(j_),v_^2);
+ % Sum_n^Nmax Kn x Ln x Lj
+ KLL = zeros(size(kp));
+ for n_ = 0:Jmax
+ KLL = KLL + Kn(n_,kp).*polyval(LaguerrePoly(n_),v_^2);
+ end
+ KLL = KLL.*polyval(LaguerrePoly(j_),v_^2);
+ % Sum_n^Nmax Kn Sum_s^Smax dnjs x Ls
+ KdL1_ = zeros(size(kp));
+ for n_ = 0:Jmax-j_
+ tmp_ = 0;
+ for s_ = 0:n_+j_
+ tmp_ = tmp_ + dnjs_(n_,j_,s_)*polyval(LaguerrePoly(s_),v_^2);
+ end
+ KdL1_ = KdL1_ + Kn(n_,kp).*tmp_;
+ end
+ % Sum_n^Nmax Kn Sum_s^Smax dnjs x Ls
+ KdL2_ = zeros(size(kp));
+ for n_ = 0:Jmax
+ tmp_ = 0;
+ for s_ = 0:min(Jmax,n_+j_)
+ tmp_ = tmp_ + dnjs_(n_,j_,s_)*polyval(LaguerrePoly(s_),v_^2);
+ end
+ KdL2_ = KdL2_ + Kn(n_,kp).*tmp_;
+ end
+ f_ = @(x) max(x);
+ kp2 = kp.^2+0.001; intJ0Lj = f_(abs(J0Lj_));
+ err_KLL_(id2) = f_(abs(KLL -J0Lj_)./kp2)./f_(intJ0Lj./kp2);
+ err_KdL1_(id2) = f_(abs(KdL1_-J0Lj_)./kp2)./f_(intJ0Lj./kp2);
+ err_KdL2_(id2) = f_(abs(KdL2_-J0Lj_)./kp2)./f_(intJ0Lj./kp2);
+ err_T1_(id2) = f_(abs(T1_ -J0Lj_)./kp2)./f_(intJ0Lj./kp2);
+ J0Lj(:,id2) = J0Lj_;
+ KdL1(:,id2) = KdL1_;
+ KdL2(:,id2) = KdL2_;
+ T1 (:,id2) = T1_;
+ end
+ err_KLL(id1) = sum(err_KLL_)./numel(wperp);
+ err_KdL1(id1) = sum(err_KdL1_)./numel(wperp);
+ err_KdL2(id1) = sum(err_KdL2_)./numel(wperp);
+ err_T1(id1) = sum(err_T1_)./numel(wperp);
+ end
+ clr_ = lines(10);
+ figure
+ plot(js,err_KLL,'-',...
+ 'DisplayName','$\sum_{n=0}^{J}\mathcal K_n L_n L_j$'); hold on
+ plot(js,err_KdL1,':x',...
+ 'DisplayName','$\sum_{n=0}^{J-j}\mathcal K_n \sum_{s=0}^{n+j} d_{njs} L_s$',...
+ 'Color',clr_(5,:)); hold on
+ plot(js,err_KdL2,':o',...
+ 'DisplayName','$\sum_{n=0}^{J}\mathcal K_n \sum_{s=0}^{\min(J,n+j)} d_{njs} L_s$',...
+ 'Color',clr_(2,:)); hold on
+ plot(js,err_T1,':s',...
+ 'DisplayName','$(1 - l_\perp) L_j$',...
+ 'Color',clr_(4,:)); hold on
+ set(gca,'YScale','log');
+ xlabel('$j$'); ylabel('$\varepsilon_J$')
+ if 1
+ %%
+ figure
+ nc = 25
+ subplot(131)
+ contourf(kp,wperp,J0Lj',nc); xlabel('$l_{\perp a}$'); ylabel('$w_{\perp a}$')
+ clim_ = clim;
+ subplot(132)
+ contourf(kp,wperp,KdL1', nc); xlabel('$l_{\perp a}$'); ylabel('$w_{\perp a}$')
+ clim(clim_);
+ % subplot(132)
+ subplot(133)
+ contourf(kp,wperp,KdL2', nc); xlabel('$l_{\perp a}$'); ylabel('$w_{\perp a}$')
+ % contourf(kp,wperp,T1 ,20); xlabel('$k_\perp$'); ylabel('$w_\perp$')
+ clim(clim_);
+ end
+end
diff --git a/matlab/plot/plot_metric.m b/matlab/plot/plot_metric.m
index 1267999b25d0121511e955858f9b6b53c6509e4f..9a88d3b5dae77018ff8478197e9bd7e90b43ef9b 100644
--- a/matlab/plot/plot_metric.m
+++ b/matlab/plot/plot_metric.m
@@ -96,9 +96,14 @@ if NPLOT > 0
end
end
%outputs
-arrays = squeeze(geo_arrays(:,:));
+% arrays = squeeze(geo_arrays(:,:));
R = [geo_arrays(:,12);geo_arrays(1,12)];
Z = [geo_arrays(:,13);geo_arrays(1,13)];
+
+for i_ = 1:numel(names)
+ namae = names{i_};
+ arrays.(namae) = geo_arrays(:,i_);
+end
try
if ~options.SHOWFIG
close
diff --git a/matlab/plot/plot_radial_transport_and_spacetime.m b/matlab/plot/plot_radial_transport_and_spacetime.m
index 8fdcbc3d2cb4250cb560cd2de67a429ea9032e70..a33d3612a46247dd21d15ac9e292f418a18fbbe8 100644
--- a/matlab/plot/plot_radial_transport_and_spacetime.m
+++ b/matlab/plot/plot_radial_transport_and_spacetime.m
@@ -56,7 +56,8 @@ 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)]});
-%% radial shear radial profile
+ axis tight
+ %% radial shear radial profile
% computation
Ns3D = numel(DATA.Ts3D);
[KX, KY] = meshgrid(DATA.grids.kx, DATA.grids.ky);
@@ -68,18 +69,33 @@ mvm = @(x) movmean(x,OPTIONS.NMVA);
OPTIONS.NAME = OPTIONS.ST_FIELD;
OPTIONS.PLAN = 'xy';
OPTIONS.COMP = 'avg';
+ OPTIONS.TAVG = 0;
OPTIONS.TIME = DATA.Ts3D;
OPTIONS.POLARPLOT = 0;
- toplot = process_field(DATA,OPTIONS);
- f2plot = toplot.FIELD;
+ try
+ if DATA.fort_00.DIAGNOSTICS.dtsave_2d > 0
+ opt_.PLAY = 0;
+ [MOVIE] = cinecita_2D(DATA.dir,OPTIONS.ST_FIELD,opt_);
+ f2plot = MOVIE.F;
+ TIME = MOVIE.T;
+ else
+ toplot = process_field(DATA,OPTIONS);
+ f2plot = toplot.FIELD;
+ TIME = DATA.Ts3D(toplot.FRAMES);
+ end
+ catch
+ toplot = process_field(DATA,OPTIONS);
+ f2plot = toplot.FIELD;
+ TIME = DATA.Ts3D(toplot.FRAMES);
+ end
dframe = ite3D - its3D;
- clim = max(max(max(abs(plt(f2plot(:,:,:))))));
+ clim_ = max(max(max(abs(plt(f2plot(:,:,:))))));
FIGURE.ax3 = subplot(2,1,2,'parent',FIGURE.fig);
- [TY,TX] = meshgrid(DATA.grids.x,DATA.Ts3D(toplot.FRAMES));
+ [TY,TX] = meshgrid(DATA.grids.x,TIME);
pclr = pcolor(TX,TY,squeeze(plt(f2plot))');
set(pclr, 'edgecolor','none');
legend(['$\langle ',OPTIONS.ST_FIELD,' \rangle_{y,z}$']) %colorbar;
- caxis(clim*[-1 1]);
+ clim(clim_*[-1 1]);
cmap = bluewhitered(256);
colormap(cmap)
xlabel('$t c_s/R$'), ylabel('$x/\rho_s$');
diff --git a/matlab/plot/show_moments_spectrum.m b/matlab/plot/show_moments_spectrum.m
index ab9d3aea55bd247e410291c6560aa38f823805d4..f193f118b32725c1187cee06d4fe97f8b1bf215e 100644
--- a/matlab/plot/show_moments_spectrum.m
+++ b/matlab/plot/show_moments_spectrum.m
@@ -6,16 +6,8 @@ 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 == 0
- OPTIONS.LOGSCALE = 0;
-end
-if OPTIONS.LOGSCALE
- logname = 'log';
- compress = @(x,ia) log(sum(abs((x(ia,:,:,:,:))),4));
-else
- logname = '';
- compress = @(x,ia) sum(abs((x(ia,:,:,:,:))),4);
-end
+logname = '';
+compress = @(x,ia) sum(abs((x(ia,:,:,:,:))),4);
for ia = 1:DATA.inputs.Na
Napjz = compress(DATA.Napjz,ia);
Napjz = reshape(Napjz,DATA.grids.Np,DATA.grids.Nj,numel(DATA.Ts3D));
@@ -71,8 +63,10 @@ for ia = 1:DATA.inputs.Na
% plots
% scaling
if OPTIONS.TAVG_2D
- nt_half = ceil(numel(Time_)/2);
- Napjz_tavg = mean(Napjz(:,:,nt_half:end),3);
+ t0 = OPTIONS.TWINDOW(1); t1 = OPTIONS.TWINDOW(2);
+ [~,it0] = min(abs(DATA.Ts3D-t0));
+ [~,it1] = min(abs(DATA.Ts3D-t1));
+ Napjz_tavg = mean(Napjz(:,:,it0:it1),3);
x_ = DATA.grids.Parray;
y_ = DATA.grids.Jarray;
if OPTIONS.TAVG_2D_CTR
@@ -84,7 +78,7 @@ for ia = 1:DATA.inputs.Na
xlabel('$p$'); ylabel('$j$');
end
clb = colorbar; colormap(jet);
- clb.Label.String = '$\log\langle E(p,j) \rangle_t$';
+ clb.Label.String = '$\langle E(p,j) \rangle_t$';
clb.TickLabelInterpreter = 'latex';
clb.Label.Interpreter = 'latex';
clb.Label.FontSize= 18;
@@ -95,7 +89,7 @@ for ia = 1:DATA.inputs.Na
plt = @(x,i) reshape(x(i,:),numel(p2j),numel(Time_));
end
Nlabelmax = 15;
- nskip = floor(numel(p2j)/Nlabelmax);
+ nskip = max(1,floor(numel(p2j)/Nlabelmax));
if OPTIONS.ST
imagesc(Time_,p2j,plt(Na_ST,1:numel(p2j)));
set(gca,'YDir','normal')
@@ -105,10 +99,13 @@ for ia = 1:DATA.inputs.Na
yticklabels(ticks_labels(1:nskip:end))
end
if OPTIONS.FILTER
- caxis([0 0.2]);
- title('Relative fluctuation from the average');
+ clim([0 0.2]);
+ title('Relative fluctuation from the average');
end
colorbar
+ if OPTIONS.LOGSCALE
+ set(gca,'ColorScale','log');
+ end
else
colors_ = jet(numel(p2j));
for i = 1:numel(p2j)
@@ -120,6 +117,7 @@ for ia = 1:DATA.inputs.Na
end
end
top_title(DATA.paramshort)
+ title(plotname)
end
diff --git a/matlab/plot/tight_layout_generator.m b/matlab/plot/tight_layout_generator.m
new file mode 100644
index 0000000000000000000000000000000000000000..f85c09ea5da8cb84dd92c719614570d0a22ecdf7
--- /dev/null
+++ b/matlab/plot/tight_layout_generator.m
@@ -0,0 +1,11 @@
+function [] = tight_layout_generator(m,n)
+% Generate a layout with tight subplots
+% m = 3; % number lines
+% n = 3; % number of columns
+figure
+t = tiledlayout(m,n,'TileSpacing','Compact','Padding','Compact');
+for i = 1:(m*n)
+nexttile
+end
+
+end
diff --git a/matlab/process_field.m b/matlab/process_field.m
index 9390b1b2cd33e40186f458a292373df32dad64bd..9ead884c631372d4790cf003240b3e5aebbad2df 100644
--- a/matlab/process_field.m
+++ b/matlab/process_field.m
@@ -8,45 +8,44 @@ for i = 1:numel(OPTIONS.TIME)
[~,FRAMES(i)] =min(abs(OPTIONS.TIME(i)-DATA.Ts3D));
end
FRAMES = unique(FRAMES);
+TIME = DATA.Ts3D(FRAMES);
%% Setup the plot geometry
[KX, KY] = meshgrid(DATA.grids.kx, DATA.grids.ky);
directions = {'y','x','z'};
Nx = DATA.grids.Nx; Ny = DATA.grids.Ny; Nz = DATA.grids.Nz; Nt = numel(FRAMES);
-POLARPLOT = OPTIONS.POLARPLOT;
LTXNAME = OPTIONS.NAME;
+SKIP_COMP = 0; % Turned on only for kin. distr. func. plot;
+Z = 0;
switch OPTIONS.PLAN
case 'xy'
- XNAME = '$x$'; YNAME = '$y$';
+ XNAME = '$x/\rho_s$'; YNAME = '$y/\rho_s$';
[X,Y] = meshgrid(DATA.grids.x,DATA.grids.y);
REALP = 1; COMPDIM = 3; POLARPLOT = 0; SCALE = 1;
case 'xz'
- XNAME = '$x$'; YNAME = '$z$';
+ XNAME = '$x/\rho_s$'; YNAME = '$z/L_\parallel$';
[Y,X] = meshgrid(DATA.grids.z,DATA.grids.x);
- REALP = 1; COMPDIM = 1; SCALE = 0;
+ REALP = 1; COMPDIM = 1; SCALE = 0; POLARPLOT = 0;
case 'yz'
- XNAME = '$y$'; YNAME = '$z$';
+ XNAME = '$y/\rho_s$'; YNAME = '$z/L_\parallel$';
[Y,X] = meshgrid(DATA.grids.z,DATA.grids.y);
- REALP = 1; COMPDIM = 2; SCALE = 0;
+ REALP = 1; COMPDIM = 2; POLARPLOT = 0; SCALE = 0;
case 'kxky'
- XNAME = '$k_x$'; YNAME = '$k_y$';
+ XNAME = '$k_x\rho_s$'; YNAME = '$k_y\rho_s$';
[X,Y] = meshgrid(DATA.grids.kx,DATA.grids.ky);
REALP = 0; COMPDIM = 3; POLARPLOT = 0; SCALE = 1;
case 'kxz'
- XNAME = '$k_x$'; YNAME = '$z$';
+ XNAME = '$k_x\rho_s$'; YNAME = '$z/L_\parallel$';
[Y,X] = meshgrid(DATA.grids.z,DATA.grids.kx);
REALP = 0; COMPDIM = 1; POLARPLOT = 0; SCALE = 0;
case 'kyz'
- XNAME = '$k_y$'; YNAME = '$z$';
+ XNAME = '$k_y\rho_s$'; YNAME = '$z/L_\parallel$';
[Y,X] = meshgrid(DATA.grids.z,DATA.grids.ky);
REALP = 0; COMPDIM = 2; POLARPLOT = 0; SCALE = 0;
- case 'sx'
- XNAME = '$v_\parallel$'; YNAME = '$\mu$';
- [Y,X] = meshgrid(OPTIONS.XPERP,OPTIONS.SPAR);
- REALP = 1; COMPDIM = 3; POLARPLOT = 0; SCALE = 0;
- for i = 1:numel(OPTIONS.TIME)
- [~,FRAMES(i)] =min(abs(OPTIONS.TIME(i)-DATA.Ts5D));
- end
- FRAMES = unique(FRAMES); Nt = numel(FRAMES);
+ case 'xyz'
+ XNAME = '$x/\rho_s$'; YNAME = '$y/\rho_s$';
+ [X,Y] = meshgrid(DATA.grids.x,DATA.grids.y);
+ REALP = 1; COMPDIM = 3; POLARPLOT = 0; SCALE = 1;
+ SKIP_COMP = 1; OPTIONS.COMP = 3;
end
Xmax_ = max(max(abs(X))); Ymax_ = max(max(abs(Y)));
Lmin_ = min([Xmax_,Ymax_]);
@@ -73,29 +72,6 @@ else
FIELD = zeros(sz(1),sz(2),Nt);
end
%% Process the field to plot
-% short term writing --
-% b_e = DATA.sigma_e*sqrt(2*DATA.tau_e)*sqrt(KX.^2+KY.^2);
-% adiab_e = kernel(0,b_e);
-% pol_e = kernel(0,b_e).^2;
-% for n = 1:DATA.Jmaxe
-% adiab_e = adiab_e + kernel(n,b_e).^2;
-% pol_e = pol_e + kernel(n,b_e).^2;
-% end
-% adiab_e = DATA.q_e/DATA.tau_e .* adiab_e;
-% pol_e = DATA.q_e^2/DATA.tau_e * (1 - pol_e);
-%
-% b_i = DATA.sigma_i*sqrt(2*DATA.tau_i)*sqrt(KX.^2+KY.^2);
-% adiab_i = kernel(0,b_i);
-% pol_i = kernel(0,b_i).^2;
-% for n = 1:DATA.Jmaxi
-% adiab_i = adiab_i + kernel(n,b_i).^2;
-% pol_i = pol_i + kernel(n,b_i).^2;
-% end
-% pol_i = DATA.q_i^2/DATA.tau_i * (1 - pol_i);
-% adiab_i = DATA.q_i/DATA.tau_i .* adiab_i;
-% poisson_op = (pol_i + pol_e);
-adiab_e =0; adiab_i =0; poisson_op=1;
-SKIP_COMP = 0; % Turned on only for kin. distr. func. plot
% --
switch OPTIONS.COMP
case 'avg'
@@ -112,18 +88,18 @@ switch OPTIONS.COMP
i = OPTIONS.COMP;
compr = @(x) x(i,:,:);
if REALP
- COMPNAME = sprintf(['y=','%2.1f'],DATA.grids.x(i));
+ COMPNAME = sprintf(['y=','%2.1f'],DATA.grids.y(i));
else
- COMPNAME = sprintf(['k_y=','%2.1f'],DATA.grids.kx(i));
+ COMPNAME = sprintf(['k_y=','%2.1f'],DATA.grids.ky(i));
end
FIELDNAME = [LTXNAME,'(',COMPNAME,')'];
case 2
i = OPTIONS.COMP;
compr = @(x) x(:,i,:);
if REALP
- COMPNAME = sprintf(['x=','%2.1f'],DATA.grids.y(i));
+ COMPNAME = sprintf(['x=','%2.1f'],DATA.grids.x(i));
else
- COMPNAME = sprintf(['k_x=','%2.1f'],DATA.grids.ky(i));
+ COMPNAME = sprintf(['k_x=','%2.1f'],DATA.grids.kx(i));
end
FIELDNAME = [LTXNAME,'(',COMPNAME,')'];
case 3
@@ -188,6 +164,10 @@ switch OPTIONS.NAME
NAME = 'sx';
FLD_ = DATA.PHI(:,:,:,FRAMES);
OPE_ = KY.^2;
+ case 'w_{Ez}' % x-comp of ExB shear
+ NAME = 'wEz';
+ FLD_ = DATA.PHI(:,:,:,FRAMES);
+ OPE_ = (KX.^2+KY.^2).*(abs(KX)<=2).*(abs(KY)<=2);
case '\omega_z' % ES pot vorticity
NAME = 'vorticity';
FLD_ = DATA.PHI(:,:,:,FRAMES);
@@ -206,11 +186,11 @@ switch OPTIONS.NAME
OPE_ = 1;
case 'n_i' % ion density
NAME = 'ni';
- FLD_ = DATA.DENS_I(:,:,:,FRAMES) - adiab_i.* DATA.PHI(:,:,:,FRAMES);
+ FLD_ = DATA.DENS_I(:,:,:,FRAMES);
OPE_ = 1;
case 'n_e' % electron density
NAME = 'ne';
- FLD_ = DATA.DENS_E(:,:,:,FRAMES) - adiab_e.* DATA.PHI(:,:,:,FRAMES);
+ FLD_ = DATA.DENS_E(:,:,:,FRAMES);
OPE_ = 1;
case 'k^2n_e' % electron vorticity
NAME = 'k2ne';
@@ -219,12 +199,24 @@ switch OPTIONS.NAME
case 'n_i-n_e' % polarisation
NAME = 'pol';
OPE_ = 1;
- FLD_ = ((DATA.DENS_I(:,:,:,FRAMES)- adiab_i.* DATA.PHI(:,:,:,FRAMES))...
- -(DATA.DENS_E(:,:,:,FRAMES)- adiab_e.* DATA.PHI(:,:,:,FRAMES)));
+ FLD_ = DATA.DENS_I(:,:,:,FRAMES)...
+ -DATA.DENS_E(:,:,:,FRAMES);
+ case 'u_i' % ion density
+ NAME = 'ui';
+ FLD_ = DATA.UPAR_I(:,:,:,FRAMES);
+ OPE_ = 1;
+ case 'u_e' % electron density
+ NAME = 'ue';
+ FLD_ = DATA.UPAR_E(:,:,:,FRAMES);
+ OPE_ = 1;
case 'T_i' % ion temperature
NAME = 'Ti';
FLD_ = DATA.TEMP_I(:,:,:,FRAMES);
OPE_ = 1;
+ case 'T_e' % ion temperature
+ NAME = 'Te';
+ FLD_ = DATA.TEMP_E(:,:,:,FRAMES);
+ OPE_ = 1;
case 'G_x' % ion particle flux
NAME = 'Gx';
FLD_ = 0.*DATA.PHI(:,:,:,FRAMES);
@@ -238,24 +230,55 @@ switch OPTIONS.NAME
% tmp(:,:,iz) = abs(fftshift((squeeze(fft2(gx_,DATA.grids.Ny,DATA.grids.Nx))),2));
tmp(:,:,iz) = abs((squeeze(fft2(gx_,DATA.grids.Ny,DATA.grids.Nx))));
end
- FLD_(:,:,it)= squeeze(compr(tmp(1:DATA.grids.Nky,1:DATA.grids.Nkx,:)));
- end
- case 'Q_x' % ion heat flux
- NAME = 'Qx';
- % FLD_ = 0.*DATA.PHI(:,:,:,FRAMES);
+ FLD_(:,:,:,it)= squeeze((tmp(1:DATA.grids.Nky,1:DATA.grids.Nkx,:)));
+ end
+ case 'n_i T_i' % ion heat flux
+ NAME = 'niTi';
+ FLD_ = 0.*DATA.PHI(:,:,:,FRAMES);
OPE_ = 1;
for it = 1:numel(FRAMES)
- tmp = zeros(DATA.grids.Ny,DATA.grids.Nx,Nz);
+ qx_c = zeros(DATA.grids.Ny,DATA.grids.Nx,Nz);
for iz = 1:DATA.grids.Nz
- vx_ = real((ifourier_GENE(-1i*KY.*(DATA.PHI (:,:,iz,FRAMES(it))))));
- ni_ = real((ifourier_GENE( DATA.DENS_I(:,:,iz,FRAMES(it)))));
- Ti_ = real((ifourier_GENE( DATA.TEMP_I(:,:,iz,FRAMES(it)))));
- qx_ = vx_.*ni_.*Ti_;
-% tmp(:,:,iz) = abs(fftshift((squeeze(fft2(gx_,DATA.grids.Ny,DATA.grids.Nx))),2));
- tmp(:,:,iz) = abs((squeeze(fft2(qx_,DATA.grids.Ny,DATA.grids.Nx))));
+ ni_ = real((ifourier_GENE( DATA.DENS_I(:,:,iz,FRAMES(it)))));
+ Ti_ = real((ifourier_GENE( DATA.TEMP_I(:,:,iz,FRAMES(it)))));
+ cx_r = ni_.*Ti_;
+ cx_ = fft(cx_r,[],1);
+ cx_c(:,:,iz) = fft(cx_ ,[],2);
end
- FLD_(:,:,:,it)= squeeze(tmp(1:DATA.grids.Nky,1:DATA.grids.Nkx,:));
- end
+ FLD_(:,:,:,it)= squeeze(cx_c(1:DATA.grids.Nky,1:DATA.grids.Nkx,:))./numel(qx_c);
+ end
+ case 'Q_{xi}' % ion heat flux
+ NAME = 'Qxi';
+ FLD_ = 0.*DATA.PHI(:,:,:,FRAMES);
+ OPE_ = 1;
+ for it = 1:numel(FRAMES)
+ qx_c = zeros(DATA.grids.Ny,DATA.grids.Nx,Nz);
+ for iz = 1:DATA.grids.Nz
+ vx_ = real((ifourier_GENE(-1i*KY.*(DATA.PHI (:,:,iz,FRAMES(it))))));
+ ni_ = real((ifourier_GENE( DATA.DENS_I(:,:,iz,FRAMES(it)))));
+ Ti_ = real((ifourier_GENE( DATA.TEMP_I(:,:,iz,FRAMES(it)))));
+ qx_r = vx_.*ni_.*Ti_;
+ qx_ = fft(qx_r,[],1);
+ qx_c(:,:,iz) = fft(qx_ ,[],2);
+ end
+ FLD_(:,:,:,it)= squeeze(qx_c(1:DATA.grids.Nky,1:DATA.grids.Nkx,:))./numel(qx_c);
+ end
+ case 'Q_{xe}' % electron heat flux
+ NAME = 'Qxe';
+ FLD_ = 0.*DATA.PHI(:,:,:,FRAMES);
+ OPE_ = 1;
+ for it = 1:numel(FRAMES)
+ qx_c = zeros(DATA.grids.Ny,DATA.grids.Nx,Nz);
+ for iz = 1:DATA.grids.Nz
+ vx_ = real((ifourier_GENE(-1i*KY.*(DATA.PHI (:,:,iz,FRAMES(it))))));
+ ne_ = real((ifourier_GENE( DATA.DENS_E(:,:,iz,FRAMES(it)))));
+ Te_ = real((ifourier_GENE( DATA.TEMP_E(:,:,iz,FRAMES(it)))));
+ qx_r = vx_.*ne_.*Te_;
+ qx_ = fft(qx_r,[],1);
+ qx_c(:,:,iz) = fft(qx_ ,[],2);
+ end
+ FLD_(:,:,:,it)= squeeze(qx_c(1:DATA.grids.Nky,1:DATA.grids.Nkx,:))./numel(qx_c);
+ end
case 'f_i'
SKIP_COMP = 1;
shift_x = @(x) x; shift_y = @(y) y;
@@ -284,29 +307,34 @@ switch OPTIONS.NAME
end
% Process the field according to the 2D plane and the space (real/cpx)
if ~SKIP_COMP
-if COMPDIM == 3
- for it = 1:numel(FRAMES)
- tmp = squeeze(compr(OPE_.*FLD_(:,:,:,it)));
- FIELD(:,:,it) = squeeze(process(tmp));
- end
-else
- if REALP
- tmp = zeros(Ny,Nx,Nz);
+ if COMPDIM == 3
+ for it = 1:numel(FRAMES)
+ tmp = squeeze(compr(OPE_.*FLD_(:,:,:,it)));
+ FIELD(:,:,it) = squeeze(process(tmp));
+ end
else
- tmp = zeros(DATA.grids.Nky,DATA.grids.Nkx,Nz);
- end
- for it = 1:numel(FRAMES)
- for iz = 1:numel(DATA.grids.z)
- tmp(:,:,iz) = squeeze(process(OPE_.*FLD_(:,:,iz,it)));
+ if REALP
+ tmp = zeros(Ny,Nx,Nz);
+ else
+ tmp = zeros(DATA.grids.Nky,DATA.grids.Nkx,Nz);
end
- FIELD(:,:,it) = squeeze(compr(tmp));
- end
-end
+ for it = 1:numel(FRAMES)
+ for iz = 1:numel(DATA.grids.z)
+ tmp(:,:,iz) = squeeze(process(OPE_.*FLD_(:,:,iz,it)));
+ end
+ FIELD(:,:,it) = squeeze(compr(tmp));
+ end
+ end
+else
+ clear FIELD;
+ FIELD = squeeze(process(FLD_));
end
TOPLOT.FIELD = FIELD;
+TOPLOT.TIME = TIME;
TOPLOT.FRAMES = FRAMES;
TOPLOT.X = shift_x(X);
TOPLOT.Y = shift_y(Y);
+TOPLOT.Z = Z;
TOPLOT.FIELDNAME = FIELDNAME;
TOPLOT.XNAME = XNAME;
TOPLOT.YNAME = YNAME;
@@ -315,5 +343,6 @@ TOPLOT.DIMENSIONS= DIMENSIONS;
TOPLOT.ASPECT = ASPECT;
TOPLOT.FRAMES = FRAMES;
TOPLOT.INTERP = INTERP;
+
end
diff --git a/matlab/read_flux_out_XX.m b/matlab/read_flux_out_XX.m
index 2c7737ea4c28c895ca065fcae4df59af0b02c592..22f4ce581ac07712ec5492602e0ff5f354c28395 100644
--- a/matlab/read_flux_out_XX.m
+++ b/matlab/read_flux_out_XX.m
@@ -1,4 +1,4 @@
-function [t_all, Pxi_all, Qxi_all, Pxe_all, Qxe_all] = read_flux_out_XX(folderPath,PLOT)
+function [out] = read_flux_out_XX(folderPath,PLOT)
% Check if the prompt_string is provided as an argument
if nargin < 1
% If not provided, prompt the user for input
@@ -90,4 +90,11 @@ function [t_all, Pxi_all, Qxi_all, Pxe_all, Qxe_all] = read_flux_out_XX(folderPa
legend('show');
title('Radial heat flux')
end
+
+
+ out.t = t_all;
+ out.Pxi = Pxi_all;
+ out.Qxi = Qxi_all;
+ out.Pxe = Pxe_all;
+ out.Qxe = Qxe_all;
end