function [ TOPLOT ] = process_field( DATA, OPTIONS )
%UNTITLED4 Summary of this function goes here
% Detailed explanation goes here
%% Setup time
%%
FRAMES = zeros(size(OPTIONS.TIME));
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);
LTXNAME = OPTIONS.NAME;
SKIP_COMP = 0; % Turned on only for kin. distr. func. plot;
Z = 0;
switch OPTIONS.PLAN
case 'xy'
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/\rho_s$'; YNAME = '$z/L_\parallel$';
[Y,X] = meshgrid(DATA.grids.z,DATA.grids.x);
REALP = 1; COMPDIM = 1; SCALE = 0; POLARPLOT = 0;
case 'yz'
XNAME = '$y/\rho_s$'; YNAME = '$z/L_\parallel$';
[Y,X] = meshgrid(DATA.grids.z,DATA.grids.y);
REALP = 1; COMPDIM = 2; POLARPLOT = 0; SCALE = 0;
case 'kxky'
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\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\rho_s$'; YNAME = '$z/L_\parallel$';
[Y,X] = meshgrid(DATA.grids.z,DATA.grids.ky);
REALP = 0; COMPDIM = 2; POLARPLOT = 0; SCALE = 0;
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_]);
Rx = Xmax_/Lmin_ * SCALE + (1-SCALE)*1.2;
Ry = Ymax_/Lmin_ * SCALE + (1-SCALE)*1.2;
ASPECT = [Rx, Ry, 1];
DIMENSIONS = [500, 1000, OPTIONS.RESOLUTION*Rx, OPTIONS.RESOLUTION*Ry];
% Polar grid
POLARNAME = [];
if POLARPLOT
POLARNAME = 'polar';
X__ = (X+DATA.a).*cos(Y);
Y__ = (X+DATA.a).*sin(Y);
X = X__;
Y = Y__;
XNAME='X';
YNAME='Z';
DIMENSIONS = [100, 100, OPTIONS.RESOLUTION, OPTIONS.RESOLUTION];
ASPECT = [1,1,1];
sz = size(X);
FIELD = zeros(sz(1),sz(2),Nt);
else
sz = size(X);
FIELD = zeros(sz(1),sz(2),Nt);
end
%% Process the field to plot
% --
switch OPTIONS.COMP
case 'avg'
compr = @(x) mean(x,COMPDIM);
COMPNAME = ['avg',directions{COMPDIM}];
FIELDNAME = ['\langle ',LTXNAME,'\rangle_',directions{COMPDIM}];
case 'max'
compr = @(x) max(x,[],COMPDIM);
COMPNAME = ['max',directions{COMPDIM}];
FIELDNAME = ['\max_',directions{COMPDIM},' ',LTXNAME];
otherwise
switch COMPDIM
case 1
i = OPTIONS.COMP;
compr = @(x) x(i,:,:);
if REALP
COMPNAME = sprintf(['y=','%2.1f'],DATA.grids.y(i));
else
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.x(i));
else
COMPNAME = sprintf(['k_x=','%2.1f'],DATA.grids.kx(i));
end
FIELDNAME = [LTXNAME,'(',COMPNAME,')'];
case 3
i = OPTIONS.COMP;
compr = @(x) x(:,:,i);
COMPNAME = sprintf(['z=','%2.1f'],DATA.grids.z(i));
FIELDNAME = [LTXNAME,'(',COMPNAME,')'];
end
end
switch REALP
case 1 % Real space plot
INTERP = OPTIONS.INTERP;
process = @(x) real(fftshift(ifourier_GENE(x)));
% process = @(x) real(fftshift(ifft2(x,sz(1),sz(2))));
shift_x = @(x) x;
shift_y = @(x) x;
case 0 % Frequencies plot
INTERP = 0;
switch COMPDIM
case 1
process = @(x) abs(fftshift(x,2));
shift_x = @(x) fftshift(x,1);
shift_y = @(x) fftshift(x,1);
case 2
process = @(x) abs((x));
shift_x = @(x) (x);
shift_y = @(x) (x);
case 3
process = @(x) abs(fftshift(x,2));
shift_x = @(x) fftshift(x,2);
shift_y = @(x) fftshift(x,2);
end
end
switch OPTIONS.NAME
case '\phi' %ES pot
NAME = 'phi';
FLD_ = DATA.PHI(:,:,:,FRAMES); % data to plot
OPE_ = 1; % Operation on data
case '\psi' %EM pot
NAME = 'psi';
FLD_ = DATA.PSI(:,:,:,FRAMES);
OPE_ = 1;
case '\phi^{NZ}' % non-zonal ES pot
NAME = 'phiNZ';
FLD_ = DATA.PHI(:,:,:,FRAMES);
OPE_ = (KY~=0);
case 'v_{Ey}' % y-comp of ExB velocity
NAME = 'vy';
FLD_ = DATA.PHI(:,:,:,FRAMES);
OPE_ = -1i*KX;
case 'v_{Ex}' % x-comp of ExB velocity
NAME = 'vx';
FLD_ = DATA.PHI(:,:,:,FRAMES);
OPE_ = -1i*KY;
case 's_{Ey}' % y-comp of ExB shear
NAME = 'sy';
FLD_ = DATA.PHI(:,:,:,FRAMES);
OPE_ = KX.^2;
case 's_{Ex}' % x-comp of ExB shear
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);
OPE_ = -(KX.^2+KY.^2);
case 'N_i^{00}' % first ion gyromoment
NAME = 'Ni00';
FLD_ = DATA.Ni00(:,:,:,FRAMES);
OPE_ = 1;
case 'N_e^{00}' % first electron gyromoment
NAME = 'Ne00';
FLD_ = DATA.Ne00(:,:,:,FRAMES);
OPE_ = 1;
case 'N_i^{00}-N_e^{00}' % first electron gyromoment
NAME = 'Ni00-Ne00';
FLD_ = (DATA.Ni00(:,:,:,FRAMES)-DATA.Ne00(:,:,:,FRAMES))./(poisson_op+(poisson_op==0));
OPE_ = 1;
case 'n_i' % ion density
NAME = 'ni';
FLD_ = DATA.DENS_I(:,:,:,FRAMES);
OPE_ = 1;
case 'n_e' % electron density
NAME = 'ne';
FLD_ = DATA.DENS_E(:,:,:,FRAMES);
OPE_ = 1;
case 'k^2n_e' % electron vorticity
NAME = 'k2ne';
FLD_ = DATA.DENS_E(:,:,:,FRAMES);
OPE_ = -(KX.^2+KY.^2);
case 'n_i-n_e' % polarisation
NAME = 'pol';
OPE_ = 1;
FLD_ = DATA.DENS_I(:,:,:,FRAMES)...
-DATA.DENS_E(:,:,:,FRAMES);
case 'upar_i' % ion density
NAME = 'upari';
FLD_ = DATA.UPAR_I(:,:,:,FRAMES);
OPE_ = 1;
case 'upar_e' % electron density
NAME = 'upare';
FLD_ = DATA.UPAR_E(:,:,:,FRAMES);
OPE_ = 1;
case 'uper_i' % ion density
NAME = 'uperi';
FLD_ = DATA.UPER_I(:,:,:,FRAMES);
OPE_ = 1;
case 'uper_e' % electron density
NAME = 'upere';
FLD_ = DATA.UPER_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);
OPE_ = 1;
for it = 1:numel(FRAMES)
tmp = 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)))));
gx_ = vx_.*ni_;
% 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((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)
qx_c = zeros(DATA.grids.Ny,DATA.grids.Nx,Nz);
for iz = 1:DATA.grids.Nz
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(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;
NAME = 'fi'; OPTIONS.SPECIE = 'i';
for it = 1:numel(FRAMES)
OPTIONS.T = DATA.Ts5D(FRAMES(it));
OPTIONS.Z = OPTIONS.COMP;
[~,~,FIELD(:,:,it)] = compute_fa(DATA,OPTIONS);
end
case 'f_e'
SKIP_COMP = 1;
shift_x = @(x) x; shift_y = @(y) y;
NAME = 'fe'; OPTIONS.SPECIE = 'e';
[~,it0_] =min(abs(OPTIONS.TIME(1)-DATA.Ts5D));
[~,it1_]=min(abs(OPTIONS.TIME(end)-DATA.Ts5D));
dit_ = 1;%ceil((it1_-it0_+1)/10);
FRAMES = it0_:dit_:it1_;
iz = 1;
for it = 1:numel(FRAMES)
OPTIONS.T = DATA.Ts5D(FRAMES(it));
[~,~,FIELD(:,:,it)] = compute_fa(DATA,OPTIONS);
end
otherwise
disp('Fieldname not recognized');
return
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);
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)));
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;
TOPLOT.FILENAME = [NAME,'_',OPTIONS.PLAN,'_',COMPNAME,'_',POLARNAME];
TOPLOT.DIMENSIONS= DIMENSIONS;
TOPLOT.ASPECT = ASPECT;
TOPLOT.FRAMES = FRAMES;
TOPLOT.INTERP = INTERP;
end