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

fennecs/matlab/@espic2dhdf5/displaytotcurrevol_geom.m
f39b0ef38640master

Last Change
Raw File

displaytotcurrevol_geom.m
View Options

function displaytotcurrevol_geom(obj,timesteps,toptitle,scalet,dens,subdiv,nmean)
% Computes and display the time evolution of the outgoing
% currents at time obj.t2d(timesteps)
%scalet=true scales the time by the ellastic collision
%frequency
%dens = true plot the time evolution of the maximum electron
%density in the simulation domain otherwise plot the total
%number of electrons in the domain
% also plot in a subplot the color coded boundary corresponding
% to each current
if nargin<2
timesteps=1:length(obj.t2d);
end
if nargin<3
scalet=true;
end
if nargin <4
dens=true;
end
if nargin<5
subdiv=1;
end
if nargin<6
nmean=1;
end
if scalet
if obj.neutcol.present
vexb0=(obj.Ezxt(:,:,1).*obj.Br'-obj.Erxt(:,:,1).*obj.Bz')./(obj.B'.^2);
vexb0(obj.geomweight(:,:,1)<=0)=0;
potwell=obj.PotentialWell(0);
vexb0(isnan(potwell))=NaN;
vexb0=mean(abs(vexb0(:)),'omitnan');
E=0.5*obj.msim/obj.weight*vexb0^2/obj.qe;
taucol=1/(obj.neutcol.neutdens*vexb0*(obj.sigio(E)+obj.sigmela(E)+obj.sigmio(E)));
try
Sio_S=1e17*(obj.neutcol.neutdens*vexb0*obj.sigio(E))/(obj.maxwellsrce.frequency*obj.weight/(pi*(obj.maxwellsrce.rlim(2)^2-obj.maxwellsrce.rlim(1)^2)*diff(obj.maxwellsrce.zlim)))
catch
end
tlabel='t/\tau_d [-]';
else
taucol=2*pi/obj.omece;
tlabel='t/\tau_ce [-]';
end
else
taucol=1e-9;
tlabel='t [ns]';
end
if dens
N=obj.N(:,:,timesteps);
geomw=obj.geomweight(:,:,1);
geomw(geomw<0)=0;
geomw(geomw>0)=1;
N=N.*geomw;
%[~,idl]=max(N(:,:,end),[],'all','linear');
%[ir,iz]=ind2sub(size(geomw),idl);
%nmax=squeeze(max(max(N,[],1),[],2));
tn=(obj.t2d(timesteps));
nrhalf=find(obj.rgrid>0.07,1,'first');
if(nrhalf<length(obj.rgrid))
nmax=zeros(2,length(tn));
nmax(1,:)=squeeze(max(max(N(1:nrhalf,:,:),[],1),[],2));
nmax(2,:)=squeeze(max(max(N(nrhalf+1:end,:,:),[],1),[],2));
else
nmax(:)=squeeze(max(max(N(:,:,:),[],1),[],2));
end
nlabel='n_{e,max} [m^{-3}]';
ndlabel='n_{e,max}';
else
t0dst=find(obj.t0d>=obj.t2d(timesteps(1)),1,'first');
t0dlst=find(obj.t0d<=obj.t2d(timesteps(end)),1,'last');
tn=obj.t0d(t0dst:t0dlst);
nmax=obj.npart(t0dst:t0dlst)*obj.weight;
nlabel='Nb e^-';
ndlabel='Nb e^-';
end
[currents,pos]=obj.OutCurrents(timesteps,subdiv);
P=1;
if obj.neutcol.present
P=obj.neutcol.neutdens*obj.kb*300/100;% pressure at room temperature in mbar
end
currents=currents/P;
f=figure('Name',sprintf('%s Charges',obj.name));
tiledlayout(2,1)
nexttile
% Plot the evolution of nb of particles
yyaxis right
if size(nmax,2)>1
nmax=nmax';
end
for i=1:size(nmax,2)
p=plot(tn/taucol,nmax(:,i),'b','linewidth',1.8,'Displayname',sprintf('%s, %d',ndlabel,i));
hold on
end
ylabel(nlabel)
axl=gca;
axl.YAxis(2).Color=p.Color;
ylim([0 inf])
if(obj.B(1,1)>obj.B(end,1))
lname='HFS';
rname='LFS';
else
lname='LFS';
rname='HFS';
end
yyaxis( 'left');
map=colormap(lines);
set(axl,'linestyleorder',{'-',':','--','*','+'},...
'ColorOrder',map(2:7,:), 'NextPlot','replacechildren')
p(1)=plot(axl,obj.t2d(timesteps)/taucol,movmean(currents(1,:),nmean),'Displayname',lname,'linewidth',1.8);
hold on
p(2)=plot(axl,obj.t2d(timesteps)/taucol,movmean(currents(2,:),nmean),'Displayname',rname,'linewidth',1.8);
% Plot the currents
for i=3:size(currents,1)
p(i)=plot(axl,obj.t2d(timesteps)/taucol,movmean(currents(i,:),nmean),'Displayname',sprintf('border %i',i-2),'linewidth',1.8);
end
plot(axl,obj.t2d(timesteps)/taucol,movmean(sum(currents(:,:),1,'omitnan'),nmean),'k-','Displayname','total','linewidth',1.8);
xlabel(tlabel)
if obj.neutcol.present
ylabel('I/p_n [A/mbar]')
else
ylabel('I [A]')
end
grid on
set(gca,'fontsize',12)
ax.YAxis(1).Color='black';
%legend('Orientation','horizontal','location','north','numcolumns',3)
if ~isempty(toptitle)
title(toptitle)
end
ax2=nexttile;
geomw=obj.geomweight(:,:,1);
geomw(geomw<0)=0;
geomw(geomw>0)=NaN;
[c1,hContour]=contourf(ax2,obj.zgrid*1000,obj.rgrid*1000,geomw, [0 0]);
hold on
drawnow;
grid on;
for i=1:length(pos)
plot(ax2,pos{i}(1,:)*1000,pos{i}(2,:)*1000,'linestyle',p(i+2).LineStyle,...
'color',p(i+2).Color,'marker',p(i+2).Marker,...
'displayname',sprintf('border %i',i),'linewidth',1.8)
hold on
end
title('Domain')
plot(ax2,ones(size(obj.rgrid))*obj.zgrid(1)*1000,obj.rgrid*1000,'linestyle',p(1).LineStyle,...
'color',p(1).Color,'marker',p(1).Marker,...
'displayname',lname,'linewidth',1.8)
plot(ax2,ones(size(obj.rgrid))*obj.zgrid(end)*1000,obj.rgrid*1000,'linestyle',p(2).LineStyle,...
'color',p(2).Color,'marker',p(2).Marker,...
'displayname',rname,'linewidth',1.8)
xlabel('z [mm]')
ylabel('r [mm]')
grid on
set(gca,'fontsize',12)
hFills=hContour.FacePrims;
[hFills.ColorType] = deal('truecoloralpha'); % default = 'truecolor'
try
hFills(1).ColorData = uint8([150;150;150;255]);
for idx = 2 : numel(hFills)
hFills(idx).ColorData(4) = 0; % default=255
end
catch
end
%legend('Orientation','horizontal','location','north','numcolumns',4)
fprintf('mean total current: %f [A/mbar]\n',mean(sum(currents(:,max(1,size(currents,2)-30):end),1,'omitnan')));
%if nargin <3
% sgtitle(sprintf('\\phi_b-\\phi_a=%.2g kV, B=%f T',(obj.potout-obj.potinn)*obj.phinorm/1e3,mean(obj.B(:))))
%elseif ~isempty(toptitle)
% sgtitle(toptitle)
%end
if length(subdiv)>1
obj.savegraph(f,sprintf('%s/%s_totIEvol%i%i_div%i',obj.folder,obj.name,scalet,dens,nmean),[16 14]);
else
obj.savegraph(f,sprintf('%s/%s_totIEvol%i%i_%i',obj.folder,obj.name,scalet,dens,nmean),[16 14]);
end
end
c4science · Help