Homec4science
Diffusion FENNECS: Finite Element Non-Neutral Electron Cloud Simulator (master) Browse

fennecs/matlab/fft_espic2d.m
ddfe457d9de0master

Last Change
Raw File

fft_espic2d.m
View Options

file='stable_13_1058400.h5';
%file='stable_13_decoupled_fft.h5';
%file='teststablegauss_13_traject2.h5';
%file='stable_13_decoupled_fft.h5';
%file='stablegauss_13_traject.h5';
%close all hidden;
if (~exist('M','var'))
M.file=file;
end
M=load_espic2d(file,M,'all');
%close all hidden;
fig=figure;
ax1=gca;
pcolor(ax1,M.zgrid,M.rgrid,M.N(:,:,end));
hold on
[r,z]=find(M.N(:,:,end)~=0);
xlim(ax1,[M.zgrid(min(z)) M.zgrid(max(z))])
ylim(ax1,[M.rgrid(min(r)) M.rgrid(max(r))])
xlabel(ax1,'Z [m]')
ylabel(ax1,'R [m]')
c = colorbar(ax1);
c.Label.String= 'n [m^{-3}]';
view(ax1,2)
%set(ax1,'colorscale','log')
stepmin=4000;
stepend=size(M.N,3);
papsize=[16 14];
% [x,y]=ginput(1);
% nz=find(x>M.zgrid,1,'last');
% nr=find(y>M.rgrid,1,'last');
% decoupled
%nr=18; nz=127;
%nr=18; nz=76;
%stable
nr=18; nz=63;
%nr=18; nz=38;
Z=(M.zgrid(nz)+M.zgrid(nz+1))/2;
R=(M.rgrid(nr)+M.rgrid(nr+1))/2;
plot(Z,R,'rx','Markersize',12);
%decoupled
% nr=20;nz=129;
%nr=18; nz=127;
%nr=18; nz=76;
%nr=19; nz=127;
%nr=16;nz=128;
% stable
%nr=18; nz=63;
nr=18; nz=38;
Z=(M.zgrid(nz)+M.zgrid(nz+1))/2;
R=(M.rgrid(nr)+M.rgrid(nr+1))/2;
plot(Z,R,'rx','Markersize',12);
fig.PaperOrientation='landscape';
[~, name, ~] = fileparts(M.file);
fig.PaperUnits='centimeters';
fig.PaperSize=papsize;
print(fig,sprintf('%s_fftpoints',name),'-dpdf','-fillpage')
%% Compute standard deviation
figure;
ax1=gca;
Nmean=mean(M.N,3);
Nmean=1./Nmean;
Nmean(isnan(Nmean))=0;
surfc(ax1,M.zgrid,M.rgrid,std(M.N(:,:,stepmin:end),0,3));
hold on
[r,z]=find(M.N(:,:,end)~=0);
xlim(ax1,[M.zgrid(min(z)) M.zgrid(max(z))])
ylim(ax1,[M.rgrid(min(r)) M.rgrid(max(r))])
xlabel(ax1,'Z [m]')
ylabel(ax1,'R [m]')
c = colorbar(ax1);
c.Label.String= '\sigma_n [m^{-3}]';
%set(ax1,'colorscale','log')
view(ax1,2)
figure;
ax1=gca;
Nmean=mean(M.N,3);
Nmean=1./Nmean;
Nmean(isnan(Nmean))=0;
surfc(ax1,M.zgrid,M.rgrid,std(M.pot(:,:,stepmin:end),0,3));
hold on
xlim(ax1,[M.zgrid(min(z)) M.zgrid(max(z))])
ylim(ax1,[M.rgrid(min(r)) M.rgrid(max(r))])
xlabel(ax1,'Z [m]')
ylabel(ax1,'R [m]')
c = colorbar(ax1);
c.Label.String= '\sigma_\phi [V]';
%set(ax1,'colorscale','log')
view(ax1,2)
%% Start the fourrier transform
dt=M.dt*double(M.it1);
FS=1/dt;
L=stepend-stepmin;
f=FS*double(0:L/2)/L/1e9;
%S=M.N(nr,nz,stepmin:end);
S=M.pot(nr,nz,stepmin:stepend);
S=detrend(S(:),0);
S=S/max(S);
S=lowpass(S(:),.499*FS,FS);
P2=abs(fft(S(:)))/L;
P1=P2(1:floor(L/2)+1);
P1(2:end-1) = 2*P1(2:end-1);
tmin=stepmin*dt;
tmax=dt*size(M.N,3);
%% Compute the relevant frequencies
fce=M.qe/M.me/2/pi*M.Bz(nz,nr)/1e9;
N=M.N(:,:,end);
fpe=sqrt(mean(M.N(nr,nz,stepmin:end))*M.qe^2/(M.eps_0*M.me))/(2*pi*1e9);
fpemin=sqrt(min(min(N(N>0.01e14)))*M.qe^2/(M.eps_0*M.me))/(2*pi*1e9);
fpemax=sqrt(max(max(N))*M.qe^2/(M.eps_0*M.me))/(2*pi*1e9);
freplus=-sign(M.qsim)*fce/2*(1+sqrt(1-2*fpe^2/fce^2));
freminus=-sign(M.qsim)*fce/2*(1-sqrt(1-2*fpe^2/fce^2));
zIndices=mean(M.Presstens(6,nr,:,stepmin:end),4)>0;
lengths=Blocklength(mean(M.N(nr,:,stepmin:end),3)>1e13);
%zlength=M.width*2;
zlength=max(lengths)*M.dz;
fz=sqrt(mean(mean(M.Presstens(6,nr,zIndices,stepmin:end))))*M.vlight/zlength/1e9;
fig=figure();
%% Plot the fft
papsize=[14 16];
subplot(3,1,1)
pl=plot(f(1:end),P1(1:end),'linewidth',1.2,'Displayname','Spectrum');
title(sprintf('r=%3.2e[m],z=%3.2e[m], t=[%3.2e,%3.2e][s]',R,Z,tmin,tmax))
xlabel('f [GHz]','Interpreter','latex','fontsize',14)
ylabel('$\left|\hat{\Phi}(f)\right|$ [a.u.]','Interpreter','latex','fontsize',14)
grid on
ylimits=ylim();
hold on
plot(freminus*[1 1],ylimits,'k-.','Displayname','f_{re}^-','linewidth',1.2);
plot(freplus*[1 1],ylimits,'k-','Displayname','f_{re}^+','linewidth',1.2);
plot(fce*[1 1],ylimits,'r--','Displayname','f_{ce}','linewidth',1.2);
plot(fz*[1 1],ylimits,'k:','Displayname','f_{z}','linewidth',1.2);
plot(fpe*[1 1],ylimits,'r:','Displayname','local f_{pe}','linewidth',1.2);
fill([fpemin fpemax fpemax fpemin],[ylimits(1) ylimits(1) ylimits(2) ylimits(2)],'r','FaceColor',[1 0 0]*0.5,'FaceAlpha',.3,'Displayname','f_{pe} range');
legend('location','eastoutside')
uistack(pl,'top');
ylim(ylimits)
xlim([min(f),1.5*max(fz,fpemax)])
%%
subplot(3,1,2)
pl=plot(f(1:end),P1(1:end),'linewidth',1.2,'Displayname','Spectrum');
xlabel('f [GHz]','Interpreter','latex','fontsize',14)
ylabel('$\left|\hat{\Phi}(f)\right|$ [a.u.]','Interpreter','latex','fontsize',14)
grid on
%xlim([0 10])
ylimits=ylim();
hold on
plot(freminus*[1 1],ylimits,'k-.','Displayname','f_{re}^-','linewidth',1.2);
plot(freplus*[1 1],ylimits,'k-','Displayname','f_{re}^+','linewidth',1.2);
plot(fce*[1 1],ylimits,'r--','Displayname','f_{ce}','linewidth',1.2);
plot(fz*[1 1],ylimits,'k:','Displayname','f_{z}','linewidth',1.2);
plot(fpe*[1 1],ylimits,'r:','Displayname','f_{pe}','linewidth',1.2);
fill([fpemin fpemax fpemax fpemin],[ylimits(1) ylimits(1) ylimits(2) ylimits(2)],'r','FaceColor',[1 0 0]*0.5,'FaceAlpha',.3,'Displayname','f_{pe}');
uistack(pl,'top');
ylim(ylimits)
xlim([min(f),max(f)])
subplot(3,1,3)
Signal=M.pot(nr,nz,1:stepend);
%Signal=M.N(nr,nz,stepmin:end);
pl=plot(dt*((1:length(Signal))-1),Signal(:));
xlabel('t [s]','Interpreter','latex','fontsize',14)
ylabel('$\Phi(t)$ [V]','Interpreter','latex','fontsize',14)
ylimits=ylim();
%ylimits=[-56 -54];
hold on
x=[stepmin-1 length(Signal)-1 length(Signal)-1 stepmin-1]*dt;
y=[ylimits(1) ylimits(1) ylimits(2) ylimits(2)];
fill(x,y,'r','FaceColor',[1 1 1]*0.5,'FaceAlpha',.3)
ylim(ylimits)
uistack(pl,'top');
%fig.PaperOrientation='landscape';
[~, name, ~] = fileparts(M.file);
fig.PaperUnits='centimeters';
fig.PaperSize=papsize;
print(fig,sprintf('%s_fftPhi_R%dZ%d',name,nr,nz),'-dpdf','-fillpage')
S=M.N(nr,nz,stepmin:stepend);
S=detrend(S(:),0);
S=S/max(S);
S=lowpass(S(:),.499*FS,FS);
P2=abs(fft(S(:)))/L;
P1=P2(1:floor(L/2)+1);
P1(2:end-1) = 2*P1(2:end-1);
fig=figure();
%%
subplot(3,1,1)
pl=plot(f(1:end),P1(1:end),'linewidth',1.2,'Displayname','Spectrum');
title(sprintf('r=%3.2e[m],z=%3.2e[m], t=[%3.2e,%3.2e][s]',R,Z,tmin,tmax))
xlabel('f [GHz]','Interpreter','latex','fontsize',14)
ylabel('$\left|\hat{n}(f)\right|$ [a.u.]','Interpreter','latex','fontsize',14)
grid on
%xlim([0 10])
ylimits=ylim();
hold on
plot(freminus*[1 1],ylimits,'k-.','Displayname','f_{re}^-','linewidth',1.2);
plot(freplus*[1 1],ylimits,'k-','Displayname','f_{re}^+','linewidth',1.2);
plot(fce*[1 1],ylimits,'r--','Displayname','f_{ce}','linewidth',1.2);
plot(fz*[1 1],ylimits,'k:','Displayname','f_{z}','linewidth',1.2);
plot(fpe*[1 1],ylimits,'r:','Displayname','local f_{pe}','linewidth',1.2);
fill([fpemin fpemax fpemax fpemin],[ylimits(1) ylimits(1) ylimits(2) ylimits(2)],'r','FaceColor',[1 0 0]*0.5,'FaceAlpha',.3,'Displayname','range f_{pe}');
legend('location','eastoutside')
uistack(pl,'top');
xlim([min(f),1.5*max(fz,fpemax)])
%%
subplot(3,1,2)
pl=plot(f(1:end),P1(1:end),'linewidth',1.2,'Displayname','Spectrum');
xlabel('f [GHz]','Interpreter','latex','fontsize',14)
ylabel('$\left|\hat{n}(f)\right|$ [a.u.]','Interpreter','latex','fontsize',14)
grid on
%xlim([0 10])
ylimits=ylim();
hold on
plot(freminus*[1 1],ylimits,'k-.','Displayname','f_{re}^-','linewidth',1.2);
plot(freplus*[1 1],ylimits,'k-','Displayname','f_{re}^+','linewidth',1.2);
plot(fce*[1 1],ylimits,'r--','Displayname','f_{ce}','linewidth',1.2);
plot(fz*[1 1],ylimits,'k:','Displayname','f_{z}','linewidth',1.2);
plot(fpe*[1 1],ylimits,'r:','Displayname','f_{pe}','linewidth',1.2);
fill([fpemin fpemax fpemax fpemin],[ylimits(1) ylimits(1) ylimits(2) ylimits(2)],'r','FaceColor',[1 0 0]*0.5,'FaceAlpha',.3,'Displayname','f_{pe}');
uistack(pl,'top');
xlim([min(f),max(f)])
%%
subplot(3,1,3)
Signal=M.N(nr,nz,1:stepend);
%Signal=M.N(nr,nz,stepmin:end);
pl=plot(dt*((1:length(Signal))-1),Signal(:));
xlabel('t [s]','Interpreter','latex','fontsize',14)
ylabel('$n(t)$ [m$^{-3}$]','Interpreter','latex','fontsize',14)
ylimits=ylim();
hold on
x=[stepmin-1 length(Signal)-1 length(Signal)-1 stepmin-1]*dt;
y=[ylimits(1) ylimits(1) ylimits(2) ylimits(2)];
fill(x,y,'r','FaceColor',[1 1 1]*0.5,'FaceAlpha',.3)
uistack(pl,'top');
%fig.PaperOrientation='landscape';
[~, name, ~] = fileparts(M.file);
fig.PaperUnits='centimeters';
fig.PaperSize=papsize;
ylim(ylimits)
print(fig,sprintf('%s_fftdens_R%dZ%d',name,nr,nz),'-dpdf','-fillpage')
c4science · Help