diff --git a/matlab/espic2dhdf5.m b/matlab/espic2dhdf5.m
index a35e76c..145917c 100644
--- a/matlab/espic2dhdf5.m
+++ b/matlab/espic2dhdf5.m
@@ -1,1399 +1,1421 @@
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 particles
%% 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
phi % Electric potential in spline form
Er % Radial electric field
Ez % Axial electric field
Presstens % Pressure tensor
%% 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
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
% 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;
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.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
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);
else
obj.Presstens=splinepressure(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, obj.vnorm^2*obj.msim, 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
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
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;
for i=1:obj.nbcelldiag
nbparts=h5read(obj.fullpath,sprintf('%s/Nparts',sprintf('/data/celldiag/%02d',i)));
if (sum(nbparts)>0)
obj.celldiag(i)=h5parts(obj.fullpath,sprintf('/data/celldiag/%02d',i),obj);
end
end
end
end
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 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
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 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 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.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.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 SimParticles(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)
xlabel('t [s]')
ylabel('N particles')
obj.savegraph(f,sprintf('%s/%sntrapped',obj.folder,obj.name));
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 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 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 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
I=zeros(2,length(timestep));
flux=obj.Axialflux(timestep,[1 obj.nz+1]);
I=squeeze(trapz(obj.rgrid,flux.*obj.rgrid)*2*pi*obj.qsim/obj.weight);
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',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]')
grid on
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
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
Phi=-obj.pot(:,:,timestep);
lvls=sort(unique([obj.rAthet(:,end)' obj.rAthet(1,:)]));
contpoints=contourc(obj.zgrid,obj.rgrid,obj.rAthet,lvls(:));
[x,y,zcont]=C2xyz(contpoints);
k=1;
zdiff=[diff(zcont),0];
potfinal=zeros(numel(cell2mat(x)),length(timestep));
pot=cell(1,size(x,2));
rmin=obj.rgrid(1);
rmax=obj.rgrid(end);
rathet=x;
for i=1:length(timestep)
locPhi=Phi(:,:,i);
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)>=1)
pot{j}=interp2(obj.zgrid,obj.rgrid,locPhi,xloc,yloc,'makima');
pot{j}=pot{j}-max(pot{j});
end
k=k+(zdiff(j)~=0);
end
potfinal(:,i)=cell2mat(pot);
end
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 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;
title(sprintf('Potential well t=%1.2f [ns]',obj.t2d(timestep)*1e9))
N0=obj.N(:,:,1);
Nend=obj.N(:,:,timestep);
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);
N=mean(obj.N(:,:,timestep),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 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 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 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(end),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 dispV(obj,V,Vend,label,t, perp)
+ function dispV(obj,V,Vend,label,t, dist)
if nargin<6
- perp=false
+ 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,'binwidth',binwidth);
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 perp
+ 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))
- else
+ 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 displayVdistribParPer(obj,timestep)
+ 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(end),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], true)
+ obj.dispV(Vperp,Vperpend,'v_\perp',[1,timesteppart], 'None')
subplot(1,2,2)
- obj.dispV(Vpar,Vparend,'v_{par}',[1,timesteppart])
+ obj.dispV(Vpar,Vparend,'v_{par}',[1,timesteppart],'None')
- sgtitle(sprintf('dt=%1.2e[ns]',obj.dt*1e9))
+ 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
end
end