diff --git a/matlab/compute/compute_fa_2D.m b/matlab/compute/compute_fa_2D.m
index 01cf310596f9fb43e54fca71a6b80c3e218b805c..50b3593fa32ac761714a468f94a09c04722ab527 100644
--- a/matlab/compute/compute_fa_2D.m
+++ b/matlab/compute/compute_fa_2D.m
@@ -13,6 +13,7 @@ FaM = @(s,x) exp(-s.^2-x);
[~, ikx0] = min(abs(data.kx));
[~, iky0] = min(abs(data.ky));
+kx_ = data.kx; ky_ = data.ky;
switch options.SPECIE
case 'e'
@@ -28,9 +29,11 @@ end
switch options.Z
case 'avg'
Napj_ = mean(Napj_,5);
+ phi_ = mean(data.PHI,3);
otherwise
- [~,iz] = min(abs(options.Z-data.z));
+ iz = options.Z;
Napj_ = Napj_(:,:,:,:,iz,:);
+ phi_ = data.PHI(:,:,iz);
end
Napj_ = squeeze(Napj_);
@@ -38,6 +41,7 @@ frames = options.T;
for it = 1:numel(options.T)
[~,frames(it)] = min(abs(options.T(it)-data.Ts5D));
end
+frames = unique(frames);
Napj_ = mean(Napj_(:,:,:,:,frames),5);
@@ -46,32 +50,58 @@ Napj_ = squeeze(Napj_);
Np = numel(parray); Nj = numel(jarray);
-FF = zeros(data.Nkx,data.Nky,numel(options.XPERP),numel(options.SPAR));
-% FF = zeros(numel(options.XPERP),numel(options.SPAR));
-
-FAM = FaM(SS,XX);
-for ip_ = 1:Np
- p_ = parray(ip_);
- HH = Hp(p_,SS);
+if options.non_adiab
for ij_ = 1:Nj
- j_ = jarray(ij_);
- LL = Lj(j_,XX);
-% FF = FF + Napj_(ip_,ij_,ikx0,iky0)*HH.*LL.*FAM;
- HLF = HH.*LL.*FAM;
for ikx = 1:data.Nkx
- for iky = 1:data.Nky
- FF(ikx,iky,:,:) = squeeze(FF(ikx,iky,:,:)) + Napj_(ip_,ij_,ikx,iky)*HLF;
+ for iky = 1:data.Nky
+ kp_ = sqrt(kx_(ikx)^2 + ky_(iky)^2);
+ Napj_(1,ij_,ikx,iky) = Napj_(1,ij_,ikx,iky) + kernel(ij_,kp_)*phi_(ikx,iky);
end
end
end
end
+
+if options.RMS
+ FF = zeros(data.Nkx,data.Nky,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);
+ HLF = HH.*LL.*FAM;
+ for ikx = 1:data.Nkx
+ for iky = 1:data.Nky
+ FF(ikx,iky,:,:) = squeeze(FF(ikx,iky,:,:)) + Napj_(ip_,ij_,ikx,iky)*HLF;
+ end
+ 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|
-FF = sqrt(squeeze(mean(mean(FF,1),2))); %sqrt(<|f_a|>kx,ky)
+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 = FF./max(max(FF));
FF = FF';
-% FF = FF.*conj(FF);
% FF = sqrt(FF);
-% FF = FF./max(max(FF));
% FF = FF';
end
\ No newline at end of file
diff --git a/matlab/plot/plot_fa_gene.m b/matlab/plot/plot_fa_gene.m
new file mode 100644
index 0000000000000000000000000000000000000000..5a228a1f609a87c4f9fef44ea9a3d54b87292f19
--- /dev/null
+++ b/matlab/plot/plot_fa_gene.m
@@ -0,0 +1,97 @@
+function plot_fa_gene(OPTIONS)
+folder = OPTIONS.folder;
+TIMES = OPTIONS.times;
+specie = OPTIONS.specie;
+PLT_FCT= OPTIONS.PLT_FCT;
+
+file = 'coord.dat.h5';
+vp = h5read([folder,file],'/coord/vp');
+mu = h5read([folder,file],'/coord/mu');
+z = h5read([folder,file],'/coord/z');
+[XX,SS] = meshgrid(mu,vp);
+
+switch OPTIONS.Z
+ case 'avg'
+ zcomp_name = ' z-avg';
+ zcomp = @(x) squeeze(mean(x,1));
+ otherwise
+ zcomp_name = [' z=',sprintf('%2.2f',z(OPTIONS.Z))];
+ zcomp = @(x) squeeze(x(OPTIONS.Z,:,:));
+end
+
+[~,iv0] = min(abs(vp));
+[~,im0] = min(abs(mu));
+
+file = 'vsp.dat.h5';
+time = h5read([folder,file],'/vsp/time');
+
+fig = figure;
+
+G_t = [];
+Gdata = 0;
+for T = TIMES
+[~, it] = min(abs(time-T));
+tmp = h5read([folder,file],['/vsp/',OPTIONS.FIELD,'/',sprintf('%10.10d',it-1)]);
+Gdata = Gdata + tmp;
+G_t = [G_t time(it)];
+end
+Gdata = Gdata ./ numel(TIMES);
+
+if OPTIONS.ONED
+ switch specie
+ case 'e'
+ FFa = squeeze(Gdata(1,:,:,2));
+ FFa = abs(FFa)./max(max(abs(FFa)));
+ subplot(1,2,1)
+ plot(vp,FFa(:,im0)); hold on;
+ subplot(1,2,2)
+ plot(mu,FFa(iv0,:)); hold on;
+ case 'i'
+ FFa = squeeze(Gdata(1,:,:,1));
+ FFa = abs(FFa)./max(max(abs(FFa)));
+ end
+
+ subplot(1,2,1)
+ plot(vp,FFa(:,im0)); hold on;
+ legend(specie)
+ xlabel('$v_\parallel, (\mu=0)$'); ylabel('$\langle |f_a|^2\rangle_{xy}^{1/2}$');
+
+ subplot(1,2,2)
+ plot(mu,FFa(iv0,:)); hold on;
+ legend(specie)
+ xlabel('$\mu, (v_\parallel=0)$'); ylabel('$\langle |f_a|^2\rangle_{xy}^{1/2}$');
+else
+switch specie
+ case 'e'
+ name = '$f_e(v_\parallel,\mu_p)$';
+ FFa = zcomp(squeeze(Gdata));
+ FFa = abs(FFa)./max(max(abs(FFa)));
+ switch PLT_FCT
+ case 'contour'
+ contour(SS,XX,FFa,128);
+ case 'contourf'
+ pclr = contourf(SS,XX,FFa,128); set(pclr, 'edgecolor','none')
+ case 'pcolor'
+ pclr = pcolor(SS,XX,FFa); set(pclr, 'edgecolor','none'); shading interp
+ end
+ case 'i'
+ name = '$f_i(v_\parallel,\mu_p)$';
+ FFa = zcomp(squeeze(Gdata));
+ FFa = abs(FFa)./max(max(abs(FFa)));
+ switch PLT_FCT
+ case 'contour'
+ contour(SS,XX,FFa,128);
+ case 'contourf'
+ pclr = contourf(SS,XX,FFa,128); set(pclr, 'edgecolor','none')
+ case 'pcolor'
+ pclr = pcolor(SS,XX,FFa); set(pclr, 'edgecolor','none'); shading interp
+ end
+end
+if numel(TIMES) == 1
+ title(['Gene ',name,zcomp_name,', $t = ',sprintf('%2.2f',time(it)),'$']);
+else
+ title(['Gene ',name,zcomp_name,', averaged $t\in$[',sprintf('%2.2f',G_t(1)),',',sprintf('%2.2f',G_t(end)),']']);
+end
+
+clear FFi FFe tmp Gdata
+end
\ No newline at end of file
diff --git a/wk/gene_analysis_3D.m b/wk/gene_analysis_3D.m
index 6a08d4a4b05465b10b33e928fcf51891c667adb1..ed4d356579cfc22098722e57c9d79c22ef3ec25c 100644
--- a/wk/gene_analysis_3D.m
+++ b/wk/gene_analysis_3D.m
@@ -85,3 +85,97 @@ subplot(313)
xlim([min(gene_data.z),max(gene_data.z)]); legend('show');
end
+% folder = '/misc/gene_results/HP_fig_2b_mu_5e-2/';
+folder = '/misc/gene_results/shearless_cyclone/s_alpha_output_1.0/';
+% folder = '/misc/gene_results/cyclone/s_alpha_output_1.0/';
+gene_data = load_gene_data(folder);
+gene_data = rotate_c_plane_nxnky_to_nkxny(gene_data);
+if 1
+%% Space time diagramm (fig 11 Ivanov 2020)
+TAVG_0 = 0.8*gene_data.Ts3D(end); TAVG_1 = gene_data.Ts3D(end); % Averaging times duration
+% chose your field to plot in spacetime diag (uzf,szf,Gx)
+field = 'phi';
+compz = 'avg';
+nmvm = 1;
+fig = plot_radial_transport_and_spacetime(gene_data,TAVG_0,TAVG_1,field,nmvm,compz);
+% save_figure(data,fig)
+end
+
+if 0
+%% 2D snapshots
+% Options
+options.INTERP = 1;
+options.POLARPLOT = 0;
+options.AXISEQUAL = 1;
+options.NAME = '\phi';
+% options.NAME = 'n_i';
+% options.NAME = 'T_i';
+% options.NAME = '\Gamma_x';
+% options.NAME = 'k^2n_e';
+options.PLAN = 'xz';
+% options.NAME ='f_e';
+% options.PLAN = 'sx';
+options.COMP = 'avg';
+options.TIME = [100 300 900];
+gene_data.a = data.EPS * 2000;
+fig = photomaton(gene_data,options);
+save_figure(gene_data,fig)
+end
+
+if 0
+%% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Options
+options.INTERP = 0;
+options.POLARPLOT = 0;
+options.NAME = '\phi';
+% options.NAME = 'v_y';
+% options.NAME = 'n_i^{NZ}';
+% options.NAME = '\Gamma_x';
+% options.NAME = 'n_i';
+options.PLAN = 'xy';
+% options.NAME = 'f_e';
+% options.PLAN = 'sx';
+options.COMP = 'avg';
+options.TIME = 000:170;
+gene_data.a = data.EPS * 2000;
+create_film(gene_data,options,'.gif')
+end
+
+if 0
+%% Geometry
+names = {'$g^{xx}$','$g^{xy}$','$g^{xz}$','$g^{yy}$','$g^{yz}$','$g^{zz}$',...
+ '$B_0$','$\partial_x B_0$','$\partial_y B_0$','$\partial_z B_0$',...
+ '$J$','$R$','$\phi$','$Z$','$\partial_R x$','$\partial_Z x$'};
+figure;
+subplot(311)
+ for i = 1:6
+ plot(gene_data.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on;
+ end
+ xlim([min(gene_data.z),max(gene_data.z)]); legend('show'); title('GENE geometry');
+
+subplot(312)
+ for i = 7:10
+ plot(gene_data.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on;
+ end
+ xlim([min(gene_data.z),max(gene_data.z)]); legend('show');
+
+subplot(313)
+ for i = 11:16
+ plot(gene_data.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on;
+ end
+ xlim([min(gene_data.z),max(gene_data.z)]); legend('show');
+
+end
+
+if 0
+%% Show f_i(vpar,mu)
+options.times = 200:600;
+options.specie = 'i';
+options.PLT_FCT = 'pcolor';
+options.folder = folder;
+options.Z = 'avg';
+options.FIELD = '<f_>';
+options.ONED = 0;
+% options.FIELD = 'Q_es';
+plot_fa_gene(options);
+end