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

fennecs/matlab/Continuityspline.m
118e68d1a0d5breakdown

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.6;
norm=false;
log=false;
zlims=[-inf inf];
case 3
log=false;
contourscale=0.6;
zlims=[-inf inf];
case 4
contourscale=0.6;
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;
N=M.N(:,:,it);
maxdens=max(N(:));
densitycontour=contour(M.zgrid,M.rgrid,mean(N,3),[contourscale contourscale]*maxdens);
rgridmax=sum(M.nnr(1:2));
%% Continuity
f3=figure();
t=M.t2d(it);
dndt=zeros(size(N,1),size(N,2),size(N,3)-2);
ndivu=zeros(size(N,1),size(N,2),size(N,3)-2,3);
ugradn=zeros(size(N,1),size(N,2),size(N,3)-2,2);
[R,~]=meshgrid(M.rgrid,M.zgrid);
Rinv=1./R';
Rinv(isinf(Rinv))=0;
for i=2:size(N,3)-1
dndt(:,:,i-1)=(N(:,:,i+1)-N(:,:,i-1))./((t(i+1)-t(i-1)));
ndivu(:,:,i-1,1)=N(:,:,i).*(M.fluidUR.der(:,:,it(i),[1 0]));
ndivu(:,:,i-1,2)=N(:,:,i).*(M.fluidUR(:,:,it(i)).*Rinv);
ndivu(:,:,i-1,3)=N(:,:,i).*(M.fluidUZ.der(:,:,it(i),[0 1]));
ugradn(:,:,i-1,1)= M.fluidUR(:,:,it(i)).*M.N.der(:,:,it(i),[1 0]);
ugradn(:,:,i-1,2)= M.fluidUZ(:,:,it(i)).*M.N.der(:,:,it(i),[0 1]);
end
Continuity.N=N(:,:,2:end-1);
Continuity.dndt=dndt;
Continuity.ndivu=ndivu;
Continuity.ugradn=ugradn;
Continuity.maxdens=maxdens;
Continuity.it=it;
n=mean(N,3);
dn=ones(size(n));%n/(t(3)-t(1));
ax7=subplot(4,1,1);
plotvalue(ax7,mean(Continuity.dndt,3)./dn,'\partialn/\partialt','\partialn/\partialt [m^{-3}s^{-1}]',true)
%zlim=[max([min(min(mean(sum(Continuity.ndivu(2:end,:,:),4),3))),min(min(mean(sum(Continuity.ugradn(2:end,:,:),4),3)))]) min([max(max(mean(sum(Continuity.ndivu(2:end,:,:),4),3))),max(max(mean(sum(Continuity.ugradn(2:end,:,:),4),3)))])];
zlim=[-inf inf];
ax8=subplot(4,1,2);
plotvalue(ax8,mean(sum(Continuity.ndivu,4)./dn,3),'n\nabla(u)','n\nabla(u)[m^{-3}s^{-1}]',true,zlim)
ax9=subplot(4,1,3);
plotvalue(ax9,mean(sum(Continuity.ugradn,4)./dn,3),'u\nabla(n)','u\nabla(n)[m^{-3}s^{-1}]',true,zlim)
ax10=subplot(4,1,4);
plotvalue(ax10,mean(sum(Continuity.ndivu,4)+sum(Continuity.ugradn,4)+Continuity.dndt,3)./dn,'Total','total[m^{-3}s^{-1}]',true,zlim)
sgtitle(sprintf('Continuity t=[%1.2g-%1.2g]s n_e=%1.2g m^{-3}',M.t2d(min(it)),M.t2d(max(it)),double(maxdens)))
f3.PaperOrientation='portrait';
[~, name, ~] = fileparts(M.file);
f3.PaperUnits='centimeters';
papsize=[20 19];
f3.PaperSize=papsize;
print(f3,sprintf('%scontinuity',name),'-dpdf','-fillpage')
linkaxes([ax7,ax8,ax9,ax10])
f3=figure;
ax7=subplot(3,1,1);
plotvalue(ax7,squeeze(mean(Continuity.ndivu(:,:,:,1),3))./dn,'ndivu er1','radial 1 ndivu [m^{-3}s^{-1}]',true)
ax8=subplot(3,1,2);
plotvalue(ax8,squeeze(mean(Continuity.ndivu(:,:,:,2),3))./dn,'ndivu er2','radial 2 ndivu [m^{-3}s^{-1}]',true)
ax9=subplot(3,1,3);
plotvalue(ax9,squeeze(mean(Continuity.ndivu(:,:,:,3),3))./dn,'ndivu ez','axial ndivu [m^{-3}s^{-1}]',true)
sgtitle(sprintf('Continuity ndivu t=[%1.2g-%1.2g]s n_e=%1.2g m^{-3}',M.t2d(min(it)),M.t2d(max(it)),double(maxdens)))
f3.PaperOrientation='portrait';
[~, name, ~] = fileparts(M.file);
f3.PaperUnits='centimeters';
papsize=[20 19];
f3.PaperSize=papsize;
print(f3,sprintf('%scontinuityndivu',name),'-dpdf','-fillpage')
linkaxes([ax7,ax8,ax9])
f3=figure;
ax7=subplot(2,1,1);
plotvalue(ax7,squeeze(mean(Continuity.ugradn(:,:,:,1),3))./dn,'ugradn er','radial ugradn [m^{-3}s^{-1}]',true)
ax8=subplot(2,1,2);
plotvalue(ax8,squeeze(mean(Continuity.ugradn(:,:,:,2),3))./dn,'ugradn ez','axial ugradn [m^{-3}s^{-1}]',true)
sgtitle(sprintf('Continuity ugradn t=[%1.2g-%1.2g]s n_e=%1.2g m^{-3}',M.t2d(min(it)),M.t2d(max(it)),double(maxdens)))
f3.PaperOrientation='portrait';
[~, name, ~] = fileparts(M.file);
f3.PaperUnits='centimeters';
papsize=[20 19];
f3.PaperSize=papsize;
print(f3,sprintf('%scontinuityugradn',name),'-dpdf','-fillpage')
linkaxes([ax7,ax8])
function plotvalue(axis,matrix,titl, clab, rm1strow, zlims)
if nargin < 5
rm1strow=false;
end
if nargin>5
zmin=zlims(1);
zmax=zlims(2);
else
zmin=-inf;
zmax=inf;
end
if rm1strow
if log
h=surface(axis,M.zgrid,M.rgrid(2:end),abs(matrix(2:end,:)));
else
h=surface(axis,M.zgrid,M.rgrid(2:end),matrix(2:end,:));
end
else
if log
h=surface(axis,M.zgrid,M.rgrid,abs(matrix));
else
h=surface(axis,M.zgrid,M.rgrid,matrix);
end
end
set(h,'edgecolor','none');
hold on;
if log
zpos=interp2(M.zgrid, M.rgrid, abs(matrix), densitycontour(1,2:end), densitycontour(2,2:end));
else
zpos=interp2(M.zgrid, M.rgrid, matrix, densitycontour(1,2:end), densitycontour(2,2:end));
end
plot3(axis,densitycontour(1,2:end),densitycontour(2,2:end),zpos,'-.','linewidth',2,'color',contourcolor)
xlabel(axis,'z [m]')
ylabel(axis,'r [m]')
xlim(axis,[M.zgrid(1) M.zgrid(end)])
ylim(axis,[M.rgrid(1) M.rgrid(rgridmax)])
colormap(axis,'jet')
c = colorbar(axis);
c.Label.String=clab;
caxis(axis,[zmin zmax])
if norm
caxis(axis,[Fminr Fmaxr]);
end
if log
set(axis,'colorscale','log')
end
title(axis,titl)
grid(axis, 'on')
view(axis,2)
end
end
c4science · Help