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

fennecs/matlab/Continuityspline.m
438cf8d62988master

Last Change
Raw File

Continuityspline.m
View Options

function [Continuity]=Continuityspline(M,it,norm,log,contourscale, zlims)
%ForceBalance Show the radial force balance
% Plot the three radial forces for the given time-step it or averaged
% over the range of time steps defined in it
if strcmp(it,':')
it=floor(0.95*size(M.t2d)):size(M.t2d);
end
switch nargin
case 2
contourscale=0.2;
norm=false;
log=false;
zlims=[-inf inf];
case 3
log=false;
contourscale=0.2;
zlims=[-inf inf];
case 4
contourscale=0.2;
zlims=[-inf inf];
case 5
zlims=[-inf inf];
case 6
otherwise
error("Invalid number of arguments")
end
contourcolor=[0 0 0];%[255 20 147]/255;
if(it(1)==1)
it=it(2:end);
end
if(it(end)==length(M.t2d))
it=it(1:end-1);
end
N=M.N(:,:,it);
maxdens=max(N(:));
densitycontour=contourc(M.zgrid,M.rgrid,mean(N,3),[contourscale contourscale]*maxdens);
rgridmax=sum(M.nnr(1:2));
%% Continuity
t=M.t2d(it);
dndt=zeros(size(N,1),size(N,2),size(N,3));
S=zeros(size(N,1),size(N,2),size(N,3));
ndivu=zeros(size(N,1),size(N,2),size(N,3),3);
ugradn=zeros(size(N,1),size(N,2),size(N,3),2);
[R,~]=meshgrid(M.rgrid,M.zgrid);
R=R';
Rinv=1./R;
Rinv(isinf(Rinv))=0;
geommask=M.geomweight(:,:,1)<0;
for i=1:size(N,3)
ddt=(M.N(:,:,it(i)+1)-M.N(:,:,it(i)-1))./(M.t2d(it(i)+1)-M.t2d(it(i)-1));
ddt(geommask)=0;
dndt(:,:,i)=ddt;
ndivu1=N(:,:,i).*(M.fluidUR.der(:,:,it(i),[1 0]));
ndivu1(geommask)=0;
ndivu(:,:,i,1)=ndivu1;
ndivu2=N(:,:,i).*(M.fluidUR(:,:,it(i)).*Rinv);
ndivu2(geommask)=0;
ndivu(:,:,i,2)=ndivu2;
ndivu3=N(:,:,i).*(M.fluidUZ.der(:,:,it(i),[0 1]));
ndivu3(geommask)=0;
ndivu(:,:,i,3)=ndivu3;
ugradn1= M.fluidUR(:,:,it(i)).*M.N.der(:,:,it(i),[1 0]);
ugradn1(geommask)=0;
ugradn(:,:,i,1)=ugradn1;
ugradn2= M.fluidUZ(:,:,it(i)).*M.N.der(:,:,it(i),[0 1]);
ugradn2(geommask)=0;
ugradn(:,:,i,2)=ugradn2;
if M.neutcol.present && ~isempty(M.neutcol.io_cross_sec)
% average kinetic energy in each direction
% Er=M.fluidEkin(1,:,:,it(i));
% Ethet=M.fluidEkin(2,:,:,it(i));
% Ez=M.fluidEkin(3,:,:,it(i));
% Ek=squeeze((Er+Ethet+Ez)/M.qe);
% % % corresponding particle velocity
% U=squeeze(sqrt(2*M.weight/M.msim*(Er+Ethet+Ez)));
U=-M.Er(:,:,it(i))./M.Bz';
Ek=0.5*M.msim/M.weight/M.qe*U.^2;
Ek(N(:,:,i)<=0)=0;
% ionisation cross section per cell
sigio=M.sigio(Ek);
Si=M.neutcol.neutdens.*N(:,:,i).*sigio.*U;
end
if M.maxwellsrce.present
Si=Si+M.maxwellsrce.frequency*M.weight/(pi*diff(M.maxwellsrce.zlim)*(M.maxwellsrce.rlim(2)^2-M.maxwellsrce.rlim(1)^2));
end
Si(geommask)=0;
S(:,:,i)=Si;
end
Continuity.N=N;
Continuity.dndt=dndt;
Continuity.ndivu=ndivu;
Continuity.ugradn=ugradn;
Continuity.maxdens=maxdens;
Continuity.it=it;
Continuity.S=S;
end
c4science · Help