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

fennecs/matlab/espic2dhdf5.m
d8035b6d442ageom_testmumps

Last Change
Raw File

espic2dhdf5.m
View Options

classdef espic2dhdf5
%UNTITLED Summary of this class goes here
% Detailed explanation goes here
properties
filename
name
folder
fullpath
timestamp
info
t0d
t1d
t2d
tpart
it0
it1
it2
%% Physical constants
vlight=299792458;
qe=1.60217662E-19;
me=9.109383E-31;
eps_0=8.85418781762E-12;
kb=1.38064852E-23;
%% Run parameters
dt
nrun
nlres
nlsave
nlclassical
nlPhis
nz
nnr
lz
nplasma
potinn
potout
B0
Rcurv
width
n0
temp
femorder
ngauss
plasmadim
radii
H0
P0
normalized
nbspecies
%% Frequencies
omepe
omece
%% Normalizations
tnorm
rnorm
bnorm
enorm
phinorm
vnorm
%% Grid data
rgrid
zgrid
dz
dr
Vcell
%% Magnetic field
Br
Bz
Athet
B
%% Energies
epot
ekin
etot
etot0
eerr
%% 2D time data
N
fluidUR
fluidUZ
fluidUTHET
pot
Er
Ez
CellVol
Presstens
%% Splines
knotsr
knotsz
%% Particle parameters
weight
qsim
msim
nbparts
partepot
R
Z
Rindex
Zindex
partindex
VR
VZ
VTHET
THET
species
end
methods
function file=file(obj)
file=obj.filename;
end
function obj = espic2dhdf5(filename,readparts)
filedata=dir(filename);
if (isempty(filedata))
error("File: ""%s"" doesn't exist",filename)
end
obj.folder=filedata.folder;
obj.filename=filename;
[~, obj.name, ~] = fileparts(obj.filename);
obj.fullpath=[obj.folder,'/',obj.filename];
obj.timestamp=filedata.date;
if nargin==1
readparts=true;
end
%obj.info=h5info(filename);
obj.dt = h5readatt(obj.fullpath,'/data/input.00/','dt');
obj.nrun = h5readatt(obj.fullpath,'/data/input.00/','nrun');
obj.nlres = strcmp(h5readatt(obj.fullpath,'/data/input.00/','nlres'),'y');
obj.nlsave = strcmp(h5readatt(obj.fullpath,'/data/input.00/','nlsave'),'y');
obj.nlclassical =strcmp(h5readatt(obj.fullpath,'/data/input.00/','nlclassical'),'y');
obj.nlPhis =strcmp(h5readatt(obj.fullpath,'/data/input.00/','nlPhis'),'y');
obj.nz = h5readatt(obj.fullpath,'/data/input.00/','nz');
obj.nnr = h5read(obj.fullpath,'/data/input.00/nnr');
obj.lz = h5read(obj.fullpath,'/data/input.00/lz');
obj.qsim = h5readatt(obj.fullpath,'/data/input.00/','qsim');
obj.msim = h5readatt(obj.fullpath,'/data/input.00/','msim');
try
obj.weight=h5readatt(obj.fullpath,'/data/part/','weight');
catch
obj.weight=obj.msim/obj.me;
end
obj.nplasma = h5readatt(obj.fullpath,'/data/input.00/','nplasma');
obj.potinn = h5readatt(obj.fullpath,'/data/input.00/','potinn');
obj.potout = h5readatt(obj.fullpath,'/data/input.00/','potout');
obj.B0 = h5readatt(obj.fullpath,'/data/input.00/','B0');
obj.Rcurv = h5readatt(obj.fullpath,'/data/input.00/','Rcurv');
obj.width = h5readatt(obj.fullpath,'/data/input.00/','width');
obj.n0 = h5readatt(obj.fullpath,'/data/input.00/','n0');
obj.temp = h5readatt(obj.fullpath,'/data/input.00/','temp');
try
obj.it0 = h5readatt(obj.fullpath,'/data/input.00/','it0d');
obj.it1 = h5readatt(obj.fullpath,'/data/input.00/','it2d');
obj.it2 = h5readatt(obj.fullpath,'/data/input.00/','itparts');
catch
obj.it0 = h5readatt(obj.fullpath,'/data/input.00/','it0');
obj.it1 = h5readatt(obj.fullpath,'/data/input.00/','it1');
obj.it1 = h5readatt(obj.fullpath,'/data/input.00/','it2');
end
try
obj.nbspecies=h5readatt(obj.fullpath,'/data/input.00/','nbspecies');
obj.normalized=strcmp(h5readatt(obj.fullpath,'/data/input.00/','rawparts'),'y');
catch
obj.nbspecies=1;
obj.normalized=false;
end
obj.omepe=sqrt(abs(obj.n0)*obj.qe^2/(obj.me*obj.eps_0));
obj.omece=obj.qe*obj.B0/obj.me;
obj.nbparts= h5read(obj.fullpath, '/data/var0d/nbparts');
try
obj.H0 = h5read(obj.fullpath,'/data/input.00/H0');
obj.P0 = h5read(obj.fullpath,'/data/input.00/P0');
catch
obj.H0=3.2e-14;
obj.P0=8.66e-25;
end
% Normalizations
obj.tnorm=1/obj.omepe;
obj.rnorm=obj.vlight*obj.tnorm;
obj.bnorm=obj.B0;
obj.enorm=obj.vlight*obj.bnorm;
obj.phinorm=obj.enorm*obj.rnorm;
obj.vnorm=obj.vlight;
% Grid data
obj.rgrid= h5read(obj.fullpath, '/data/var1d/rgrid')*obj.rnorm;
obj.zgrid= h5read(obj.fullpath, '/data/var1d/zgrid')*obj.rnorm;
obj.dz=(obj.zgrid(end)-obj.zgrid(1))/double(obj.nz);
rid=1;
for i=1:length(obj.nnr)
obj.dr(i)=(obj.rgrid(sum(obj.nnr(1:i))+1)-obj.rgrid(rid))/double(obj.nnr(i));
rid=rid+obj.nnr(i);
end
obj.Vcell=(obj.zgrid(2:end)-obj.zgrid(1:end-1))*((obj.rgrid(2:end).^2-obj.rgrid(1:end-1).^2)*pi)';
Br = h5read(obj.fullpath,'/data/fields/Br')*obj.bnorm;
obj.Br= reshape(Br,length(obj.zgrid),length(obj.rgrid));
Bz = h5read(obj.fullpath,'/data/fields/Bz')*obj.bnorm;
obj.Bz= reshape(Bz,length(obj.zgrid),length(obj.rgrid));
try
Atheta = h5read(obj.fullpath,'/data/fields/Athet')*obj.bnorm;
obj.Athet= reshape(Atheta,length(obj.zgrid),length(obj.rgrid));
catch
end
obj.B=sqrt(obj.Bz.^2+obj.Br.^2);
clear Br Bz
try
obj.t0d=h5read(obj.fullpath,'/data/var0d/time');
catch
obj.t0d=obj.dt.*double(0:length(obj.epot)-1);
end
obj.femorder = h5read(obj.fullpath,'/data/input.00/femorder');
obj.ngauss = h5read(obj.fullpath,'/data/input.00/ngauss');
obj.plasmadim = h5read(obj.fullpath,'/data/input.00/plasmadim');
obj.radii = h5read(obj.fullpath,'/data/input.00/radii');
obj.epot = h5read(obj.fullpath,'/data/var0d/epot');
obj.ekin = h5read(obj.fullpath,'/data/var0d/ekin');
obj.etot = h5read(obj.fullpath,'/data/var0d/etot');
try
obj.etot0 = h5read(obj.fullpath,'/data/var0d/etot0');
obj.eerr = obj.etot-obj.etot0;
catch
obj.eerr = obj.etot-obj.etot(2);
end
if(obj.normalized)
obj.pot=gridquantity(obj.fullpath,'/data/fields/pot',sum(obj.nnr)+1, obj.nz+1,1);
obj.Er=gridquantity(obj.fullpath,'/data/fields/Er',sum(obj.nnr)+1, obj.nz+1,1);
obj.Ez=gridquantity(obj.fullpath,'/data/fields/Ez',sum(obj.nnr)+1, obj.nz+1,1);
else
obj.pot=gridquantity(obj.fullpath,'/data/fields/pot',sum(obj.nnr)+1, obj.nz+1,obj.phinorm);
obj.Er=gridquantity(obj.fullpath,'/data/fields/Er',sum(obj.nnr)+1, obj.nz+1,obj.enorm);
obj.Ez=gridquantity(obj.fullpath,'/data/fields/Ez',sum(obj.nnr)+1, obj.nz+1,obj.enorm);
end
try
obj.t2d = h5read(obj.fullpath,'/data/fields/time');
catch
info=h5info(obj.fullpath,'/data/fields/partdensity');
obj.t2d=obj.dt*(0:info.objspace.Size(2)-1)*double(obj.it1);
end
try
info=h5info(obj.fullpath,'/data/fields/moments');
obj.femorder = h5read(obj.fullpath,'/data/input.00/femorder');
kr=obj.femorder(2)+1;
obj.knotsr=augknt(obj.rgrid,kr);
kz=obj.femorder(1)+1;
obj.knotsz=augknt(obj.zgrid,kz);
try
obj.CellVol= reshape(h5read(obj.fullpath,'/data/fields/volume'),length(obj.knotsz)-kz,length(obj.knotsr)-kr);
obj.CellVol=permute(obj.CellVol,[2,1,3])*obj.rnorm^3;
catch
zvol=fnder(spmak(obj.knotsz,ones(1,length(obj.knotsz)-kz)), -1 );
rvol=fnder(spmak(obj.knotsr,2*pi*[obj.rgrid' 2*obj.rgrid(end)-obj.rgrid(end-1)]), -1 );
ZVol=diff(fnval(zvol,obj.knotsz));
RVol=diff(fnval(rvol,obj.knotsr));
obj.CellVol=RVol(3:end-1)*ZVol(3:end-1)';
obj.CellVol=padarray(obj.CellVol,[1,1],'replicate','post');
end
if(obj.normalized)
obj.N=splinedensity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, 1, 1);
else
obj.N=splinedensity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, abs(obj.qsim/obj.qe), 1);
end
obj.fluidUR=splinevelocity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.vnorm, 2);
obj.fluidUTHET=splinevelocity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.vnorm, 3);
obj.fluidUZ=splinevelocity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.vnorm, 4);
if(obj.normalized)
obj.Presstens=splinepressure(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, obj.vnorm^2*obj.me);
else
obj.Presstens=splinepressure(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, obj.vnorm^2*obj.msim);
end
catch
obj.CellVol=(obj.zgrid(2:end)-obj.zgrid(1:end-1))*((obj.rgrid(2:end).^2-obj.rgrid(1:end-1).^2)*pi)';
obj.CellVol=obj.CellVol';
obj.N=griddensity(obj.fullpath, '/data/fields/partdensity', sum(obj.nnr)+1, obj.nz+1, obj.CellVol, abs(obj.qsim/obj.qe), true);
obj.fluidUR=gridquantity(obj.fullpath, '/data/fields/fluidur', sum(obj.nnr)+1, obj.nz+1, obj.vnorm, true);
obj.fluidUTHET=gridquantity(obj.fullpath, '/data/fields/fluiduthet', sum(obj.nnr)+1, obj.nz+1, obj.vnorm, true);
obj.fluidUZ=gridquantity(obj.fullpath, '/data/fields/fluiduz', sum(obj.nnr)+1, obj.nz+1, obj.vnorm, true);
end
if(readparts)
if(obj.normalized)
obj.R = h5partsquantity(obj.fullpath,'/data/part','R');
obj.Z = h5partsquantity(obj.fullpath,'/data/part','Z');
else
obj.R = h5partsquantity(obj.fullpath,'/data/part','R',obj.rnorm);
obj.Z = h5partsquantity(obj.fullpath,'/data/part','Z',obj.rnorm);
end
try
obj.THET = h5partsquantity(obj.fullpath,'/data/part','THET');
catch
clear obj.THET
end
try
obj.Rindex=h5partsquantity(obj.fullpath,'/data/part','Rindex');
obj.Zindex=h5partsquantity(obj.fullpath,'/data/part','Zindex');
catch
clear obj.Rindex obj.Zindex
end
obj.VR = h5partsquantity(obj.fullpath,'/data/part','UR',obj.vnorm);
obj.VZ = h5partsquantity(obj.fullpath,'/data/part','UZ',obj.vnorm);
obj.VTHET= h5partsquantity(obj.fullpath,'/data/part','UTHET',obj.vnorm);
if(obj.normalized)
obj.partepot = h5partsquantity(obj.fullpath,'/data/part','pot',sign(obj.qsim)*obj.qe);
else
obj.partepot = h5partsquantity(obj.fullpath,'/data/part','pot',sign(obj.qsim)*obj.qe*obj.phinorm);
end
try
obj.partindex = h5partsquantity(obj.fullpath,'/data/part/','partindex');
catch
end
end
try
obj.tpart = h5read(obj.fullpath,'/data/part/time');
catch
obj.tpart=obj.dt*(0:size(obj.R,2)-1)*double(obj.it2);
end
if(obj.nbspecies >1)
obj.species=h5parts.empty;
for i=2:obj.nbspecies
nbparts=h5read(obj.fullpath,sprintf('%s/Nparts',sprintf('/data/part/%2d',i)));
if (nbparts(1)>0)
obj.species(i-1)=h5parts(obj.fullpath,sprintf('/data/part/%2d',i),obj);
end
end
end
end
function Atheta=Atheta(obj,R,Z)
Atheta=0.5*obj.B0*(R-obj.width/pi*(obj.Rcurv-1)/(obj.Rcurv+1)...
.*besseli(1,2*pi*R/obj.width).*cos(2*pi*Z/obj.width));
end
function quantity=H(obj,indices)
%% computes the total energy for the main specie particle indices{1} at time indices{2}
if indices{1}==':'
p=1:obj.VR.nparts;
else
p=indices{1};
end
if indices{2}==':'
t=1:length(obj.tpart);
else
t=indices{2};
end
if size(indices,1)>2
track=indices{3};
else
track=false;
end
quantity=0.5*obj.me*(obj.VR(p,t,track).^2+obj.VTHET(p,t,track).^2+obj.VZ(p,t,track).^2)+obj.partepot(p,t,track);
end
function quantity=P(obj,indices)
%% computes the canonical angular momentum for the main specie particle indices{1} at time indices{2}
if indices{1}==':'
p=1:obj.R.nparts;
else
p=indices{1};
end
if indices{2}==':'
t=1:length(obj.tpart);
else
t=indices{2};
end
if size(indices,1)>2
track=indices{3};
else
track=false;
end
quantity=obj.R(p,t,track).*(obj.VTHET(p,t,track)*obj.me+sign(obj.qsim)*obj.qe*obj.Atheta(obj.R(p,t,track),obj.Z(p,t,track)));
end
function displaypsi(obj,deltat)
%% plot the initial and final radial profile at position z=0 and show the normalized enveloppe function Psi
% relevant for Davidson annular distribution function
f=figure('Name', sprintf('%s Psi',obj.name));
zpos=floor(length(obj.zgrid)/2);
tinit=1;
tend=length(obj.t2d);
deltat=cell2mat(deltat);
if(obj.R.nt<2)
H0=obj.H0;
P0=obj.P0;
else
H0=mean(H(obj,{1:obj.VR.nparts,obj.VR.nt,false}));
P0=mean(P(obj,{1:obj.VR.nparts,obj.VR.nt,false}));
end
lw=1.5;
Mirrorratio=(obj.Rcurv-1)/(obj.Rcurv+1);
locpot=mean(obj.pot(:,zpos,tend-deltat:tend),3);
psi=1+obj.qe*locpot(:)/H0-1/(2*obj.me*H0)*(P0./obj.rgrid+obj.qe*0.5*obj.B0.*(obj.rgrid-obj.width/pi*Mirrorratio*cos(2*pi*obj.zgrid(zpos)/obj.width)*besseli(1,2*pi*obj.rgrid/obj.width))).^2;
locdens=mean(obj.N(:,zpos,tend-deltat:tend),3);
[maxn,In]=max(locdens);%M.N(:,zpos,tinit));
plot(obj.rgrid,obj.N(:,zpos,tinit),'bx-','DisplayName',sprintf('t=%1.2f[ns]',obj.t2d(tinit)*1e9),'linewidth',lw)
hold on
plot(obj.rgrid,locdens,'rx-','DisplayName',sprintf('t=[%1.2f-%1.2f] [ns] averaged',obj.t2d(tend-deltat)*1e9,obj.t2d(tend)*1e9),'linewidth',lw)
plot(obj.rgrid(In-2:end),1./obj.rgrid(In-2:end)*maxn*obj.rgrid(In),'DisplayName','N=c*1/r','linewidth',lw)
xlabel('r [m]')
ylabel('n_e [m^{-3}]')
I=find(psi>0);
if (length(I)>1)
I=[I(1)-2; I(1)-1; I; I(end)+1; I(end)+2];
else
I=obj.nnr(1):length(psi);
end
rq=linspace(obj.rgrid(max(I(1),1)),obj.rgrid(I(end)),500);
psiinterp=interp1(obj.rgrid(I),psi(I),rq,'pchip');
zeroindices=find(diff(psiinterp>=0),2);
maxpsiinterp=max(psiinterp);
plot(rq,maxn*psiinterp/abs(maxpsiinterp),'Displayname','normalized \Psi [a.u.]','linewidth',lw)
ylim([0 inf])
for i=1:length(zeroindices)
border=plot([rq(zeroindices(i)) rq(zeroindices(i))],[0 obj.rgrid(In)/obj.rgrid(In-2)*maxn],'k--','linewidth',lw);
set(get(get(border,'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
end
legend
xlim([0 0.02])
grid on
title(sprintf('Radial density profile at z=%1.2e[m]',obj.zgrid(zpos)))
obj.savegraph(f,sprintf('%sPsi',obj.name),[15 10]);
end
function sref = subsref(obj,s)
% obj(i) is equivalent to obj.Data(i)
switch s(1).type
case '.'
if(strcmp(s(1).subs,'H'))
sref=H(obj,s(2).subs);
elseif(strcmp(s(1).subs,'P'))
sref=P(obj,s(2).subs);
elseif(strcmp(s(1).subs,'plotEnergy'))
plotEnergy(obj);
elseif(strcmp(s(1).subs,'trappedParticles'))
trappedParticles(obj);
elseif(strcmp(s(1).subs,'displayHP'))
displayHP(obj,s(2).subs);
elseif(strcmp(s(1).subs,'displaypsi'))
displaypsi(obj,s(2).subs);
elseif(strcmp(s(1).subs,'Wells'))
Wells(obj,s(2).subs{1},s(2).subs{2},s(2).subs{3});
elseif(strcmp(s(1).subs,'displayaveragetemp'))
displayaveragetemp(obj);
elseif(strcmp(s(1).subs,'displayCurrentsevol'))
displayCurrentsevol(obj,s(2).subs{1});
else
sref=builtin('subsref',obj,s);
end
case '()'
sref=subsref(obj(s(1).subs{1}),s(2:end));
case '{}'
error('MYDataClass:subsref',...
'Not a supported subscripted reference')
end
end
function changed=ischanged(obj)
%ischanged Check if the file has been changed and some data
%must be reloaded
try
filedata=dir(obj.fullpath);
checkedtimestamp=filedata.date;
if (max(checkedtimestamp > obj.timestamp) )
changed=true;
return
end
changed=false;
return
catch
changed=true;
return
end
end
function plotEnergy(obj)
%% Plot the time evolution of the system energy and number of simulated macro particles
tmin=2;
tmax=length(obj.ekin);
f=figure('Name', sprintf('%s Energy',obj.name));
subplot(2,1,1)
plot(obj.t0d(tmin:tmax),obj.ekin(tmin:tmax),'o-',...
obj.t0d(tmin:tmax),obj.epot(tmin:tmax),'d-',...
obj.t0d(tmin:tmax),obj.etot(tmin:tmax),'h-',...
obj.t0d(tmin:tmax),obj.eerr(tmin:tmax),'x--')
legend('ekin', 'epot', 'etot','eerr')
xlabel('Time [s]')
ylabel('Energies [J]')
grid on
xlimits=xlim();
subplot(2,1,2)
try
semilogy(obj.t0d(tmin:tmax),abs(obj.eerr(tmin:tmax)./obj.etot0(tmin:tmax)),'h-')
catch
semilogy(obj.t0d(tmin:tmax),abs(obj.eerr(tmin:tmax)/obj.etot(2)),'h-')
end
xlabel('t [s]')
ylabel('E_{err}/E_{tot}')
xlim(xlimits)
grid on
try
yyaxis right
plot(obj.t0d(tmin:tmax),abs(obj.nbparts(tmin:tmax)./obj.nbparts(1)*100),'d--')
ylabel('Nparts %')
%ylim([0 110])
catch
end
obj.savegraph(f,sprintf('%s/%sEnergy',obj.folder,obj.name));
end
function trappedParticles(obj)
%% Plot the time evolution of the number of physical trapped particles in the main specie
f=figure('Name', sprintf('%s Trapped particles',obj.name));
plot(obj.t0d,obj.nbparts*obj.weight)
xlabel('t [s]')
ylabel('N particles')
obj.savegraph(f,sprintf('%s/%sntrapped',obj.folder,obj.name));
end
function savegraph(obj, figure, name, papsize)
%% Saves the given figure as a pdf and a .fig in the
figure.PaperUnits='centimeters';
if (nargin < 4)
papsize=[14 16];
end
figure.PaperSize=papsize;
print(figure,name,'-dpdf','-fillpage')
savefig(figure,name)
end
function displayHP(obj,tstart)
%% Plot the histogramm of the total energy and canonical angular momentum at the beginning and
% end of the simulation over the full simulation space for the main specie
if(iscell(tstart))
tstart=cell2mat(tstart);
end
if(obj.R.nt>=2)
f=figure('Name', sprintf('%s HP',obj.name));
legtext=sprintf("t=%2.1f - %2.1f [ns]",obj.tpart(tstart)*1e9,obj.tpart(end)*1e9);
subplot(1,2,1)
partsmax=min(obj.nbparts(end),obj.R.nparts);
Hloc=H(obj,{1:obj.nbparts(1),1,false});
h1=histogram(Hloc,20,'BinLimits',[min(Hloc(:)) max(Hloc(:))],'DisplayName',sprintf("t=%2.3d [ns]",obj.tpart(1)*1e9));
hold on
Hloc=mean(H(obj,{1:partsmax,tstart:obj.R.nt,false}),2);
%,'Binwidth',h1.BinWidth
h1=histogram(Hloc,20,'BinLimits',[min(Hloc(:)) max(Hloc(:))],'DisplayName',legtext);
ylabel('counts')
xlabel('H [J]')
legend
subplot(1,2,2)
Ploc=P(obj,{1:obj.nbparts(1),1,false});
h2=histogram(Ploc,50,'BinLimits',[min(Ploc(:)) max(Ploc(:))],'DisplayName',sprintf("t=%2.3d [ns]",obj.tpart(1)*1e9));
hold on
Ploc=mean(P(obj,{1:partsmax,tstart:obj.R.nt,false}),2);
histogram(Ploc,50,'BinLimits',[min(Ploc(:)) max(Ploc(:))],'DisplayName',legtext);
ylabel('counts')
xlabel('P [kg\cdotm^2\cdots^{-1}]')
%clear P
%clear H
legend
%xlim([0.95*h2.BinLimits(1) 1.05*h2.BinLimits(2)])
obj.savegraph(f,sprintf('%s/%sParts_HP',obj.folder,obj.name));
end
end
function Wells(obj,pot,fieldtime,maxdens)
%% Plot the effective potential along the magnetic field lines given the potential pot at time fieldtime
rpos=zeros(size(pot));
[r,~]=meshgrid(obj.rgrid,obj.zgrid);
rathet=(obj.Athet.*r)';
poteff=zeros(size(pot));
potdepth=zeros(size(pot));
for j=1:size(rpos,2)
for i=1:size(rpos,1)
if(isempty(find(rathet(i,j)>rathet(:,1),1,'last')))
rpos(i,j)=1;
rinterp=0;
else
rpos(i,j)=find(rathet(i,j)>rathet(:,1),1,'last');
rinterp=(rathet(i,j)-rathet(rpos(i,j),1))./(rathet(rpos(i,j)+1,1)-rathet(rpos(i,j),1));
end
bmin=obj.B(1,rpos(i,j))+rinterp*(obj.B(1,rpos(i,j)+1)-obj.B(1,rpos(i,j)));
potmin=pot(rpos(i,j),1)+rinterp*(pot(rpos(i,j)+1,1)-pot(rpos(i,j),1));
poteff(i,j)=-(pot(i,j)-potmin)*bmin/(obj.B(j,i)-bmin);
end
end
zzero=ceil(length(obj.zgrid)/2);
for j=1:size(rpos,2)
for i=1:size(rpos,1)
if(isempty(find(rathet(i,j)<rathet(:,zzero),1,'first')))
rpos(i,j)=length(obj.rgrid);
else
rpos(i,j)=max(find(rathet(i,j)<rathet(:,zzero),1,'first'),2);
end
rinterp=(rathet(i,j)-rathet(rpos(i,j)-1,zzero))./(rathet(rpos(i,j),zzero)-rathet(rpos(i,j)-1,zzero));
potmin=pot(rpos(i,j)-1,zzero)+rinterp*(pot(rpos(i,j),zzero)-pot(rpos(i,j)-1,zzero));
potdepth(i,j)=pot(i,j)-potmin;
end
end
f=figure('Name', sprintf('%s Wells',obj.name));
ax1=gca;
h=surface(ax1,obj.zgrid(2:end-1),obj.rgrid(2:end-1),poteff(2:end-1,2:end-1),'edgecolor','none','facecolor','interp');
xlim(ax1,[obj.zgrid(1) obj.zgrid(end)])
ylim(ax1,[obj.rgrid(1) obj.rgrid(end)])
legend
xlabel(ax1,'z [m]')
ylabel(ax1,'r [m]')
if size(fieldtime)<2
title(ax1,sprintf('Ion Well t=%1.2g s n_e=%1.2gm^{-3}',obj.t2d(fieldtime),double(maxdens)))
figurename=sprintf('%sWell%d',obj.name,fieldtime);
else
title(ax1,sprintf('Ion Well t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',obj.t2d(fieldtime(1)),obj.t2d(fieldtime(2)),double(maxdens)))
figurename=sprintf('%sWell%d_%d',obj.name,fieldtime);
end
c = colorbar(ax1);
c.Label.String= 'Effective Potential [eV]';
view(ax1,2)
%set(h,'edgecolor','none');
grid on;
obj.savegraph(f,sprintf('%s/%s',obj.folder,figurename),[16 14]);
f=figure('Name', sprintf('%s Well depth',obj.name));
ax1=gca;
h=surface(ax1,obj.zgrid(2:end-1),obj.rgrid(2:end-1),potdepth(2:end-1,2:end-1),'edgecolor','none','facecolor','interp');
xlim(ax1,[obj.zgrid(1) obj.zgrid(end)])
ylim(ax1,[obj.rgrid(1) obj.rgrid(end)])
legend
xlabel(ax1,'z [m]')
ylabel(ax1,'r [m]')
if size(fieldtime)<2
title(ax1,sprintf('Potential Well t=%1.2g s n_e=%1.2gm^{-3}',obj.t2d(fieldtime),double(maxdens)))
figurename=sprintf('%sPotWell%d',obj.name,fieldtime);
else
title(ax1,sprintf('Potential Well t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',obj.t2d(fieldtime(1)),obj.t2d(fieldtime(2)),double(maxdens)))
figurename=sprintf('%sPotWell%d_%d',obj.name,fieldtime);
end
c = colorbar(ax1);
c.Label.String= 'Potential [eV]';
view(ax1,2)
%set(h,'edgecolor','none');
grid on;
obj.savegraph(f,sprintf('%s/%s',obj.folder,figurename),[16 14]);
end
function displayaveragetemp(obj)
f=figure('Name',sprintf('%s potinn=%f part temperature',obj.name,obj.potinn*obj.phinorm));
vr2=obj.VR(:,:);
vr2=mean(vr2.^2,1)-mean(vr2,1).^2;
vz2=obj.VZ(:,:);
vz2=mean(vz2.^2,1)-mean(vz2,1).^2;
vthet2=obj.VTHET(:,:);
vthet2=mean(vthet2.^2,1)-mean(vthet2,1).^2;
plot(obj.tpart,0.5*obj.me*vr2/obj.qe,'displayname','T_r')
hold on
plot(obj.tpart,0.5*obj.me*vz2/obj.qe,'displayname','T_z')
plot(obj.tpart,0.5*obj.me*vthet2/obj.qe,'displayname','T_{thet}')
xlabel('time [s]')
ylabel('T [eV]')
title(sprintf('\\phi_a=%.1f kV \\phi_b=%.1f kV R=%.1f',obj.potinn*obj.phinorm/1e3,obj.potout*obj.phinorm/1e3,obj.Rcurv))
legend
grid
obj.savegraph(f,sprintf('%s/%s_partstemp',obj.folder,obj.name));
end
function Gamma=Axialflux(obj,timestep,zpos)
Gamma=obj.fluidUZ(:,zpos,timestep).*obj.N(:,zpos,timestep);
end
function I=OutCurrents(obj,timestep)
I=zeros(2,length(timestep));
for j=1:length(timestep)
I(1,j)=trapz(obj.rgrid,obj.Axialflux(timestep(j),1)*2*pi)*obj.qsim/obj.weight;
I(2,j)=trapz(obj.rgrid,obj.Axialflux(timestep(j),obj.nz)*2*pi)*obj.qsim/obj.weight;
end
end
function displayCurrentsevol(obj,timesteps)
if nargin<2
timesteps=1:length(obj.t2d);
end
currents=obj.OutCurrents(timesteps);
f=figure('Name',sprintf('%s Currents',obj.name));
plot(obj.t2d(timesteps),currents(1,:),'Displayname',sprintf('z=%.2f cm',obj.zgrid(1)*100));
hold on
plot(obj.t2d(timesteps),currents(2,:),'Displayname',sprintf('z=%.2f cm',obj.zgrid(end)*100));
legend('location','east')
xlabel('time [s]')
ylabel('I [A]')
title(sprintf('\\phi_b-\\phi_a=%.2g kV, R=%.1f',(obj.potout-obj.potinn)*obj.phinorm/1e3,obj.Rcurv))
obj.savegraph(f,sprintf('%s/%s_outCurrents',obj.folder,obj.name));
end
end
end
c4science · Help