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

fennecs/matlab/@espic2dhdf5/espic2dhdf5.m
118e68d1a0d5breakdown

Last Change
Raw File

espic2dhdf5.m
View Options

classdef espic2dhdf5
%espic2dhdf5 General class used to treat hdf5 result files of espic2d code
% A result file is loaded with a call to M=espic2dhdf5(filename) where filename is the relative or absolute file path
% after loading, several quantities and composite diagnostics such as moments of the distribution function or individual particles
% quantities can be accessed.
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 % simulation time step
nrun % number of time steps simulated
nlres
nlsave
nlclassical % Was the equation of motion solved in the classical framework
nlPhis % Was the self-consistent electric field computed
nz % number of intervals in the z direction for the grid
nnr % number of intervals in the r direction for the grid for each of the 3 mesh regions
lz % physical axial dimension of the simulation space
nplasma % Number of initial macro particles
potinn % Normalized electric potential at the coaxial insert
potout % Normalized electric potential at the cylinder surface
B0 % Normalization for the magnetic field
Rcurv % Magnetic mirror ratio
width % Magnetic mirror length
n0 % Initial particle density in case of old particle loading
temp % Initial particle temperature in case of old particle loading
femorder % finite element method order in z and r direction
ngauss % Order of the Gauss integration method for the FEM
plasmadim % initial dimensions of the plasma for the old particle loading system
radii % Radial limits of the three mesh regions coarse,fine,coarse
H0 % Initial particle Energy for Davidsons distribution function
P0 % Initial particle Angular momentum for Davidsons distribution function
normalized % Are the parts quantities normalized in the h5 file
nbspecies % Number of species simulated
%% Frequencies
omepe % Reference plasma frequency used for normalization
omece % Reference cyclotronic frequency for normalization
%% Normalizations
tnorm % Time normalization
rnorm % Dimension normalization
bnorm % Magnetic field normalization
enorm % Electric field normalization
phinorm % Electric potential normalization
vnorm % Velocity normalization
%% Grid data
rgrid % Radial grid position points
zgrid % Axial grid position points
dz % Axial grid step
dr % Radial grid step for the three mesh regions
CellVol % Volume of the cell used for density calculation
%% Magnetic field
Br % Radial magnetic field
Bz % Axial magnetic field
Athet % Azimuthal component of the Magnetic potential vector
rAthet % r*Athet used for the representation of magnetic field lines
B % Magnetic field amplitude
sinthet % ratio to project quantities along the magnetic field lines
costhet % ratio to project quantities along the magnetic field lines
%% Energies
epot % Time evolution of the particles potential energy
ekin % Time evolution of the particles kinetic energy
etot % Time evolution of the particles total energy
etot0 % Time evolution of the reference particle total energy
eerr % Time evolution of the error on the energy conservation
npart % Time evolution of the number of simulated
%% 2D time data evaluated on grid points
N % main specie Density
fluidUR % main specie radial fluid velocity
fluidUZ % main specie axial fluid velocity
fluidUTHET % main specie azimuthal fluid velocity
pot % Electric potential evaluated at grid points
potxt % External Electric potential evaluated at grid points
phi % Electric potential in spline form
Er % Radial electric field
Ez % Axial electric field
Erxt % External Radial electric field
Ezxt % External Axial electric field
Presstens % Pressure tensor
fluidEkin % average kinetic energy in each direction
%% Splines
knotsr % Spline radial knots
knotsz % Spline axial knots
%% Particle parameters
weight % Macro particle numerical weight of the main specie
qsim % Macro particle charge
msim % Macro particle mass
nbparts % Time evolution of the number of simulated particles
partepot % Electric potential at the particles positions
R % Particles radial position
Z % Particles axial position
Rindex % Particles radial grid index
Zindex % Particles axial grid index
partindex % Particles unique id for tracing trajectories
VR % Particles radial velocity
VZ % Particles axial velocity
VTHET % Particles azimuthal velocity
THET % Particles azimuthal position
species % Array containing the other simulated species
%% Celldiag
celldiag % Array containing the cell diagnostic data
nbcelldiag % Total number of cell diagnostics
%% Curvilinear geometry
conformgeom % stores if we use the conforming or nonconforming boundary conditions
r_a
r_b
z_r
z_0
r_0
r_r
L_r
L_z
Interior
above1
above2
interior
walltype
geomweight
%% Maxwell source parameters
maxwellsrce
%% Collision with neutral parameters
neutcol
nudcol % effective momentum collision frequency
end
methods
function file=file(obj)
% returns the h5 file name
file=obj.filename;
end
function obj = espic2dhdf5(filename,readparts,old)
% Reads the new result file filename and read the parts data if readparts==true
% adds the helper_classes folder to the path
matlabfuncpath = dir([mfilename('fullpath'),'.m']);
addpath(sprintf('%s/../helper_classes',matlabfuncpath.folder));
% Try catch are there for compatibility with older simulation files
filedata=dir(filename);
if (isempty(filedata))
error("File: ""%s"" doesn't exist",filename)
end
obj.folder=filedata.folder;
obj.filename=filename;
[~, obj.name, ext] = fileparts(obj.filename);
obj.filename=[obj.name,ext];
obj.fullpath=[obj.folder,'/',obj.filename];
obj.timestamp=filedata.date;
if nargin==1
readparts=true;
end
if nargin<3
old=false;
end
%obj.info=h5info(filename);
%% Read the run parameters
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.r_a=h5readatt(obj.fullpath,'/data/input.00/geometry','r_a');
obj.r_b=h5readatt(obj.fullpath,'/data/input.00/geometry','r_b');
obj.z_r=h5readatt(obj.fullpath,'/data/input.00/geometry','z_r');
obj.r_r=h5readatt(obj.fullpath,'/data/input.00/geometry','r_r');
obj.r_0=h5readatt(obj.fullpath,'/data/input.00/geometry','r_0');
obj.z_0=h5readatt(obj.fullpath,'/data/input.00/geometry','z_0');
obj.above1=h5readatt(obj.fullpath,'/data/input.00/geometry','above1');
obj.above2=h5readatt(obj.fullpath,'/data/input.00/geometry','above2');
obj.interior=h5readatt(obj.fullpath,'/data/input.00/geometry','interior');
obj.walltype=h5readatt(obj.fullpath,'/data/input.00/geometry','walltype');
obj.conformgeom=false;
catch
obj.conformgeom=true;
obj.walltype=0;
obj.r_a=obj.rgrid(1);
obj.r_b=obj.rgrid(end);
obj.above1=1;
obj.above2=-1;
end
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
try
obj.nbspecies=h5readatt(obj.fullpath,'/data/part/','nbspecies');
catch
obj.nbspecies=h5readatt(obj.fullpath,'/data/input.00/','nbspecies');
end
obj.normalized=strcmp(h5readatt(obj.fullpath,'/data/input.00/','rawparts'),'y');
catch
obj.nbspecies=1;
obj.normalized=false;
end
try
obj.nbcelldiag=h5readatt(obj.fullpath,'/data/celldiag/','nbcelldiag');
catch
obj.nbcelldiag=0;
end
obj.omepe=sqrt(abs(obj.n0)*obj.qe^2/(obj.me*obj.eps_0));
obj.omece=obj.qe*obj.B0/obj.me;
obj.npart= h5read(obj.fullpath, '/data/var0d/nbparts');
try
obj.nudcol= h5read(obj.fullpath, '/data/var0d/nudcol');
catch
end
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
if old
obj.tnorm=abs(1/obj.omepe);
else
obj.tnorm=min(abs(1/obj.omepe),abs(1/obj.omece));
end
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
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));
[rmeshgrid,~]=meshgrid(obj.rgrid,obj.zgrid);
obj.rAthet=(rmeshgrid.*obj.Athet)';
catch
end
obj.B=sqrt(obj.Bz.^2+obj.Br.^2);
obj.costhet=(obj.Br./obj.B)';
obj.sinthet=(obj.Bz./obj.B)';
clear Br Bz
try
obj.t0d=h5read(obj.fullpath,'/data/var0d/time');
catch
obj.t0d=obj.dt.*double(0:length(obj.epot)-1);
end
try
obj.Erxt=reshape(h5read(obj.fullpath,'/data/input.00/Erxt'),length(obj.zgrid),length(obj.rgrid))'*obj.enorm;
obj.Ezxt=reshape(h5read(obj.fullpath,'/data/input.00/Ezxt'),length(obj.zgrid),length(obj.rgrid))'*obj.enorm;
obj.potxt=reshape(h5read(obj.fullpath,'/data/input.00/potxt'),length(obj.zgrid),length(obj.rgrid))'*obj.phinorm;
catch
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
try
obj.geomweight = h5read(obj.fullpath,'/data/input.00/geometry/geomweight');
obj.geomweight= reshape(obj.geomweight,length(obj.zgrid),length(obj.rgrid),[]);
obj.geomweight = permute(obj.geomweight,[2,1,3]);
catch
obj.geomweight=ones(length(obj.rgrid),length(obj.zgrid),3);
end
geomweight=ones(length(obj.rgrid),length(obj.zgrid));
if(obj.normalized)
obj.N=splinedensity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, 1, geomweight, 1);
obj.phi=splinequantity(obj.fullpath,'/data/fields/phi', obj.knotsr, obj.knotsz, obj.femorder, 1, obj.geomweight(:,:,1), -1);
else
obj.N=splinedensity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, abs(obj.qsim/obj.qe), geomweight, 1);
end
obj.fluidUR=splinevelocity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.vnorm, geomweight, 2);
obj.fluidUTHET=splinevelocity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.vnorm, geomweight, 3);
obj.fluidUZ=splinevelocity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.vnorm, geomweight, 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, geomweight);
obj.fluidEkin=splineenergy(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, obj.vnorm^2*obj.me*0.5, geomweight);
else
obj.Presstens=splinepressure(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, obj.vnorm^2*obj.msim, geomweight);
obj.fluidEkin=splineenergy(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, obj.vnorm^2*obj.msim*0.5, geomweight);
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 we have a maxwellian source, read its parameters
try
obj.maxwellsrce.rlim=h5read(obj.fullpath, '/data/input.00/maxwellsource/rlimits');
obj.maxwellsrce.zlim=h5read(obj.fullpath, '/data/input.00/maxwellsource/zlimits');
obj.maxwellsrce.frequency=h5readatt(obj.fullpath, '/data/input.00/maxwellsource','frequency');
obj.maxwellsrce.radialtype=h5readatt(obj.fullpath, '/data/input.00/maxwellsource','radialtype');
obj.maxwellsrce.temperature=h5readatt(obj.fullpath, '/data/input.00/maxwellsource','temperature');
obj.maxwellsrce.time_end=h5readatt(obj.fullpath, '/data/input.00/maxwellsource','time_end');
obj.maxwellsrce.time_start=h5readatt(obj.fullpath, '/data/input.00/maxwellsource','time_start');
obj.maxwellsrce.vth=h5readatt(obj.fullpath, '/data/input.00/maxwellsource','vth');
obj.maxwellsrce.present=true;
catch
obj.maxwellsrce.present=false;
end
%% load neutcol parameters
try
obj.neutcol.neutdens=double(h5readatt(obj.fullpath, '/data/input.00/neutcol','neutdens'));
obj.neutcol.neutpressure=double(h5readatt(obj.fullpath, '/data/input.00/neutcol','neutpressure'));
obj.neutcol.scatter_fac=double(h5readatt(obj.fullpath, '/data/input.00/neutcol','scatter_fac'));
obj.neutcol.Eion=double(h5readatt(obj.fullpath, '/data/input.00/neutcol','Eion'));
obj.neutcol.E0=double(h5readatt(obj.fullpath, '/data/input.00/neutcol','E0'));
obj.neutcol.Escale=double(h5readatt(obj.fullpath, '/data/input.00/neutcol','Escale'));
try
obj.neutcol.io_cross_sec=double(h5read(obj.fullpath, '/data/input.00/neutcol/io_cross_sec'));
obj.neutcol.io_cross_sec(:,2)=obj.neutcol.io_cross_sec(:,2)*obj.rnorm^2;
obj.neutcol.io_cross_sec(:,3)=[log(obj.neutcol.io_cross_sec(2:end,2)./obj.neutcol.io_cross_sec(1:end-1,2))...
./log(obj.neutcol.io_cross_sec(2:end,1)./obj.neutcol.io_cross_sec(1:end-1,1)); 0];
catch
obj.neutcol.io_cross_sec=[];
end
try
obj.neutcol.ela_cross_sec=double(h5read(obj.fullpath, '/data/input.00/neutcol/ela_cross_sec'));
obj.neutcol.ela_cross_sec(:,2)=obj.neutcol.ela_cross_sec(:,2)*obj.rnorm^2;
obj.neutcol.ela_cross_sec(:,3)=[log(obj.neutcol.ela_cross_sec(2:end,2)./obj.neutcol.ela_cross_sec(1:end-1,2))...
./log(obj.neutcol.ela_cross_sec(2:end,1)./obj.neutcol.ela_cross_sec(1:end-1,1)); 0];
catch
obj.neutcol.ela_cross_sec=[];
end
obj.neutcol.present=true;
catch
obj.neutcol.present=false;
end
% Read the main particles parameters
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
vscale=obj.vnorm;
obj.VR = h5partsquantity(obj.fullpath,'/data/part','UR',vscale);
obj.VZ = h5partsquantity(obj.fullpath,'/data/part','UZ',vscale);
obj.VTHET= h5partsquantity(obj.fullpath,'/data/part','UTHET',vscale);
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');
obj.nbparts = h5read(obj.fullpath,'/data/part/Nparts');
catch
obj.nbparts=obj.npart;
obj.tpart=obj.dt*(0:size(obj.R,2)-1)*double(obj.it2);
end
if(obj.nbspecies >1)
obj.species=h5parts.empty(obj.nbspecies,0);
for i=2:obj.nbspecies
obj.species(i-1)=h5parts(obj.fullpath,sprintf('/data/part/%2d',i),obj);
end
end
if(obj.nbcelldiag > 0)
obj.celldiag=h5parts.empty;
j=0;
for i=1:obj.nbcelldiag
nbparts=h5read(obj.fullpath,sprintf('%s/Nparts',sprintf('/data/celldiag/%02d',i)));
if (sum(nbparts)>0)
j=j+1;
obj.celldiag(j)=h5parts(obj.fullpath,sprintf('/data/celldiag/%02d',i),obj);
obj.celldiag(j).rindex=double(h5readatt(obj.fullpath, sprintf('/data/celldiag/%02d',i),'rindex'))+(1:2);
obj.celldiag(j).zindex=double(h5readatt(obj.fullpath, sprintf('/data/celldiag/%02d',i),'zindex'))+(1:2);
end
end
end
end
%------------------------------------------
% Functions for accesing secondary simulation quantities
function Atheta=Atheta(obj,R,Z)
% halflz=(obj.zgrid(end)+obj.zgrid(1))/2;
% 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-halflz)/obj.width));
Atheta=interp2(obj.rgrid,obj.zgrid,obj.Athet,R,Z);
end
function quantity=H(obj,indices)
%% computes the total energy for the main specie particle indices{1} at time indices{2}
if strcmp(indices{1},':')
p=1:obj.VR.nparts;
else
p=indices{1};
end
if strcmp(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 strcmp(indices{1},':')
p=1:obj.R.nparts;
else
p=indices{1};
end
if strcmp(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 quantity=Vpar(obj,varargin)
%Vpar Computes the parallel velocity for the main specie particle indices{1} at time indices{2}
if(~iscell(varargin))
indices=mat2cell(varargin);
else
indices=varargin;
end
if strcmp(indices{1},':')
p=1:obj.R.nparts;
else
p=indices{1};
end
if strcmp(indices{2},':')
t=1:length(obj.tpart);
else
t=indices{2};
end
if size(indices,1)>2
track=indices{3};
else
track=false;
end
Zp=obj.Z(p,t,track);
Rp=obj.R(p,t,track);
Bzp=interp2(obj.zgrid,obj.rgrid,obj.Bz',Zp,Rp,'makima');
Brp=interp2(obj.zgrid,obj.rgrid,obj.Br',Zp,Rp,'makima');
Bp=sqrt(Bzp.^2+Brp.^2);
costhet=Bzp./Bp;
sinthet=Brp./Bp;
quantity=obj.VR(p,t,track).*sinthet+obj.VZ(p,t,track).*costhet;
end
function quantity=Vperp(obj,varargin)
%Vperp Computes the perpendicular velocity in the guidind center reference frame,
% for the main specie particle indices{1} at time indices{2}
if(~iscell(varargin))
indices=mat2cell(varargin);
else
indices=varargin;
end
if strcmp(indices{1},':')
p=1:obj.R.nparts;
else
p=indices{1};
end
if strcmp(indices{2},':')
t=1:length(obj.tpart);
else
t=indices{2};
end
if size(indices,2)>2
track=indices{3};
else
track=false;
end
if size(indices,2)>3
gcs=indices{4};
else
gcs=false;
end
Zp=obj.Z(p,t,track);
Rp=obj.R(p,t,track);
Bzp=interp2(obj.zgrid,obj.rgrid,obj.Bz',Zp,Rp,'makima');
Brp=interp2(obj.zgrid,obj.rgrid,obj.Br',Zp,Rp,'makima');
Bp=sqrt(Bzp.^2+Brp.^2);
costhet=Bzp./Bp;
sinthet=Brp./Bp;
Vdrift=zeros(size(Zp));
if gcs
for j=1:length(t)
[~, tfield]=min(abs(obj.t2d-obj.tpart(t(j))));
timeEr=obj.Er(:,:,tfield);
timeEz=obj.Ez(:,:,tfield);
%posindE=sub2ind(size(timeEr),Rind(:,j),Zind(:,j));
timeErp=interp2(obj.zgrid,obj.rgrid,timeEr,Zp(:,j),Rp(:,j));
timeEzp=interp2(obj.zgrid,obj.rgrid,timeEz,Zp(:,j),Rp(:,j));
Vdrift(:,j)=(timeEzp.*Brp(:,j)-timeErp.*Bzp(:,j))./Bp(:,j).^2;
end
end
quantity=sqrt((obj.VTHET(p,t,track)-Vdrift).^2+(obj.VR(p,t,track).*costhet-obj.VZ(p,t,track).*sinthet).^2);
end
function quantity=cyclphase(obj,varargin)
%Vperp Computes the cyclotronic phase for the main specie
if(~iscell(varargin))
indices=mat2cell(varargin);
else
indices=varargin;
end
if strcmp(indices{1},':')
p=1:obj.R.nparts;
else
p=indices{1};
end
if strcmp(indices{2},':')
t=1:length(obj.tpart);
else
t=indices{2};
end
if size(indices,2)>2
track=indices{3};
else
track=false;
end
Zp=obj.Z(p,t,track);
Rp=obj.R(p,t,track);
[~, zind(1)]=min(abs(obj.zgrid-0.005262));
[~, zind(2)]=min(abs(obj.zgrid-0.006637));
[~, rind(1)]=min(abs(obj.rgrid-0.0784));
[~, rind(2)]=min(abs(obj.rgrid-0.07861));
indices=Zp<obj.zgrid(zind(2)) & Zp>=obj.zgrid(zind(1)) &...
Rp<obj.rgrid(rind(2)) & Rp>=obj.rgrid(rind(1));
%Zp=Zp(indices);
%Rp=Rp(indices);
%p=indices;
Bzp=interp2(obj.zgrid,obj.rgrid,obj.Bz',Zp,Rp,'makima');
Brp=interp2(obj.zgrid,obj.rgrid,obj.Br',Zp,Rp,'makima');
Bp=sqrt(Bzp.^2+Brp.^2);
costhet=Bzp./Bp;
sinthet=Brp./Bp;
vr=(obj.VR(p,t,track).*costhet-obj.VZ(p,t,track).*sinthet);
vperp=obj.Vperp(p,t,track,true);
vr=vr(indices);
vperp=vperp(indices);
cospsi=vr./vperp;
quantity=acos(cospsi);
end
function p=borderpoints(obj)
%gw= contourc(obj.zgrid,obj.rgrid,obj.geomweight(:,:,1),[0 0]);
gw=obj.ellipseborder;
p{1}=gw;
if obj.walltype==2
zpos=obj.zgrid(obj.zgrid<(min(gw(1,:))) | obj.zgrid>(max(gw(1,:))));
p{2}=[zpos,obj.rgrid(end)*ones(size(zpos))]';
end
if sum(obj.geomweight(1,:,1)==0)==length(obj.zgrid)
p{3}=[obj.zgrid';obj.rgrid(1)*ones(size(obj.zgrid))'];
end
end
function p=ellipseborder(obj)
z=linspace(-0.998,0.998,1000)*obj.z_r;
p=zeros(4,length(z));
for i=1:length(z)
p(1,i)=z(i)+obj.z_0;
p(2,i)=obj.r_0-obj.r_r*sqrt(1-(z(i)/obj.z_r)^2);
p(3,i)=2/(obj.z_r^2)*(z(i));
p(4,i)=2/(obj.r_r^2)*(p(2,i)-obj.r_0);
end
norm=sqrt(p(3,:).^2+p(4,:).^2);
p(3,:)=double(obj.interior)*p(3,:)./norm;
p(4,:)=double(obj.interior)*p(4,:)./norm;
end
function charge=totcharge(obj,fieldstep)
% Integrates the density profile over the full volume to obtain
% the total number of electrons in the volume
N=splinedensity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder,ones(size(obj.CellVol)), 1, 1);
charge=sum(sum(N(:,:,end)));
end
function Gamma=Axialflux(obj,timestep,zpos)
% Computes the axial particle flux n*Uz at timestep timestep and axial position zpos
Gamma=obj.fluidUZ(:,zpos,timestep).*obj.N(:,zpos,timestep);
end
function Gamma=Metallicflux(obj,timestep)
% Computes the particle flux at timestep timestep on the
% metallic boundaries
% We find the borderpoints
p=obj.borderpoints;
gamma=cell(size(p));
for i=1:length(p)
bp=p{i};
if size(bp,1)==2
% We get the normals at these positions and normalise them
Nr=-interp2(obj.zgrid,obj.rgrid,obj.geomweight(:,:,3),bp(1,:),bp(2,:));
Nz=-interp2(obj.zgrid,obj.rgrid,obj.geomweight(:,:,2),bp(1,:),bp(2,:));
norm=sqrt(Nr.^2+Nz.^2);
Nr=Nr./norm;
Nz=Nz./norm;
else
Nr=bp(4,:);
Nz=bp(3,:);
end
gam=zeros(size(bp,2),length(timestep));
% we get the density and fluid velocities at the desired time
% steps
for j=1:length(timestep)
%n=interp2(obj.zgrid,obj.rgrid,obj.N(:,:,timestep(j)),bp(1,:),bp(2,:));
n=obj.N.posval({bp(2,:),bp(1,:),timestep(j)})';
%uz=interp2(obj.zgrid,obj.rgrid,obj.fluidUZ(:,:,timestep(j)),bp(1,:),bp(2,:));
%ur=interp2(obj.zgrid,obj.rgrid,obj.fluidUR(:,:,timestep(j)),bp(1,:),bp(2,:));
uz=obj.fluidUZ.posval({bp(2,:),bp(1,:),timestep(j)})';
ur=obj.fluidUR.posval({bp(2,:),bp(1,:),timestep(j)})';
gam(:,j)=n.*(ur.*Nr+uz.*Nz);
end
gamma{i}=gam;
end
Gamma.p=p;
Gamma.gamma=gamma;
end
function I=OutCurrents(obj,timestep)
% Computes the Outgoing currens at the simulation axial boundaries at timestep timestep
% This is simply the surface integral of the axial flux
%Iz=zeros(2,length(timestep));
flux=obj.Axialflux(timestep,[1 obj.nz+1]);
Iz=squeeze(trapz(obj.rgrid,flux.*obj.rgrid)*2*pi*obj.qsim/obj.weight);
Iz(1,:)=-Iz(1,:);
gamm=obj.Metallicflux(timestep);
Im=zeros(length(gamm.p),length(timestep));
for i=1:length(gamm.p)
p=gamm.p{i};
flux=gamm.gamma{i};
for j=1:length(timestep)
Im(i,j)=pi/2*sum((p(2,1:end-1)+p(2,2:end)).*(flux(2:end,j)+flux(1:end-1,j))'...
.*(p(1,2:end)-p(1,1:end-1)));
end
end
I=-cat(1,Iz,Im*obj.qsim/obj.weight);
end
function model=potentialwellmodel(obj,timestep)
% Computes the potential well at the given timestep and return the model to be able to
% interpolate either using grid coordinates or magnetic field line coordinates
if iscell(timestep)
timestep=cell2mat(timestep);
end
% if one of the timesteps is 0 we take the external potential
id=find(timestep==0);
if(~isempty(id))
timestep(id)=[];
end
Phi=-obj.pot(:,:,timestep);
if(~isempty(id))
phiext=-obj.potxt;
Phi=cat(3,Phi(:,:,1:id-1),phiext,Phi(:,:,id:end));
end
% We get the magnetic field lines rA_\theta values
lvls=sort(unique([obj.rAthet(:,1)',obj.rAthet(:,end)', obj.rAthet(1,:),obj.rAthet(end,:)]));
% We obtaine the magnetic field lines coordinates for each
% level in r and z
contpoints=contourc(obj.zgrid,obj.rgrid,obj.rAthet,lvls(:));
[x,y,zcont]=C2xyz(contpoints);
zdiff=[diff(zcont),0];
% memory allocation for spped
potfinal=zeros(numel(cell2mat(x)),length(timestep));
pot=cell(1,size(x,2));
rathet=x;
% for each timestep we calculate the well
for i=1:length(timestep)+~isempty(id)
locPhi=Phi(:,:,i);
% We delete points outside of the domain
locPhi(obj.geomweight(:,:,1)<0)=NaN;
% For each field line we remove the lowest maximum
for j=1:size(x,2)
%for i=1:size(pot,1)
xloc=x{j};
yloc=y{j};
% xloc=xloc(yloc<rmax & yloc>rmin);
% yloc=yloc(yloc<rmax & yloc>rmin);
% x{j}=xloc;
% y{j}=yloc;
rathet{j}=zcont(j)*ones(1,length(xloc));
if(length(xloc)>=3)
pot{j}=interp2(obj.zgrid,obj.rgrid,locPhi,xloc,yloc);
locpot=pot{j};
%locpot=locpot-min(locpot);
for k=2:length(locpot)-1
left=max(locpot(1:k-1));
right=max(locpot(k+1:end));
pot{j}(k)=locpot(k)-min(left,right);
end
pot{j}(1)=0;
pot{j}(end)=0;
%pot{j}=pot{j}-max(pot{j});
else
pot{j}=zeros(size(xloc));
end
end
potfinal(:,i)=cell2mat(pot);
end
%potfinal(potfinal==0)=NaN;
model.z=cell2mat(x);
model.r=cell2mat(y);
model.pot=potfinal;
model.rathet=cell2mat(rathet);
end
function [pot] = PotentialWell(obj,timestep)
% PotentialWell Computes the potential well at the given timestep on the grid points
model=obj.potentialwellmodel(timestep);
z=model.z;
r=model.r;
modpot=model.pot;
rathet=model.rathet;
[Zmesh,Rmesh]=meshgrid(obj.zgrid,obj.rgrid);
pot=zeros(length(obj.zgrid),length(obj.rgrid),length(fieldstep));
for i=1:length(fieldstep)
pot(:,:,i)=griddata(z,r,modpot(i,:),Zmesh,Rmesh);
end
end
function Epar = Epar(obj,fieldstep)
% Computes the electric field component parallel to the magnetic field line
Epar=obj.Er(:,:,fieldstep).*(obj.Br./obj.B)' + (obj.Bz./obj.B)'.*obj.Ez(:,:,fieldstep);
end
function Eperp = Eperp(obj,fieldstep)
% Computes the electric field component perpendicular to the magnetic field line
Eperp=obj.Er(:,:,fieldstep).*(obj.Bz./obj.B)' - (obj.Br./obj.B)'.*obj.Ez(:,:,fieldstep);
end
function Ekin = Ekin(obj,varargin)
%Vpar Computes the parallel velocity for the main specie particle indices{1} at time indices{2}
if(~iscell(varargin))
indices=mat2cell(varargin);
else
indices=varargin;
end
if strcmp(indices{1},':')
p=1:obj.R.nparts;
else
p=indices{1};
end
if strcmp(indices{2},':')
t=1:length(obj.tpart);
else
t=indices{2};
end
if size(indices,1)>2
track=indices{3};
else
track=false;
end
Vr=obj.VR(p,t,track);
Vthet= obj.VTHET(p,t,track);
Vz=obj.VZ(p,t,track);
Ekin=0.5*obj.msim/obj.weight*(Vr.^2+Vthet.^2+Vz.^2);
end
function sigm=sigmio(obj,E)
if(~obj.neutcol.present || isempty(obj.neutcol.io_cross_sec))
sigm=zeros(size(E));
return
end
Ebar=obj.neutcol.scatter_fac;
Ei=obj.neutcol.Eion;
E0=obj.neutcol.E0;
nE=length(E(:));
nchi=100;
nalpha=200;
ngamma=200;
alpha=linspace(0,pi,nalpha);
gamma=linspace(0,pi,ngamma);
chi=linspace(0,0.5,nchi);
sigm2=zeros(nE,nchi);
sigm=zeros(size(E));
for i=1:nE
for j=1:nchi
sigm2(i,j)=trapz(alpha,trapz(gamma,dsigmaio(obj,E(i),Ebar,Ei,E0,chi(j),alpha,gamma),2));
end
sigm(ind2sub(size(E),i))=trapz(chi,sigm2(i,:))*obj.sigio(E(i));
end
end
function sig=sigio(obj,E)
sig=zeros(size(E));
if(~obj.neutcol.present || isempty(obj.neutcol.io_cross_sec))
return
end
for i=1:length(E(:))
if(E(i)>obj.neutcol.Eion)
sig(ind2sub(size(E),i))=obj.fit_cross_sec(E(ind2sub(size(E),i)),obj.neutcol.io_cross_sec);
end
end
end
function sigm=sigmela(obj,E)
sigm=zeros(size(E));
if(~obj.neutcol.present || isempty(obj.neutcol.ela_cross_sec))
return
end
for i=1:length(E(:))
sigm(ind2sub(size(E),i))=obj.fit_cross_sec(E(ind2sub(size(E),i)),obj.neutcol.ela_cross_sec);
end
end
function sig=sigela(obj,E)
E0=obj.neutcol.E0;
chi=E./(0.25*E0+E);
sig=(2*chi.^2)./((1-chi).*((1+chi).*log((1+chi)./(1-chi))-2*chi)).*obj.sigmela(E);
end
function [Forces, Density]=ForceBalancespline(obj,it,fdens)
%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(obj.t2d)):size(obj.t2d);
end
if nargin<3
fdens=true;
end
fluiduThet=obj.fluidUTHET(:,:,it);
Density.fluiduThet=fluiduThet;
m_e=obj.me;
q_e=obj.qe;
n=obj.N(:,:,it);
[r,~]=meshgrid(obj.rgrid,obj.zgrid);
Rinv=1./r';
Rinv(isinf(Rinv))=0;
Density.N=n;
% Calculate electric forces
Eforcer=-obj.qe*obj.Er(:,:,it);
Eforcez=-obj.qe*obj.Ez(:,:,it);
% Allocate memory
Bforcer=zeros(size(n,1),size(n,2),size(n,3));
Bforcez=zeros(size(n,1),size(n,2),size(n,3));
Bforcethet=zeros(size(n,1),size(n,2),size(n,3));
inertforcer=zeros(size(n,1),size(n,2),size(n,3));
inertforcez=zeros(size(n,1),size(n,2),size(n,3));
inertforcethet=zeros(size(n,1),size(n,2),size(n,3));
Pressforcer=zeros(size(n,1),size(n,2),size(n,3));
Pressforcethet=zeros(size(n,1),size(n,2),size(n,3));
Pressforcez=zeros(size(n,1),size(n,2),size(n,3));
durdt=zeros(size(n,1),size(n,2),size(n,3));
duthetdt=zeros(size(n,1),size(n,2),size(n,3));
duzdt=zeros(size(n,1),size(n,2),size(n,3));
if( obj.neutcol.present)
Dragforcer=zeros(size(n,1),size(n,2),size(n,3));
Dragforcethet=zeros(size(n,1),size(n,2),size(n,3));
Dragforcez=zeros(size(n,1),size(n,2),size(n,3));
end
time=obj.t2d(it);
parfor j=1:size(n,3)
Bforcer(:,:,j)=-q_e.*fluiduThet(:,:,j).*obj.Bz';
Bforcethet(:,:,j)=-q_e.*(obj.fluidUZ(:,:,it(j)).*obj.Br'-obj.fluidUR(:,:,it(j)).*obj.Bz');
Bforcez(:,:,j)=q_e.*fluiduThet(:,:,j).*obj.Br';
inertforcer(:,:,j)=-m_e.*(-fluiduThet(:,:,j).^2.*Rinv+obj.fluidUR(:,:,it(j)).*obj.fluidUR.der(:,:,it(j),[1 0])+obj.fluidUZ(:,:,it(j)).*obj.fluidUR.der(:,:,it(j),[0 1]));
inertforcethet(:,:,j)=-m_e.*(obj.fluidUR(:,:,it(j)).*fluiduThet(:,:,j).*Rinv+obj.fluidUR(:,:,it(j)).*obj.fluidUTHET.der(:,:,it(j),[1 0])+obj.fluidUZ(:,:,it(j)).*obj.fluidUTHET.der(:,:,it(j),[0 1]));
inertforcez(:,:,j)=-m_e.*(obj.fluidUR(:,:,it(j)).*obj.fluidUZ.der(:,:,it(j),[1 0])+obj.fluidUZ(:,:,it(j)).*obj.fluidUZ.der(:,:,it(j),[0 1]));
Pressforcer(:,:,j)=-( squeeze(obj.Presstens.der(1,:,:,it(j),[1 0]))...
+ squeeze(obj.Presstens(1,:,:,it(j)) - obj.Presstens(4,:,:,it(j))).*Rinv...
+ squeeze(obj.Presstens.der(3,:,:,it(j),[0 1])) )...
./n(:,:,j);
Pressforcethet(:,:,j)=-( squeeze(obj.Presstens.der(2,:,:,it(j),[1 0]))...
+ 2*squeeze(obj.Presstens(2,:,:,it(j))).*Rinv...
+ squeeze(obj.Presstens.der(5,:,:,it(j),[0 1])) )...
./n(:,:,j);
Pressforcez(:,:,j)=-( squeeze(obj.Presstens.der(3,:,:,it(j),[1 0]))...
+ squeeze(obj.Presstens(3,:,:,it(j))).*Rinv...
+ squeeze(obj.Presstens.der(6,:,:,it(j),[0 1])) )...
./n(:,:,j);
if( obj.neutcol.present)
Ek=0.5*m_e*(obj.fluidUR(:,:,it(j)).^2+obj.fluidUTHET(:,:,it(j)).^2+obj.fluidUZ(:,:,it(j)).^2);
sigm=obj.sigmela(Ek/q_e);
dragfreq=obj.neutcol.neutdens.*sigm.*sqrt(2/obj.me*Ek);
Dragforcer(:,:,j)=-m_e*dragfreq.*obj.fluidUR(:,:,it(j));
Dragforcethet(:,:,j)=-m_e*dragfreq.*obj.fluidUTHET(:,:,it(j));
Dragforcez(:,:,j)=-m_e*dragfreq.*obj.fluidUZ(:,:,it(j));
end
if (it(j)>1 && it(j)<length(obj.t2d))
cdt=(obj.t2d(it(j)+1)-obj.t2d(it(j)-1));
durdt(:,:,j)=m_e*(obj.fluidUR(:,:,it(j)+1)-obj.fluidUR(:,:,it(j)-1))/cdt;
duthetdt(:,:,j)=m_e*(obj.fluidUTHET(:,:,it(j)+1)-obj.fluidUTHET(:,:,it(j)-1))/cdt;
duzdt(:,:,j)=m_e*(obj.fluidUZ(:,:,it(j)+1)-obj.fluidUZ(:,:,it(j)-1))/cdt;
end
end
Forces.it=it;
Pressforcer(isnan(Pressforcer))=0;
Pressforcethet(isnan(Pressforcethet))=0;
Pressforcez(isnan(Pressforcez))=0;
if(~fdens)
Forces.Eforcer=Eforcer;
Forces.Eforcez=Eforcez;
Forces.Bforcer=Bforcer;
Forces.Bforcethet=Bforcethet;
Forces.Bforcez=Bforcez;
Forces.inertforcer=inertforcer;
Forces.inertforcethet=inertforcethet;
Forces.inertforcez=inertforcez;
Forces.Pressforcer=Pressforcer;
Forces.Pressforcethet=Pressforcethet;
Forces.Pressforcez=Pressforcez;
Forces.durdt=durdt;
Forces.duthetdt=duthetdt;
Forces.duzdt=duzdt;
Forces.Dragforcer=Dragforcer;
Forces.Dragforcethet=Dragforcethet;
Forces.Dragforcez=Dragforcez;
else
Forces.Eforcer=Eforcer.*n;
Forces.Eforcez=Eforcez.*n;
Forces.Bforcer=Bforcer.*n;
Forces.Bforcethet=Bforcethet.*n;
Forces.Bforcez=Bforcez.*n;
Forces.inertforcer=inertforcer.*n;
Forces.inertforcethet=inertforcethet.*n;
Forces.inertforcez=inertforcez.*n;
Forces.Pressforcer=Pressforcer.*n;
Forces.Pressforcethet=Pressforcethet.*n;
Forces.Pressforcez=Pressforcez.*n;
Forces.durdt=durdt.*n;
Forces.duthetdt=duthetdt.*n;
Forces.duzdt=duzdt.*n;
Forces.Dragforcer=Dragforcer.*n;
Forces.Dragforcethet=Dragforcethet.*n;
Forces.Dragforcez=Dragforcez.*n;
end
end
%------------------------------------------
% Functions for plotting evolving quantities
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));
f.Name= sprintf('%s Psi',obj.name);
zpos=floor(length(obj.zgrid)/2);
tinit=1;
tend=length(obj.t2d);
if iscell(deltat)
deltat=cell2mat(deltat);
end
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)
plot(obj.rgrid(In-2:end),1./obj.rgrid(In-2:end).^2*maxn*obj.rgrid(In)^2,'DisplayName','N=c*1/r^2','linewidth',lw)
plot(obj.rgrid(In-2:end),1./obj.rgrid(In-2:end).^4*maxn*obj.rgrid(In)^4,'DisplayName','N=c*1/r^4','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 f=displayrprofile(obj,t,zpos,init)
%% plot the initial and final radial profile at the middle of the simulation space
% t: time index for the
f=figure('Name', sprintf('%s Prof',obj.name));
if nargin < 3 || length(zpos)<1
zpos=floor(length(obj.zgrid)/2);
end
if nargin<4
init=false;
end
if(iscell(t))
t=cell2mat(t);
end
lw=1.5;
if init
tinit=t(1);
t=t(2:end);
end
locdens=mean(obj.N(:,zpos,t),3);
Rinv=1./obj.rgrid;
Rinv(isnan(Rinv))=0;
VTHET=mean(obj.fluidUTHET(:,zpos,t),3);
omegare=(VTHET.*Rinv);
if(init)
plot(obj.rgrid,obj.N(:,zpos,tinit),'bx-','DisplayName',sprintf('t=%1.2f[ns]',obj.t2d(tinit)*1e9),'linewidth',lw)
end
hold on
plot(obj.rgrid,locdens,'rx-','DisplayName',sprintf('t=[%1.2f-%1.2f] [ns] averaged',obj.t2d(t(1))*1e9,obj.t2d(t(end))*1e9),'linewidth',lw)
xlabel('r [m]')
ylabel('n_e [m^{-3}]')
legend('location','Northwest')
if obj.conformgeom
xlim([obj.rgrid(1) obj.rgrid(sum(obj.nnr(1:2)))])
else
xlim([obj.rgrid(1) obj.rgrid(end)])
end
grid on
ylimits=ylim();
if obj.conformgeom
plot(obj.rgrid(1)*[1 1],ylimits,'k--')
plot(obj.rgrid(end)*[1 1],ylimits,'k--')
else
plot(obj.r_a*[1 1],ylimits,'k--')
if obj.walltype==0
plot(obj.r_b*[1 1],ylimits,'k--')
elseif obj.walltype==1
rmax=obj.r_0-obj.r_r*sqrt(1-(obj.zgrid(zpos)-obj.z_0)^2/obj.z_r^2);
plot(rmax*[1 1],ylimits,'k--')
end
end
yyaxis right
plot(obj.rgrid,omegare,'DisplayName',sprintf('<\\omega_{re}> t=[%1.2f-%1.2f] [ns] averaged',obj.t2d(t(1))*1e9,obj.t2d(t(end))*1e9),'linewidth',lw)
ylabel('\omega_{re} [1/s]')
title(sprintf('Radial density profile at z=%1.2e[m]',obj.zgrid(zpos)))
obj.savegraph(f,sprintf('%srProf',obj.name),[15 10]);
end
function displayenergy(obj)
%% Plot the time evolution of the system energy and number of simulated macro particles
tmin=1;
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.etot0(tmin:tmax),'h-',...
obj.t0d(tmin:tmax),obj.eerr(tmin:tmax),'x--')
legend('E_{kin}', 'E_{pot}', 'E_{tot}','E_{ref}','E_{err}')
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.npart(tmin:tmax)./obj.npart(1)*100),'d--')
ylabel('Nparts %')
%ylim([0 110])
catch
end
obj.savegraph(f,sprintf('%s/%sEnergy',obj.folder,obj.name));
end
function f=displaycharge(obj,f,linelegend)
%% Plot the time evolution of the system charge of electrons
% f: figure handle if you want to stack several such curves
% on the same figure
% linelegend: legend for this charge evolution
tmin=1;
tmax=length(obj.ekin);
if nargin<2
f=figure('Name', sprintf('%s Charge',obj.name));
end
if nargin<3
linelegend='';
end
ax=f.CurrentAxes;
if isempty(ax)
ax=axes(f);
end
try
plot(ax,obj.t0d(tmin:tmax),abs(obj.npart(tmin:tmax)*obj.qsim),'linewidth',1.5,'displayname',linelegend)
hold on
ylabel(ax,'Total charge [C]')
xlabel(ax,'t [s]')
grid on
if(nargin>2)
legend
end
set(ax,'fontsize',12)
catch
end
if nargin < 2
obj.savegraph(f,sprintf('%s/%scharge',obj.folder,obj.name));
end
end
function displaySimParticles(obj)
%% Plot the time evolution of the number of simulated physical particles in the main specie
f=figure('Name', sprintf('%s Trapped particles',obj.name));
plot(obj.t0d,obj.npart*obj.weight,'linewidth',1.5)
xlabel('t [s]')
ylabel('N particles')
set(gca,'fontsize',12)
obj.savegraph(f,sprintf('%s/%sntrapped',obj.folder,obj.name),[10 12]);
end
function displayLarmorRad(obj,time2d)
if nargin<2
time2d=length(obj.t2d);
end
% Plot the larmor radius for created particles with low energy
rl=abs(obj.me/obj.qe*(-obj.Er(:,:,time2d).*obj.Bz'+obj.Ez(:,:,time2d).*obj.Br')./(obj.B.^3)');
figure
rl(obj.geomweight(:,:,1)<0)=0;
contourf(obj.zgrid,obj.rgrid,rl)
hold on
contour(obj.zgrid,obj.rgrid,obj.geomweight(:,:,1),[0 0],'r--')
n=obj.N(:,:,time2d);
maxN=max(n (:));
n=n/maxN*mean(rl(:));
contour(obj.zgrid,obj.rgrid,n,mean(rl(:))*linspace(0,1,6),'r:')
c=colorbar;
xlabel('z [m]')
ylabel('r [m]')
c.Label.String='r_L [m]';
end
function displayHP(obj,tstart)
% Plot the histogramm of the total energy and canonical angular momentum at time tstart and
% end time of the simulation over the full simulation space for the main specie
if(iscell(tstart))
tstart=cell2mat(tstart);
end
if(obj.R.nt>=2)
tstart=obj.R.nt;
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=H(obj,{1:partsmax,obj.R.nt,false});
%,'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=P(obj,{1:partsmax,obj.R.nt,false});
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 displayaveragetemp(obj)
% Computes and show the particles average temperature as a function of time
f=figure('Name',sprintf('%s potinn=%f part temperature',obj.name,obj.potinn*obj.phinorm));
vr2=obj.VR(:,:,false);
vr2=mean(vr2.^2,1)-mean(vr2,1).^2;
vz2=obj.VZ(:,:,false);
vz2=mean(vz2.^2,1)-mean(vz2,1).^2;
vthet2=obj.VTHET(:,:,false);
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 displayCurrentsevol(obj,timesteps)
% Computes and display the time evolution of the outgoing currents at timesteps 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','left','linewidth',1.8);
hold on
plot(obj.t2d(timesteps),currents(2,:),'Displayname','right','linewidth',1.8);
plot(obj.t2d(timesteps),currents(3,:),'Displayname','ellipse','linewidth',1.8);
plot(obj.t2d(timesteps),currents(4,:),'Displayname','outer cylinder','linewidth',1.8);
plot(obj.t2d(timesteps),currents(5,:),'Displayname','inner cylinder','linewidth',1.8);
legend('location','Northeast')
xlabel('time [s]')
ylabel('I [A]')
grid on
set(gca,'fontsize',12)
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),[16 12]);
end
function displayChargeLossevol(obj,timesteps,toptitle)
% Computes and display the time evolution of the outgoing currents at timesteps timesteps
if nargin<2
timesteps=1:length(obj.t2d);
end
currents=obj.OutCurrents(timesteps);
f=figure('Name',sprintf('%s Charges',obj.name));
t0dst=find(obj.t0d>=obj.t2d(timesteps(1)),1,'first');
t0dlst=find(obj.t0d<=obj.t2d(timesteps(end)),1,'last');
plot(obj.t2d(timesteps),currents(1,:),'Displayname','left','linewidth',1.8);
hold on
plot(obj.t2d(timesteps),currents(2,:),'Displayname','right','linewidth',1.8);
plot(obj.t2d(timesteps),currents(3,:),'Displayname','ellipse','linewidth',1.8);
plot(obj.t2d(timesteps),currents(4,:),'Displayname','outer cylinder','linewidth',1.8);
plot(obj.t2d(timesteps),currents(5,:),'Displayname','inner cylinder','linewidth',1.8);
xlabel('time [s]')
ylabel('I [A]')
grid on
set(gca,'fontsize',12)
yyaxis right
plot(obj.t0d(t0dst:t0dlst),obj.npart(t0dst:t0dlst)*obj.weight,'r:','linewidth',1.8,'Displayname','total nb e^-')
ylabel('nb e^-')
ax=gca;
ax.YAxis(2).Color='red';
legend('location','eastoutside')
if nargin <3
title(sprintf('\\phi_b-\\phi_a=%.2g kV',(obj.potout-obj.potinn)*obj.phinorm/1e3))
elseif ~isempty(toptitle)
title(toptitle)
end
obj.savegraph(f,sprintf('%s/%s_ChargeEvol',obj.folder,obj.name),[16 14]);
end
function display2Dpotentialwell(obj,timestep,rcoord)
% Display the 2D potential well at timestep timestep
% if rcoord is true, the potential is evaluated at grid points in r,z coordinates
% if false, the potential is evaluated at grid points in magnetic field line coordinates
if iscell(timestep)
timestep=cell2mat(timestep);
end
if nargin <3
rcoord=true;
end
f=figure('Name',sprintf('%s Potential well',obj.name));
model=obj.potentialwellmodel(timestep);
z=model.z;
r=model.r;
pot=model.pot;
rathet=model.rathet;
if (timestep==0)
title(sprintf('Potential well Vacuum'))
else
title(sprintf('Potential well t=%1.2f [ns]',obj.t2d(timestep)*1e9))
end
N0=obj.N(:,:,1);
id=find(timestep==0);
timestep(id)=[];
Nend=obj.N(:,:,timestep);
if(~isempty(id))
N0=zeros(obj.N.nr,obj.N.nz);
Nend=cat(3,Nend(:,:,1:id-1),N0,Nend(:,:,id:end));
end
Nend=mean(Nend,3);
geomw=obj.geomweight(:,:,1);
if rcoord
[Zmesh,Rmesh]=meshgrid(obj.zgrid,obj.rgrid);
pot=griddata(z,r,pot,Zmesh,Rmesh);
pot(obj.geomweight(:,:,1)<=0)=0;
surface(obj.zgrid(1:end),obj.rgrid(1:end),pot(1:end,1:end),'edgecolor','none','Displayname','Well')
xlabel('z [m]')
ylabel('r [m]')
xlim([obj.zgrid(1) obj.zgrid(end)])
ylim([obj.rgrid(1) obj.rgrid(end)])
hold(gca, 'on')
rdisp=obj.rgrid;
else
rdisp=sort(obj.rAthet(:,end));
[Zmesh,Rmesh]=meshgrid(obj.zgrid,rdisp);
pot=griddata(z,rathet,pot,Zmesh,Rmesh);
[Zinit,~]=meshgrid(obj.zgrid,obj.rAthet(:,1));
% if isempty(obj.maxwellsrce)
% end
N0=griddata(Zinit,obj.rAthet,N0,Zmesh,Rmesh);
Nend=griddata(Zinit,obj.rAthet,Nend,Zmesh,Rmesh);
geomw=griddata(Zinit,obj.rAthet,geomw,Zmesh,Rmesh);
pot(geomw<=0)=0;
surface(obj.zgrid(1:end),rdisp,pot(1:end,1:end),'edgecolor','none')
ylabel('rA_\theta [Tm^2]')
xlabel('z [m]')
xlim([obj.zgrid(1) obj.zgrid(end)])
ylim([min(rdisp) max(rdisp)])
hold(gca, 'on')
end
maxdensend=max(Nend(:));
contourscale=0.1;
Nend=(Nend-contourscale*maxdensend)/maxdensend;
maxdens0=max(N0(:));
contourscale=0.1;
N0=(N0-contourscale*maxdens0)/maxdens0;
contour(obj.zgrid,rdisp,Nend,linspace(0,1-contourscale,5),'k--','linewidth',1.5,'Displayname','Cloud Boundaries');
contour(obj.zgrid,rdisp,N0,[0 0],'k-.','linewidth',1.5,'Displayname','Source boundaries');
contour(obj.zgrid,rdisp,geomw,[0 0],'w-.','linewidth',1.5,'Displayname','Vessel Boundaries');
c=colorbar;
colormap('jet');
c.Label.String= 'depth [eV]';
if rcoord
obj.savegraph(f,sprintf('%s/%s_wellr',obj.folder,obj.name));
else
obj.savegraph(f,sprintf('%s/%s_wellpsi',obj.folder,obj.name));
end
end
function display1Dpotentialwell(obj,timestep,rpos)
% Display the potential well along the magentic field line
% passing by rgrid(rpos) at the center of the simulation space
if iscell(timestep)
timestep=cell2mat(timestep);
end
f=figure('Name',sprintf('%s 1D Potential well',obj.name));
model=obj.potentialwellmodel(timestep);
z=model.z;
r=model.r;
pot=model.pot;
rathet=model.rathet;
if (mod(rpos, 1) ~= 0)
[~,rpos]=min(abs(M.rgrid-rpos));
end
crpos=obj.rgrid(rpos);
id=find(timestep==0);
timestep(id)=[];
N=obj.N(:,:,timestep);
if(~isempty(timestep==0))
N0=zeros(obj.N.nr+1,obj.N.nz+1);
N=cat(3,N(:,:,1:id-1),N0,N(:,:,id:end));
end
N=mean(N,3);
linepot=zeros(length(obj.zgrid),length(timestep));
rathetpos=obj.rAthet(rpos,ceil(length(obj.zgrid)/2));
F=scatteredInterpolant(z',rathet',pot(:,1));
for i=1:length(timestep)
F=scatteredInterpolant(z',rathet',pot(:,i));
linepot(:,i)=F(obj.zgrid,rathetpos*ones(size(obj.zgrid)));
%linepot(:,i)=griddata(z,rathet,pot(:,i),obj.zgrid,rathetpos);
end
linepot=mean(linepot,2);
[Zinit,~]=meshgrid(obj.zgrid,obj.rAthet(:,1));
N=griddata(Zinit,obj.rAthet,N,obj.zgrid,rathetpos);
plot(obj.zgrid,linepot)
ylabel('Potentiel [eV]')
xlabel('z [m]')
xlim([obj.zgrid(1) obj.zgrid(end)])
hold(gca, 'on')
yyaxis right
plot(obj.zgrid,N)
ylabel('n [m^{-3}]')
if length(timestep)==1
title(sprintf('Potential well t=%1.2f [ns] r=%1.2f [mm]',obj.t2d(timestep)*1e9,1e3*crpos))
else
title(sprintf('Potential well t=[%1.2f-%1.2f] [ns] r=%1.2f [mm]',obj.t2d(timestep(1))*1e9,obj.t2d(timestep(end))*1e9,1e3*crpos))
end
obj.savegraph(f,sprintf('%s/%s_well1Dr_%d',obj.folder,obj.name,rpos));
end
function displayVdistribRThetZ(obj,timestep, rpos, zpos)
if(obj.R.nt>=2)
if nargin<2
timesteppart=length(obj.tpart);
else
timesteppart=timestep;
end
if nargin<3 || isempty(rpos)
rspan=[obj.rgrid(1) obj.rgrid(end)];
else
r=[obj.rgrid(1);(obj.rgrid(1:end-1)+obj.rgrid(2:end))*0.5;obj.rgrid(end)];
rspan=[r(rpos) r(rpos+1)];
end
if nargin<4 || isempty(zpos)
zspan=[obj.zgrid(1) obj.zgrid(end)];
else
z=[obj.zgrid(1);(obj.zgrid(1:end-1)+obj.zgrid(2:end))*0.5;obj.zgrid(end)];
zspan=[z(zpos) z(zpos+1)];
end
nbp=min(obj.nbparts(1),obj.R.nparts);
R=obj.R(1:nbp,1,false);
Z=obj.Z(1:nbp,1,false);
Vr=obj.VR(1:nbp,1,false);
Vz=obj.VZ(1:nbp,1,false);
Vthet=obj.VTHET(1:nbp,1,false);
ids=R>=rspan(1) & R<=rspan(2) & Z>=zspan(1) & Z<=zspan(2);
Vr=Vr(ids);
Vz=Vz(ids);
Vthet=Vthet(ids);
vTr=std(Vr,1);
vTz=std(Vz,1);
vTthet=std(Vthet,1);
nbp=min(obj.nbparts(timesteppart),obj.R.nparts);
Rend=obj.R(1:nbp,timesteppart,false);
Zend=obj.Z(1:nbp,timesteppart,false);
Vrend=obj.VR(1:nbp,timesteppart,false);
Vzend=obj.VZ(1:nbp,timesteppart,false);
Vthetend=obj.VTHET(1:nbp,timesteppart,false);
ids=Rend>=rspan(1) & Rend<=rspan(2) & Zend>=zspan(1) & Zend<=zspan(2);
Vrend=Vrend(ids);
Vzend=Vzend(ids);
Vthetend=Vthetend(ids);
vTrend=std(Vrend,1);
vTzend=std(Vzend,1);
vTthetend=std(Vthetend,1);
binwidth=abs(max(Vrend)-min(Vrend))/sqrt(length(Vrend));
f=figure('Name',sprintf("%s vrz distrib",obj.file));
subplot(1,3,1);
obj.dispV(Vr,Vrend,'V_r',[1,timesteppart])
subplot(1,3,2);
obj.dispV(Vthet,Vthetend,'V\theta',[1,timesteppart])
subplot(1,3,3);
obj.dispV(Vz,Vzend,'Vz',[1,timesteppart])
sgtitle(sprintf('dt=%1.2e[ns] r=[%2.1f, %2.1f] [mm] z=[%2.1f, %2.1f] [mm]',obj.dt*1e9, rspan*1e3, zspan*1e3))
obj.savegraph(f,sprintf('%s/%sParts_V_RZ',obj.folder,obj.name));
end
end
function displayEkin(obj,timestep, rpos, zpos)
if(obj.R.nt>=2)
if nargin<2
timesteppart=[1 length(obj.tpart)];
else
if length(timestep)<2
timesteppart=[1 timestep];
else
timesteppart=[timestep(1) timestep(end)];
end
end
if nargin<3 || isempty(rpos)
rspan=[obj.rgrid(1) obj.rgrid(end)];
else
r=[obj.rgrid(1);(obj.rgrid(1:end-1)+obj.rgrid(2:end))*0.5;obj.rgrid(end)];
rspan=[r(rpos) r(rpos+1)];
end
if nargin<4 || isempty(zpos)
zspan=[obj.zgrid(1) obj.zgrid(end)];
else
z=[obj.zgrid(1);(obj.zgrid(1:end-1)+obj.zgrid(2:end))*0.5;obj.zgrid(end)];
zspan=[z(zpos) z(zpos+1)];
end
nbp=min(obj.nbparts(timesteppart(1)),obj.R.nparts);
R=obj.R(1:nbp,timesteppart(1),false);
Z=obj.Z(1:nbp,timesteppart(1),false);
Ekin=obj.Ekin(1:nbp,timesteppart(1),false);
ids=R>=rspan(1) & R<=rspan(2) & Z>=zspan(1) & Z<=zspan(2);
Ekin=Ekin(ids)/obj.qe;
nbp=min(obj.nbparts(timesteppart(2)),obj.R.nparts);
Rend=obj.R(1:nbp,timesteppart(2),false);
Zend=obj.Z(1:nbp,timesteppart(2),false);
Ekinend=obj.Ekin(1:nbp,timesteppart(2),false);
ids=Rend>=rspan(1) & Rend<=rspan(2) & Zend>=zspan(1) & Zend<=zspan(2);
Ekinend=Ekinend(ids)/obj.qe;
f=figure('Name',sprintf("%s E_k distrib",obj.file));
obj.dispV(Ekin,Ekinend,'E_k',timesteppart,'maxwell')
sgtitle(sprintf('dt=%1.2e[ns] r=[%2.1f, %2.1f] [mm] z=[%2.1f, %2.1f] [mm]',obj.dt*1e9, rspan*1e3, zspan*1e3))
obj.savegraph(f,sprintf('%s/%sParts_V_RZ',obj.folder,obj.name));
end
end
function displayVdistribParPer(obj,timestep, rpos, zpos)
if(obj.R.nt>=2)
if nargin<2
timesteppart=length(obj.tpart);
else
timesteppart=timestep;
end
if nargin<3 || isempty(rpos)
rspan=[obj.rgrid(1) obj.rgrid(end)];
else
r=[obj.rgrid(1);(obj.rgrid(1:end-1)+obj.rgrid(2:end))*0.5;obj.rgrid(end)];
rspan=[r(rpos) r(rpos+1)];
end
if nargin<4 || isempty(zpos)
zspan=[obj.zgrid(1) obj.zgrid(end)];
else
z=[obj.zgrid(1);(obj.zgrid(1:end-1)+obj.zgrid(2:end))*0.5;obj.zgrid(end)];
zspan=[z(zpos) z(zpos+1)];
end
nbp=min(obj.nbparts(1),obj.R.nparts);
R=obj.R(1:nbp,1,false);
Z=obj.Z(1:nbp,1,false);
ids=R>=rspan(1) & R<=rspan(2) & Z>=zspan(1) & Z<=zspan(2);
Vperp=obj.Vperp(1:nbp,1,false);
Vpar=obj.Vpar(1:nbp,1,false);
Vperp=Vperp(ids);
Vpar=Vpar(ids);
nbp=min(obj.nbparts(timesteppart),obj.R.nparts);
R=obj.R(1:nbp,timesteppart,false);
Z=obj.Z(1:nbp,timesteppart,false);
ids=R>=rspan(1) & R<=rspan(2) & Z>=zspan(1) & Z<=zspan(2);
Vperpend=obj.Vperp(1:nbp,timesteppart,false);
Vparend=obj.Vpar(1:nbp,timesteppart,false);
Vperpend=Vperpend(ids);
Vparend=Vparend(ids);
%binwidth=abs(max(Vparend)-min(Vparend))/500;
f=figure('Name',sprintf("%s v parper distrib",obj.file));
subplot(1,2,1)
obj.dispV(Vperp,Vperpend,'v_\perp',[1,timesteppart], 'None')
subplot(1,2,2)
obj.dispV(Vpar,Vparend,'v_{par}',[1,timesteppart],'None')
sgtitle(sprintf('dt=%1.2e[ns] r=[%2.1f, %2.1f] [mm] z=[%2.1f, %2.1f] [mm]',obj.dt*1e9, rspan*1e3, zspan*1e3))
obj.savegraph(f,sprintf('%s/%sParts_V_parper',obj.folder,obj.name));
end
end
function [p, maxnb, c]=displayPhaseSpace(obj,type,partsstep, Rindex, Zindex,legendtext, figtitle, f, maxnb, c, gcs)
if nargin<8
f=figure;
f=gca;
end
if nargin<7
figtitle=sprintf('r=%1.2f [mm] z=%1.2f [mm] \\Delta\\phi=%1.1f[kV] R=%1.1f',obj.rgrid(Rindex)*1e3,obj.zgrid(Zindex)*1e3,(obj.potout-obj.potinn)*obj.phinorm/1e3,obj.Rcurv);
end
if nargin <6
legendtext=sprintf('t=%1.3g [s]',obj.tpart(partsstep));
end
fieldstep=find(obj.tpart(partsstep(end))==obj.t2d,1);
if nargin>=10
ctemp=c;
n=zeros(length(c{1}),length(c{2}));
else
nbins=15;
n=zeros(nbins);
end
if nargin <11
gcs=true;
end
for i=1:length(partsstep)
odstep=find(obj.tpart(partsstep(i))==obj.t0d);
nbp=min(obj.R.nparts,obj.nbparts(partsstep(i)));
Rp=obj.R(1:nbp,partsstep(i),false);
Zp=obj.Z(1:nbp,partsstep(i),false);
deltar=obj.dr(2)/2;
deltarm=obj.rgrid(Rindex)-sqrt(obj.rgrid(Rindex)^2-deltar^2-2*obj.rgrid(Rindex)*deltar);
deltaz=obj.dz/2;
Indices=Rp>=obj.rgrid(Rindex)-deltarm & Rp<obj.rgrid(Rindex)+deltar & Zp>=obj.zgrid(Zindex)-deltaz & Zp<obj.zgrid(Zindex)+deltaz;
if strcmp(type,'parper')
Vperp=obj.Vperp(Indices,partsstep(i),false,gcs);
Vpar=obj.Vpar(Indices,partsstep(i),false,gcs);
if exist('ctemp','var')
[ntemp,ctemp] = hist3([Vpar, Vperp],'Ctrs',ctemp );
else
[ntemp,ctemp]=hist3([Vpar, Vperp],nbins*[1 1]);
end
else
Vr=obj.VR(Indices,partsstep(i),false);
Vthet=obj.VTHET(Indices,partsstep(i),false);
if exist('ctemp','var')
[ntemp,ctemp] = hist3([Vr, Vthet],'Ctrs',ctemp );
else
[ntemp,ctemp]=hist3([Vr, Vthet],nbins*[1 1]);
end
end
n=(n+ntemp)/2;
end
c=ctemp;
if nargin<8 || maxnb < 0
maxnb=max(max(n));
end
sumn=sum(n(:))
contourf(f,c{1},c{2},n'/maxnb,'Displayname',legendtext);
colorbar(f)
hold(f,'on')
%caxis(f,[0 1])
axis(f, 'equal')
if strcmp(type,'parper')
xlabel(f,'v_{par} [m/s]')
if gcs
ylabel(f,'v_\perp^* [m/s]')
else
ylabel(f,'v_\perp [m/s]')
end
else
xlabel(f,'v_r [m/s]')
ylabel(f,'v_\theta [m/s]')
end
xlimits=xlim(f);
ylimits=ylim(f);
xtickformat(f,'%.2g')
ytickformat(f,'%.2g')
if strcmp(type,'parper')
b=obj.rgrid(end);
a=obj.rgrid(1);
% drift velocity in vacuum vessel with bias
vd1=((obj.potout-obj.potinn)*obj.phinorm/obj.rgrid(Rindex))/obj.Bz(Zindex,Rindex)/log(b/a);
% initial particule velocity ( thermal velocity)
vinit=sqrt(obj.kb*22000/obj.me);
[Zmesh,~]=meshgrid(obj.zgrid,obj.rAthet(:,Zindex));
Zeval=obj.zgrid(1:end);
zleftlim=1;
zrightlim=length(Zeval);
Psieval=obj.rAthet(Rindex,Zindex)*ones(length(Zeval),1);
% Electrostatic potential on magnetic field line
% coordinates
phis=griddata(Zmesh,obj.rAthet,obj.pot(:,:,fieldstep),Zeval,Psieval,'natural');
% Magnetic field mirror ratio at each grid position
% compared to local position
R=griddata(Zmesh,obj.rAthet,obj.B',Zeval,Psieval,'natural')/obj.B(Zindex,Rindex);
if ~gcs
[~,tfield]=min(abs(obj.t2d-obj.tpart(partsstep(1))));
timeEr=obj.Er(:,:,tfield);
timeEz=obj.Ez(:,:,tfield);
posindE=sub2ind(size(timeEr),Rindex,Zindex);
% ExB drift velocity at considered time step
Vdrift=(timeEz(posindE).*obj.Br(Zindex,Rindex)-timeEr(posindE).*obj.Bz(Zindex,Rindex))./obj.B(Zindex,Rindex).^2;
else
Vdrift=0;
end
Rcurvl=R(zleftlim);
Rcurvr=R(zrightlim);
deltaphicompl=0.5*obj.me/obj.qe*((vd1+vinit)^2*(Rcurvl-1));
deltaphil=-phis(Zindex)+phis(zleftlim);
deltaphicompr=0.5*obj.me/obj.qe*((vd1+vinit)^2*(Rcurvr-1));
deltaphir=-phis(Zindex)+phis(zrightlim);
vpar=linspace(1,3*max(abs(Vpar)),1000);
vperstarr=sqrt((2*obj.qe/obj.me*deltaphir+vpar.^2)/(Rcurvr-1));
vperr=sqrt(vperstarr.^2+Vdrift^2-2*Vdrift*abs(vperstarr));
vparr=vpar(real(vperr)~=0& imag(vperr)==0);
vperr=vperr(real(vperr)~=0& imag(vperr)==0);
vperstarl=sqrt((2*obj.qe/obj.me*deltaphil+vpar.^2)/(Rcurvl-1));
vperl=sqrt(vperstarl.^2+Vdrift^2-2*Vdrift*abs(vperstarl));
vparl=vpar(real(vperl)~=0 & imag(vperl)==0);
vperl=vperl(real(vperl)~=0 & imag(vperl)==0);
p=plot(f,vparr,vperr,'r-','displayname','Simulation','linewidth',2);
plot(f,-vparl,vperl,'r--','linewidth',2)
% Prediction using model
vpar=linspace(1,3*max(abs(Vpar)),1000);
vperstarr=sqrt((2*obj.qe/obj.me*deltaphicompr+vpar.^2)/(Rcurvr-1));
vperr=sqrt(vperstarr.^2+Vdrift^2-2*Vdrift*abs(vperstarr));
vparr=vpar(real(vperr)~=0);
vperr=vperr(real(vperr)~=0);
vperstarl=sqrt((2*obj.qe/obj.me*deltaphicompl+vpar.^2)/(Rcurvl-1));
vperl=sqrt(vperstarl.^2+Vdrift^2-2*Vdrift*abs(vperstarl));
vparl=vpar(real(vperl)~=0);
vperl=vperl(real(vperl)~=0);
p=plot(f,vparr,vperr,'r-.','displayname','Prediction','linewidth',2);
plot(f,-vparl,vperl,'r:','linewidth',2)
else
p=gobjects(0);
end
xlim(f,xlimits)
ylim(f,ylimits)
legend(p(isgraphics(p)))
legend(p(isgraphics(p)),'location','north')
title(f,figtitle)
end
function displaydensity(obj,fieldstart,fieldend)
f=figure('Name', sprintf('%sfluid_dens',obj.name));
ax1=gca;
geomw=obj.geomweight(:,:,1);
dens=mean(obj.N(:,:,fieldstart:fieldend),3);
dens(geomw<=0)=0;
maxdens=max(dens(:));
h=surface(ax1,obj.zgrid,obj.rgrid,dens);
set(h,'edgecolor','none');
hold on
contour(ax1,obj.zgrid,obj.rgrid,obj.geomweight(:,:,1),[0 0],'r-.','linewidth',1.5,'Displayname','Boundaries');
if(obj.conformgeom)
ylim(ax1,[obj.rgrid(1) obj.rgrid(sum(obj.nnr(1:2)))])
else
ylim(ax1,[obj.rgrid(1) obj.rgrid(end)])
end
xlim(ax1,[obj.zgrid(1) obj.zgrid(end)])
xlabel(ax1,'z [m]')
ylabel(ax1,'r [m]')
title(ax1,sprintf('Density t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',obj.t2d(fieldstart),obj.t2d(fieldend),double(maxdens)))
c = colorbar(ax1);
c.Label.String= 'n[m^{-3}]';
view(ax1,2)
grid on;
hold on;
f.PaperOrientation='landscape';
f.PaperUnits='centimeters';
papsize=[16 14];
f.PaperSize=papsize;
print(f,sprintf('%sfluid_dens',obj.name),'-dpdf','-fillpage')
savefig(f,sprintf('%sfluid_dens',obj.name))
end
function displaySurfFluxes(obj,timesteps)
mflux= obj.Metallicflux(timesteps);
lflux= -squeeze(obj.Axialflux(timesteps,1))';
rflux= squeeze(obj.Axialflux(timesteps,length(obj.zgrid)))';
time=obj.t2d(timesteps);
%%
f=figure('name','fluxevol');
ax1=subplot(2,2,1);
contourf(mflux.p{1}(1,:)*100,time*1e9,mflux.gamma{1}'*obj.qe/(100^2),'linestyle','none')
title('ellipse')
xlabel('z [cm]')
c=colorbar;
c.Label.String= 'j\cdotn [A/cm^2]';
ax2=subplot(2,2,2);
contourf(mflux.p{2}(1,:)*100,time*1e9,mflux.gamma{2}'*obj.qe/(100^2),'linestyle','none')
title('top')
xlabel('z [cm]')
c=colorbar;
c.Label.String= 'j\cdotn [A/cm^2]';
ax3=subplot(2,2,3);
contourf(obj.rgrid*100,time*1e9,lflux*obj.qe/(100^2),'linestyle','none')
title('left')
xlabel('r [cm]')
c=colorbar;
c.Label.String= 'j\cdotn [A/cm^2]';
ax4=subplot(2,2,4);
contourf(obj.rgrid*100,time*1e9,rflux*obj.qe/(100^2),'linestyle','none')
title('right')
xlabel('r [cm]')
c=colorbar;
c.Label.String= 'j\cdotn [A/cm^2]';
ax=[ax1, ax2, ax3,ax4];
ylabel(ax,'t [ns]')
obj.savegraph(f,sprintf('%s/%s_surfFluxEvol',obj.folder,obj.name),[16 14]);
end
dispespicPressure(obj,logdensity,showgrid,fixed,temperature)
dispespicFields(obj,logdensity,showgrid,fixed,parper)
%------------------------------------------
% Helper functions needed for other functions
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 dispV(obj,V,Vend,label,t, dist)
if nargin<6
dist='gaussian';
end
vmean=mean(V(~isnan(V)));
vtherm=std(V(~isnan(V)),1);
vmeanend=mean(Vend(~isnan(Vend)));
vthermend=std(Vend(~isnan(Vend)),1);
if(length(V)>1)
[Counts,edges]=histcounts(V,'binmethod','scott');
binwidth=mean(diff(edges));
plot([edges(1) 0.5*(edges(2:end)+edges(1:end-1)) edges(end)],[0 Counts 0],'DisplayName',sprintf("t=%2.3d [ns]",obj.tpart(t(1))*1e9));
hold on
end
hold on
[Counts,edges]=histcounts(Vend,'binmethod','scott');
plot([edges(1) 0.5*(edges(2:end)+edges(1:end-1)) edges(end)],[0 Counts 0],'DisplayName',sprintf("t=%2.3d [ns]",obj.tpart(t(2))*1e9));
vfit=linspace(edges(1),edges(end),300);
if strcmp(dist,'maxwell')
a=vmeanend/2*sqrt(pi/2);
dist=sqrt(2/pi)*vfit.^2.*exp(-(vfit).^2/2/a^2)/a^3;
dist=dist/max(dist);
plot(vfit,max(Counts)*dist,'displayname',sprintf('Maxwellian fit mu=%2.2g sigma=%2.2g',vmeanend,vthermend))
elseif strcmp(dist,'gaussian')
dist=exp(-(vfit-vmeanend).^2/2/vthermend^2);
plot(vfit,max(Counts)*dist,'displayname',sprintf('gaussian fit mu=%2.2g sigma=%2.2g',vmeanend,vthermend))
end
ylabel('counts')
xlabel(label)
grid on
legend('location','northwest','orientation','vertical')
end
function cross_sec=fit_cross_sec(obj,energy,crosssec_table)
cross_sec=0;
if (energy<=0 || isnan(energy) || isinf(energy))
return
end
id=find(energy>crosssec_table(:,1),1,'last');
if(isempty(id))
id=1;
end
id=id+1;
id=min(size(crosssec_table,1),id);
id=max(2,id);
cross_sec=crosssec_table(id-1,2)*(energy/crosssec_table(id-1,1))^crosssec_table(id-1,3);
end
function fighandle=savegraph(obj, fighandle, name, papsize)
%% Saves the given figure as a pdf and a .fig
fighandle.PaperUnits='centimeters';
if (nargin < 4)
papsize=[14 16];
end
fighandle.PaperSize=papsize;
print(fighandle,name,'-dpdf','-fillpage')
savefig(fighandle,name)
end
function sig=dsigmaio(obj,Ekin, Ebar, Ei, E0, chi, alpha, gamma)
alpha=reshape(alpha,[],1);
gamma=reshape(gamma,1,[]);
sigalpha=(E0^2+8*chi*(Ekin-Ei)*E0)./(E0+4*chi*(Ekin-Ei)-4*chi*(Ekin-Ei)*cos(alpha)).^2/4/pi;
siggamma=(E0^2+8*(1-chi)*(Ekin-Ei)*E0)./(E0+4*(1-chi)*(Ekin-Ei)-4*(1-chi)*(Ekin-Ei)*cos(gamma)).^2/4/pi;
sigalpha=reshape(sigalpha,[],1);
siggamma=reshape(siggamma,1,[]);
sigchi=1/(Ebar*atan((Ekin-Ei)/(2*Ebar)).*(1+(chi*(Ekin-Ei)/Ebar).^2.1));
dp=1-sqrt(chi*(1-Ei/Ekin))*cos(alpha)- sqrt((1-chi)*(1-Ei/Ekin))*cos(gamma);
sig=sigalpha*sigchi*siggamma.*dp;
end
end
end
c4science · Help