Homec4science
Diffusion GYACOMO (OBSOLETE: LATEST VERSION ON gitlab.epfl.ch) (master) Browse

gyacomo/matlab/ZF_fourier_analysis.m
0d685e3b68a2master

Last Change
Raw File

ZF_fourier_analysis.m
View Options

%% Zonal flow spectral analysis
fig = figure; FIGNAME = ['zonal_flow_spectral_analysis_',PARAMS];
tend = TAVG_1; tstart = TAVG_0;
[~,its0D] = min(abs(Ts0D-tstart));
[~,ite0D] = min(abs(Ts0D-tend));
[~,its3D] = min(abs(Ts3D-tstart));
[~,ite3D] = min(abs(Ts3D-tend));
TAVG = Ts3D(ite3D)-Ts3D(its3D);
set(gcf, 'Position', [100, 100, 800, 400])
% Time series analysis (burst period and time frequencies spectrum)
subplot(121)
samplerate = Ts0D(2)-Ts0D(1);
Y = log(PGAMMA_RI(its0D:ite0D)*(2*pi/Nx/Ny)^2);
[n,~] = size(Y);
Yy= fft(Y);
Pot = Yy .* conj(Yy) / n;
Freq = (samplerate / n * (1:n))';
Pot(1) = 0;
nmax = min(20,round(n/2));
[amax, itmax] = max(Pot);
% plot((0:nmax-1) , Pot(1:nmax)/amax,'DisplayName','$\Gamma_x(\omega)$');hold on;
% plot([itmax-1,itmax-1],[0,1],'--k', 'DisplayName',['$T_{per}\approx',num2str(round(1/Freq(itmax))),'L_\perp/c_s$']);
semilogx(TAVG./(0:nmax-1) , Pot(1:nmax)/amax,'o-','DisplayName','$\Gamma_x(\omega)$');hold on;
semilogx(TAVG./[itmax-1,itmax-1],[0,1],'--k', 'DisplayName',['$T_{max}\approx',num2str(round(1/Freq(itmax-1))),'L_\perp/c_s$']);
legend('show'); grid on; box on; xlabel('Period $Tc_s/R$'); yticks([]);
title('$\Gamma_x$ temporal spectrum')
% Space analysis (spatial period of ZF)
subplot(122)
nmax = 20; n = numel(x);
[TT,NN] = meshgrid(Ts3D(its3D:ite3D),0:n-1);
Pot = NN;
for it = 1:ite3D-its3D+1
Y = mean(real(dxphi(:,:,it)),2);
Yy = fft(Y);
[n,~] = size(Yy);
Pot(:,it) = Yy .* conj(Yy) / n;
end
[amax, ikZF] = max(mean(Pot,2));
% pclr = pcolor(NN(1:nmax,:),TT(1:nmax,:),Pot(1:nmax,:)); set(pclr, 'edgecolor','none'); hold on;
% plot(0:nmax,mean(Pot(1:nmax+1,:),2)/amax,'DisplayName','$\langle\partial_x\phi\rangle_y (k_x)$'); hold on;
% plot([ikZF-1,ikZF-1],[0,1],'--k', 'DisplayName',['$L_x=',num2str(2*pi/kx(ikZF)),'\rho_s$']);
semilogx(Lx./(0:nmax),mean(Pot(1:nmax+1,:),2)/amax,'o-','DisplayName','$\langle\partial_x\phi\rangle_y (k_x)$'); hold on;
semilogx(Lx./[ikZF-1,ikZF-1],[0,1],'--k', 'DisplayName',['$L_x=',num2str(2*pi/kx(ikZF)),'\rho_s$']);
grid on; box on;
title('ZF spatial spectrum')
xlabel('Period $\lambda/\rho_s$'); yticks([]); legend('show')
save_figure
%% Pred-Pray phase space (A Zonal Flow review, Diamond 2005, Fig 15, Kobayashi 2015)
% Time evol. of the turbulence energy (Pred in Kobayashi 2015, N = sum phi_k^2 (1+k^2) Non zonal)
E_turb = zeros(1,Ns3D);
% Time evol. of the ZF energy (Pray in Kobayashi 2015, Ev = phi_q^2 q^2)
E_ZF = zeros(1,Ns3D);
for it = 1:numel(Ts3D)
E_turb(it) = sum(sum(((KY~=0).*(1+KX.^2+KY.^2).*abs(PHI(:,:,1,it)).^2)));
% E_ZF(it) = kx(ikZF)^2*abs(PHI(ikZF,1,1,it)).^2;
E_ZF(it) = sum(sum(((KY==0).*(1+KX.^2+KY.^2).*abs(PHI(:,:,1,it)).^2)));
end
fig = figure; FIGNAME = ['phi_shear_phase_space_',PARAMS];
set(gcf, 'Position', [100, 100, 700, 500])
scatter(E_ZF*SCALE,E_turb*SCALE,80,Ts3D,'.',...
'DisplayName',PARAMS); cbar = colorbar;ylabel(cbar,'$t c_s/\rho_s$','Interpreter','LaTeX')
hold on
% xlabel('$\langle \phi \rangle_z^r$'); ylabel('$\langle dV_E/dr \rangle_z^r$')
xlabel('$E_v$'); ylabel('$N$')
grid on; title('ES pot. vs Shear phase space')
% plot(phi_avgr_maxz(iTs3D:ite2D),shear_avgr_maxz(iTs3D:ite2D),'-')
% plot(phi_maxr_maxz(iTs3D:ite2D),shear_maxr_maxz(iTs3D:ite2D),'-')
% plot(phi_avgr_avgz(iTs3D:ite2D),shear_avgr_avgz(iTs3D:ite2D),'-')
save_figure
clear x_ y_
c4science · Help