fennecslitehdf5.m
cdf25dab51fdmaster

fennecslitehdf5.m
View Options

			classdef fennecslitehdf5
			% General class used to treat hdf5 result files of fennecs code
			% A result file is loaded with a call to M=fennecslitehdf5(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
			restartsteps
			restarttimes
			fileversion

			%% 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
			nthet % number of intervals in the azimuthal 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 % electric potential at the coaxial insert
			potout % 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
			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
			thetgrid % azimuthal 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
			celltype % type of cell -1 outside 1 inside 0 border
			linked_s % location of linked spline
			bsplinetype

			%% 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

			%% 2D time data evaluated on grid points
			N % main species density used to compute currents
			fluidUR % main species radial fluid velocity
			fluidUZ % main species axial fluid velocity
			fluidUTHET % main species azimuthal fluid velocity
			pot % electric potential evaluated at grid points
			potxt % external electric potential evaluated at grid points
			phi % electric potential in spline form
			phireal % Fourier modes of the potential (real part)
			phiimag % Fourier modes of the potential (imaginary part)
			rhs % RHS of the Poisson equation
			rhsreal % Fourier modes of the RHS (real part)
			rhsimag % Fourier modes of the RHS (imaginary part)
			Er % radial electric field
			Ethet % azimuthal electric field
			Ez % axial electric field
			Erxt % external radial electric field
			Ezxt % external axial electric field

			%% 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
			npart % time evolution of the number of simulated particles
			losses % number of particles lost on each boundary
			species % array containing the species
			species_moments % array containing moments for species (parts)

			%% 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
			dirichletweight
			gtilde
			fsolution
			spl_bound

			%% Maxwell source parameters
			maxwellsrce

			%% Collision with neutral parameters
			neutcol

			%% Non ideal power supply
			psupply

			end

			methods
			function file=file(obj)
			% returns the h5 file name
			file=obj.filename;
			end

			function obj = fennecslitehdf5(filename,readparts)
			% Read the new result file filename and read the parts data if readparts==true

			% Add the helper_classes folder to the path
			matlabfuncpath = dir([mfilename('fullpath'),'.m']);
			addpath(sprintf('%s/../helper_classes',matlabfuncpath.folder));
			addpath(sprintf('%s/../helper_functions',matlabfuncpath.folder));
			addpath(sprintf('%s/../export_fig',matlabfuncpath.folder));
			rehash path
			filedata=dir(filename);
			if (isempty(filedata))
			fullpath=which(filename);
			filedata=dir(fullpath);
			if (isempty(filedata))
			error("File: ""%s"" doesn't exist",filename)
			end
			end
			obj.folder=filedata.folder;
			obj.filename=filedata.name;
			[~, 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

			%% Read the run parameters
			obj.dt = h5readatt(obj.fullpath,'/data/input.00/','dt');
			obj.fileversion = h5readatt(obj.fullpath,'/data','file_version');
			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.qsim = h5readatt(obj.fullpath,'/data/input.00/','qsim');
			obj.msim = h5readatt(obj.fullpath,'/data/input.00/','msim');
			obj.weight=h5readatt(obj.fullpath,'/data/part/','weight');

			obj.nz = h5readatt(obj.fullpath,'/data/input.00/','nz');
			obj.nthet = h5readatt(obj.fullpath,'/data/input.00/','nthet');
			obj.nnr = h5read(obj.fullpath,'/data/input.00/nnr');
			obj.lz = h5read(obj.fullpath,'/data/input.00/lz');
			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.L_r=h5readatt(obj.fullpath,'/data/input.00/geometry','L_r');
			obj.L_z=h5readatt(obj.fullpath,'/data/input.00/geometry','L_z');
			obj.conformgeom=false;

			nbrst=h5readatt(obj.fullpath,'/files','jobnum');
			obj.restartsteps(1)=0;
			obj.restarttimes(1)=0;
			grp=sprintf('/data/input.%02i/',0);
			obj.dt(1)=h5readatt(obj.fullpath,grp,'dt');
			for i=1:nbrst
			grp=sprintf('/data/input.%02i/',i);
			obj.restartsteps(i+1)= h5readatt(obj.fullpath,grp,'startstep');
			obj.restarttimes(i+1)= obj.restarttimes(i) + (obj.restartsteps(i+1)-obj.restartsteps(i))*obj.dt(i);
			obj.dt(i+1)=h5readatt(obj.fullpath,grp,'dt');
			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');
			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.femorder = h5read(obj.fullpath,'/data/input.00/femorder');

			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');

			obj.nbspecies=h5readatt(obj.fullpath,'/data/input.00/','nbspecies');
			obj.normalized=strcmp(h5readatt(obj.fullpath,'/data/input.00/','rawparts'),'y');

			obj.omepe=sqrt(abs(obj.n0)obj.qe^2/(obj.meobj.eps_0));
			obj.omece=obj.qe*obj.B0/obj.me;
			try
			obj.losses=h5read(obj.fullpath,'/data/var0d/losses');
			catch
			end

			obj.tnorm=min(abs(1/obj.omepe),abs(1/obj.omece));
			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;
			% Conversion of the bias to V
			obj.potinn=obj.potinn*obj.phinorm;
			obj.potout=obj.potout*obj.phinorm;

			% 0D data
			obj.t0d=h5read(obj.fullpath,'/data/var0d/time');
			obj.epot = h5read(obj.fullpath,'/data/var0d/epot');
			obj.ekin = h5read(obj.fullpath,'/data/var0d/ekin');
			obj.etot = h5read(obj.fullpath,'/data/var0d/etot');
			obj.etot0 = h5read(obj.fullpath,'/data/var0d/etot0');
			obj.eerr = obj.etot-obj.etot0;

			% Grid data
			obj.rgrid= h5read(obj.fullpath, '/data/var1d/rgrid')*obj.rnorm;
			obj.thetgrid= h5read(obj.fullpath, '/data/var1d/thetgrid');
			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;
			Bz = h5read(obj.fullpath,'/data/fields/Bz')*obj.bnorm;
			try
			Atheta = h5read(obj.fullpath,'/data/fields/Athet');
			obj.Athet= reshape(Atheta,length(obj.zgrid),length(obj.rgrid));
			[rmeshgrid,~]=meshgrid(obj.rgrid,obj.zgrid);
			obj.rAthet=(rmeshgrid.*obj.Athet)';
			catch
			end
			obj.Br= reshape(Br,length(obj.rgrid),length(obj.zgrid))';
			obj.Bz= reshape(Bz,length(obj.rgrid),length(obj.zgrid))';
			obj.Athet= reshape(Atheta,length(obj.rgrid),length(obj.zgrid))';
			[rmeshgrid,~]=meshgrid(obj.rgrid,obj.zgrid);
			obj.rAthet=(rmeshgrid.*obj.Athet)';
			obj.B=sqrt(obj.Bz.^2+obj.Br.^2);
			obj.costhet=(obj.Br./obj.B)';
			obj.sinthet=(obj.Bz./obj.B)';
			clear Br Bz
			for i=0:nbrst
			grp=sprintf('/data/input.%02i/',i);
			obj.Erxt(:,:,i+1)=reshape(h5read(obj.fullpath,[grp,'Erxt']),length(obj.rgrid),length(obj.zgrid))*obj.enorm;
			obj.Ezxt(:,:,i+1)=reshape(h5read(obj.fullpath,[grp,'Ezxt']),length(obj.rgrid),length(obj.zgrid))*obj.enorm;

			obj.potxt(:,:,i+1)=reshape(h5read(obj.fullpath,[grp,'potxt']),length(obj.rgrid),length(obj.zgrid))*obj.phinorm;
			end

			obj.t2d = h5read(obj.fullpath,'/data/fields/time');
			obj.pot=gridquantity3D(obj.fullpath,'/data/fields/pot',sum(obj.nnr)+1, obj.nthet, obj.nz+1,1);
			obj.Er=gridquantity3D(obj.fullpath,'/data/fields/Er',sum(obj.nnr)+1, obj.nthet, obj.nz+1,1);
			obj.Ethet=gridquantity3D(obj.fullpath,'/data/fields/Ethet',sum(obj.nnr)+1, obj.nthet, obj.nz+1,1);
			obj.Ez=gridquantity3D(obj.fullpath,'/data/fields/Ez',sum(obj.nnr)+1, obj.nthet, obj.nz+1,1);

			% Spline data
			kr=obj.femorder(2)+1;
			obj.knotsr=augknt(obj.rgrid,kr);
			kz=obj.femorder(1)+1;
			obj.knotsz=augknt(obj.zgrid,kz);
			obj.celltype=h5read(obj.fullpath,'/data/input.00/geometry/ctype')';
			obj.linked_s=h5read(obj.fullpath,'/data/input.00/geometry/linked_s');
			obj.bsplinetype=h5read(obj.fullpath,'/data/input.00/geometry/bsplinetype');
			obj.bsplinetype=reshape(obj.bsplinetype,length(obj.knotsr)-kr,length(obj.knotsz)-kz)';
			obj.geomweight = h5read(obj.fullpath,'/data/input.00/geometry/geomweight');
			obj.geomweight= reshape(obj.geomweight,length(obj.rgrid),length(obj.zgrid),[]);
			obj.dirichletweight = h5read(obj.fullpath,'/data/input.00/geometry/dirichletweight');
			obj.dirichletweight= reshape(obj.dirichletweight,length(obj.rgrid),length(obj.zgrid),[]);
			obj.gtilde = h5read(obj.fullpath,'/data/input.00/geometry/gtilde');
			obj.gtilde= reshape(obj.gtilde,length(obj.rgrid),length(obj.zgrid),[]);
			obj.CellVol= reshape(h5read(obj.fullpath,'/data/fields/volume'),length(obj.knotsr)-kr,length(obj.knotsz)-kz);
			zvol=fnder(spmak(obj.knotsz,ones(1,length(obj.knotsz)-kz)), -1 );
			rvol=fnder(spmak(obj.knotsr,2piobj.knotsr(kr:length(obj.knotsr)-kr+1)'), -1 );
			ZVol=diff(fnval(zvol,obj.knotsz));
			RVol=diff(fnval(rvol,obj.knotsr));
			CellVolfull=double(1/obj.nthet)RVol(kr:end-kr+1)ZVol(kz:end-kz+1)';
			CellVolfull=padarray(CellVolfull,double([kr-1,kz-1]),'replicate','post');
			obj.CellVol(obj.CellVol./CellVolfull<0.005)=0; % We remove the volumes which are too small and could lead to noise

			%obj.phi=splinequantity3D(obj.fullpath,'/data/fields/phi', obj.knotsr, obj.thetgrid, obj.knotsz, obj.femorder, 1, obj.geomweight(:,:,1), -1);
			obj.phireal=modequantity(obj.fullpath,'/data/fields/phireal', obj.knotsr, obj.knotsz, obj.nthet/2+1, obj.femorder, 1, obj.geomweight(:,:,1));
			obj.phiimag=modequantity(obj.fullpath,'/data/fields/phiimag', obj.knotsr, obj.knotsz, obj.nthet/2+1, obj.femorder, 1, obj.geomweight(:,:,1));
			%obj.rhs=splinequantity3D(obj.fullpath,'/data/fields/rhs', obj.knotsr, obj.thetgrid, obj.knotsz, obj.femorder, 1, obj.geomweight(:,:,1), -1);
			%obj.rhsreal=modedensity(obj.fullpath,'/data/fields/rhsreal', obj.knotsr, obj.knotsz, obj.nthet/2+1, obj.femorder, obj.CellVol, 1, obj.geomweight(:,:,1));
			%obj.rhsimag=modedensity(obj.fullpath,'/data/fields/rhsimag', obj.knotsr, obj.knotsz, obj.nthet/2+1, obj.femorder, obj.CellVol, 1, obj.geomweight(:,:,1));
			obj.N=splinedensity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.thetgrid, obj.knotsz, obj.femorder, obj.CellVol, 1, obj.geomweight(:,:,1), 1);
			obj.fluidUR=splinevelocity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.thetgrid, obj.knotsz, obj.femorder, obj.vnorm, obj.geomweight(:,:,1), 2);
			obj.fluidUTHET=splinevelocity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.thetgrid, obj.knotsz, obj.femorder, obj.vnorm, obj.geomweight(:,:,1), 3);
			obj.fluidUZ=splinevelocity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.thetgrid, obj.knotsz, obj.femorder, obj.vnorm, obj.geomweight(:,:,1), 4);

			%% Load maxwellian source 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.rate=obj.maxwellsrce.frequencyobj.weight/(pi(diff(obj.maxwellsrce.rlim.^2))*diff(obj.maxwellsrce.zlim));
			obj.maxwellsrce.current=obj.maxwellsrce.frequencyobj.weightobj.qe;
			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];
			obj.neutcol.iom_cross_sec=zeros(500,3);
			obj.neutcol.iom_cross_sec(:,1)=logspace(log10(obj.neutcol.Eion+0.001),log10(5e4),size(obj.neutcol.iom_cross_sec,1));
			obj.neutcol.iom_cross_sec(:,2)=obj.sigmiopre(obj.neutcol.iom_cross_sec(:,1),true);
			obj.neutcol.iom_cross_sec(:,3)=abs([log(obj.neutcol.iom_cross_sec(2:end,2)./obj.neutcol.iom_cross_sec(1:end-1,2))...
			./log(obj.neutcol.iom_cross_sec(2:end,1)./obj.neutcol.iom_cross_sec(1:end-1,1)); 0]);

			catch
			obj.neutcol.io_cross_sec=[];
			obj.neutcol.iom_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

			%% Load spline boundaries
			try
			obj.spl_bound.nbsplines=h5readatt(obj.fullpath, '/data/input.00/geometry_spl','nbsplines');
			for i=1:obj.spl_bound.nbsplines
			splgroup=sprintf('/data/input.00/geometry_spl/%02d',i);
			obj.spl_bound.boundary(i).knots=h5read(obj.fullpath,sprintf('%s/knots',splgroup));
			obj.spl_bound.boundary(i).Dval=h5readatt(obj.fullpath,splgroup,'Dirichlet_val');
			obj.spl_bound.boundary(i).coefs=reshape(h5read(obj.fullpath,sprintf('%s/pos',splgroup)),2,[])';
			obj.spl_bound.boundary(i).order=h5readatt(obj.fullpath,splgroup,'order');
			obj.spl_bound.boundary(i).kind=h5readatt(obj.fullpath,splgroup,'kind');
			try
			obj.spl_bound.boundary(i).type=h5readatt(obj.fullpath,splgroup,'type');
			catch
			end
			obj.spl_bound.boundary(i).fun=spmak(obj.spl_bound.boundary(i).knots,obj.spl_bound.boundary(i).coefs');
			end
			catch
			obj.spl_bound.nbsplines=0;
			end

			%% load non ideal power supply parameters
			try
			obj.psupply.targetbias=h5readatt(obj.fullpath, '/data/input.00/psupply','targetbias');
			obj.psupply.expdens=h5readatt(obj.fullpath, '/data/input.00/psupply','expdens');
			obj.psupply.PSresistor=h5readatt(obj.fullpath, '/data/input.00/psupply','PSresistor');
			obj.psupply.geomcapacitor=h5readatt(obj.fullpath, '/data/input.00/psupply','geomcapacitor');
			obj.psupply.nbhdt=h5readatt(obj.fullpath, '/data/input.00/psupply','nbhdt');
			obj.psupply.biases=h5read(obj.fullpath, '/data/var0d/biases');
			obj.psupply.current=h5read(obj.fullpath, '/data/var0d/current');
			obj.psupply.tau=obj.psupply.PSresistorobj.psupply.geomcapacitorobj.psupply.expdens/obj.neutcol.neutdens;
			obj.psupply.active=true;
			obj.psupply.bdpos=h5read(obj.fullpath, '/data/input.00/psupply/bdpos');
			catch
			obj.psupply.active=false;
			end

			%% Read the main particles parameters
			obj.species=h5parts.empty(obj.nbspecies,0);
			if(readparts)
			obj.species(1)=h5parts(obj.fullpath,'/data/part',obj,obj.normalized);

			if(obj.nbspecies >1)
			for i=2:obj.nbspecies
			obj.species(i)=h5parts(obj.fullpath,sprintf('/data/part/%2d',i),obj,true);
			try
			obj.species_moments(i-1).N=splinedensity(obj.fullpath, sprintf('/data/part/%2d/moments',i), obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, 1, geomweight, 1);
			obj.species_moments(i-1).fluidUR=splinevelocity(obj.fullpath, sprintf('/data/part/%2d/moments',i), obj.knotsr, obj.knotsz, obj.femorder, obj.vnorm, geomweight, 2);
			obj.species_moments(i-1).fluidUTHET=splinevelocity(obj.fullpath, sprintf('/data/part/%2d/moments',i), obj.knotsr, obj.knotsz, obj.femorder, obj.vnorm, geomweight, 3);
			obj.species_moments(i-1).fluidUZ=splinevelocity(obj.fullpath, sprintf('/data/part/%2d/moments',i), obj.knotsr, obj.knotsz, obj.femorder, obj.vnorm, geomweight, 4);
			catch
			% No moments stored
			obj.species_moments(i-1)=0;
			end

			end
			end
			end

			obj.npart=h5read(obj.fullpath,'/data/var0d/nbparts');
			try
			obj.tpart = h5read(obj.fullpath,'/data/part/time');
			obj.nbparts = h5read(obj.fullpath,'/data/part/Nparts');
			catch
			obj.nbparts=h5read(obj.fullpath,'/data/var0d/nbparts');
			obj.tpart=obj.t2d;
			end
			end

			%------------------------------------------
			% Functions for accesing secondary simulation quantities
			function Atheta=Atheta(obj,R,Z)
			% returns the magnetic vector potential at position R,Z interpolated from stored Athet in h5 file
			% halflz=(obj.zgrid(end)+obj.zgrid(1))/2;
			% Atheta=0.5obj.B0(R-obj.width/pi*(obj.Rcurv-1)/(obj.Rcurv+1)...
			% .besseli(1,2piR/obj.width).cos(2pi(Z-halflz)/obj.width));
			Atheta=interp2(obj.rgrid,obj.zgrid,obj.Athet,R,Z,'spline');
			end

			function quantity=H(obj,indices)
			% computes the total energy for the main specie
			% for the particle with index indices{1} at timepart step indices{2}
			% which is time obj.timepart(indices{2})
			if strcmp(indices{1},':')
			p=1:obj.species(1).VR.nparts;% if nothing is defined we load all particles
			else
			p=indices{1};
			end
			if strcmp(indices{2},':')
			t=1:length(obj.tpart); %if nothing is defined all time steps are considered
			else
			t=indices{2};
			end
			% if track is true we look at specific particles with their
			% index and follow them in time
			% if it is false we just care about the distribution function
			% and specific particles can have different positions in the
			% resulting array for each timestep
			if size(indices,1)>2
			track=indices{3};
			else
			track=false;
			end
			quantity=0.5obj.me(obj.species(1).VR(p,t,track).^2+obj.species(1).VTHET(p,t,track).^2+obj.species(1).VZ(p,t,track).^2)+obj.species(1).partepot(p,t,track);
			end

			function quantity=P(obj,indices)
			%P computes the canonical angular momentum for the main specie
			% for the particle with index indices{1} at timepart step indices{2}
			% which is time obj.timepart(indices{2})
			if strcmp(indices{1},':')
			p=1:obj.species(1).R.nparts;
			else
			p=indices{1};
			end
			if strcmp(indices{2},':')
			t=1:length(obj.tpart);
			else
			t=indices{2};
			end
			% if track is true we look at specific particles with their
			% index and follow them in time
			% if it is false we just care about the distribution function
			% and specific particles can have different positions in the
			% resulting array for each timestep
			if size(indices,1)>2
			track=indices{3};
			else
			track=false;
			end
			quantity=obj.species(1).R(p,t,track).(obj.species(1).VTHET(p,t,track)obj.me+sign(obj.qsim)obj.qeobj.Atheta(obj.species(1).R(p,t,track),obj.species(1).Z(p,t,track)));
			end

			function quantity=Vpar(obj,varargin)
			%Vpar Computes the parallel velocity for the main specie
			% for the particle with index indices{1} at timepart step indices{2}
			% which is time obj.timepart(indices{2})
			if(~iscell(varargin))
			indices=mat2cell(varargin);
			else
			indices=varargin;
			end
			if strcmp(indices{1},':')
			p=1:obj.species(1).R.nparts;
			else
			p=indices{1};
			end
			if strcmp(indices{2},':')
			t=1:length(obj.tpart);
			else
			t=indices{2};
			end
			% if track is true we look at specific particles with their
			% index and follow them in time
			% if it is false we just care about the distribution function
			% and specific particles can have different positions in the
			% resulting array for each timestep
			if size(indices,1)>2
			track=indices{3};
			else
			track=false;
			end
			Zp=obj.species(1).Z(p,t,track);% get the particle axial positon
			Rp=obj.species(1).R(p,t,track);% get the particle radial position

			% interpolate the magnetic field at the particle position
			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);
			% calculate the projection angle of the radial and axial
			% directions on the magnetic field line
			Costhet=Bzp./Bp;
			Sinthet=Brp./Bp;
			% calculate the actuale parallel velocity
			quantity=obj.species(1).VR(p,t,track).Sinthet+obj.species(1).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.species(1).R.nparts;
			else
			p=indices{1};
			end
			if strcmp(indices{2},':')
			t=1:length(obj.tpart);
			else
			t=indices{2};
			end
			% if track is true we look at specific particles with their
			% index and follow them in time
			% if it is false we just care about the distribution function
			% and specific particles can have different positions in the
			% resulting array for each timestep
			if size(indices,2)>2
			track=indices{3};
			else
			track=false;
			end

			% if gcs is true, gives the perpendicular velocity in the
			% guiding center system by substracting the EXB azimuthal
			% velocity
			% else gives the total perpendicular velocity
			if size(indices,2)>3
			gcs=indices{4};
			else
			gcs=false;
			end
			% get the particle position
			Zp=obj.species(1).Z(p,t,track);
			Rp=obj.species(1).R(p,t,track);
			% interpolate the magnetic field at the particle position
			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);
			% calculate the projecting angles
			Costhet=Bzp./Bp;
			Sinthet=Brp./Bp;
			Vdrift=zeros(size(Zp));

			if gcs
			% for each particle and each timestep
			% calculate the azimuthal ExB drift velocity
			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
			% calculate the perpendicular velocity
			quantity=sqrt((obj.species(1).VTHET(p,t,track)-Vdrift).^2+(obj.species(1).VR(p,t,track).Costhet-obj.species(1).VZ(p,t,track).Sinthet).^2);
			end

			function quantity=cyclphase(obj,varargin)
			%cyclphase Computes the cyclotronic phase for the main specie
			% for particles with indices{1} at time indices{2}

			if(~iscell(varargin))
			indices=mat2cell(varargin);
			else
			indices=varargin;
			end

			if strcmp(indices{1},':')
			p=1:obj.species(1).R.nparts;
			else
			p=indices{1};
			end
			if strcmp(indices{2},':')
			t=1:length(obj.tpart);
			else
			t=indices{2};
			end
			% if track is true we look at specific particles with their
			% index and follow them in time
			% if it is false we just care about the distribution function
			% and specific particles can have different positions in the
			% resulting array for each timestep
			if size(indices,2)>2
			track=indices{3};
			else
			track=false;
			end
			Zp=obj.species(1).Z(p,t,track);
			Rp=obj.species(1).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;

			% compute the projection of the perpendicular velocity in the
			% radial direction
			vr=(obj.species(1).VR(p,t,track).Costhet-obj.species(1).VZ(p,t,track).Sinthet);
			% Get the perpendicular velocity
			vperp=obj.Vperp(p,t,track,true);
			vr=vr(indices);
			vperp=vperp(indices);
			cospsi=vr./vperp;
			quantity=acos(cospsi);
			end

			function p=borderpoints(obj,subdiv)
			%borderpoints Return a cell array containing the curves
			%defining the boundary of the domain
			% for each boundary p(1,:) and p(2,:) give axial and radial position
			% for each boundary p(3,:) and p(4,:) give axial and radial normals

			%gw= contourc(obj.zgrid,obj.rgrid,obj.geomweight(:,:,1),[0 0])
			p=cell(0,0);
			if nargin<2
			subdiv=1;
			end

			ndiv=sum(subdiv);

			%outer cylinder
			if any(obj.geomweight(end,:,1)>=0)
			idp=ceil(length(obj.zgrid)/ndiv);
			imin=1;
			for j=1:length(subdiv)
			imax=min(imin+subdiv(j)*idp-1,length(obj.zgrid));
			p{end+1}=[obj.zgrid(imin:imax)';obj.rgrid(end)*ones(imax-imin+1,1)';
			zeros(imax-imin+1,1)';ones(imax-imin+1,1)'];
			imin=imax;
			end

			end

			%inner cylinder
			if any(obj.geomweight(1,:,1)>=0)
			idp=ceil(length(obj.zgrid)/ndiv);
			imin=1;
			for j=1:length(subdiv)
			imax=min(imin+subdiv(j)*idp-1,length(obj.zgrid));
			p{end+1}=[obj.zgrid(imin:imax)';obj.rgrid(1)*ones(imax-imin+1,1)';
			zeros(imax-imin+1,1)';-ones(imax-imin+1,1)'];
			imin=imax;
			end
			end

			if obj.walltype==2
			% We have an elliptic insert that we want to isolate
			gw=obj.ellipseborder;
			zpos=obj.zgrid(obj.zgrid<(min(gw(1,:))) \| obj.zgrid>(max(gw(1,:))));
			p{2}=[zpos,obj.rgrid(end)*ones(size(zpos))]';
			p{1}=[obj.zgrid';obj.rgrid(1)*ones(size(obj.zgrid))'];
			gw=obj.ellipseborder;
			p{3}=gw;
			elseif obj.walltype==9
			% extract all the walls
			gw=contourc(obj.zgrid,obj.rgrid,obj.geomweight(:,:,1),[0 0]);
			[x,y,~]=C2xyz(gw);
			for i=1:obj.spl_bound.nbsplines
			%subdiv=[4,1,2];
			if( obj.spl_bound.boundary(i).type>1)
			continue
			end
			s=linspace(obj.spl_bound.boundary(i).knots(1),obj.spl_bound.boundary(i).knots(end),5*length(obj.spl_bound.boundary(i).knots));
			positions=fnval(obj.spl_bound.boundary(i).fun,s);
			ndiv=sum(subdiv);
			idp=ceil(size(positions,2)/ndiv);
			imin=1;
			for j=1:length(subdiv)
			imax=min(imin+subdiv(j)*idp-1,size(positions,2));
			p{end+1}=[positions(1,imin:imax);positions(2,imin:imax)];
			imin=imax;
			end
			end
			elseif obj.walltype~=0
			% extract all the walls
			gw=contourc(obj.zgrid,obj.rgrid,obj.geomweight(:,:,1),[0 0]);
			[x,y,~]=C2xyz(gw);
			for i=1:length(x)
			%subdiv=[4,1,2];
			ndiv=sum(subdiv);
			idp=ceil(length(x{i})/ndiv);
			imin=1;
			for j=1:length(subdiv)
			imax=min(imin+subdiv(j)*idp-1,length(x{i}));
			p{end+1}=[x{i}(imin:imax);y{i}(imin:imax)];
			imin=imax;
			end
			end
			end
			% figure
			% for i=1:length(p)
			% plot(p{i}(1,:),p{i}(2,:))
			% hold on
			% end
			end

			function p=ellipseborder(obj)
			%ellipseborder returns the boundary points defining the
			%elliptic insert
			z=linspace(-0.998,0.998,600)*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 Gamma=Axialflux(obj,timestep,zpos, species_id)
			if nargin <4
			species_id=1;
			end
			% Computes the axial particle flux n*Uz at timestep timestep and axial position zpos
			if species_id ==1
			Gamma=obj.fluidUZ(:,zpos,timestep).*obj.N(:,zpos,timestep);
			else
			Gamma=obj.species_moments.fluidUZ(:,zpos,timestep).*obj.species_moments.N(:,zpos,timestep);
			end
			end

			function Gamma=Metallicflux(obj,timestep,subdiv)
			% Computes the particle flux at time obj.t2d(timestep) on the
			% metallic boundaries

			if nargin<3
			subdiv=1;
			end
			% We find the borderpoints
			p=obj.borderpoints(subdiv);
			gamma=cell(size(p));
			Nr=cell(size(p));
			Nz=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{i}=-interp2(obj.zgrid,obj.rgrid,obj.geomweight(:,:,3),bp(1,:),bp(2,:),'linear',0);
			Nz{i}=-interp2(obj.zgrid,obj.rgrid,obj.geomweight(:,:,2),bp(1,:),bp(2,:),'linear',0);
			norm=sqrt(Nr{i}.^2+Nz{i}.^2);
			norm(norm==0)=1;
			norm(isnan(norm))=1;
			Nr{i}=Nr{i}./norm;
			Nz{i}=Nz{i}./norm;

			else
			Nr{i}=bp(4,:);
			Nz{i}=bp(3,:);
			end

			gamma{i}=zeros(size(bp,2),length(timestep));
			end
			[z,r]=ndgrid(obj.zgrid,obj.rgrid);
			N=obj.N(:,:,timestep(1));
			n=griddedInterpolant(z,r,N'); % choose as function of species # (dep espece)
			uz=griddedInterpolant(z,r,obj.fluidUZ(:,:,timestep(1))');
			ur=griddedInterpolant(z,r,obj.fluidUR(:,:,timestep(1))');
			% we get the density and fluid velocities at the desired time
			% steps and interpolate them at the boundary position
			for j=1:length(timestep)
			n.Values=obj.N(:,:,timestep(j))'; % dependence esp
			%n.Values(obj.geomweight(:,:,1)<=0)=0;
			uz.Values=obj.fluidUZ(:,:,timestep(j))';
			%uz.Values(obj.geomweight(:,:,1)<=0)=0;
			ur.Values=obj.fluidUR(:,:,timestep(j))';
			%ur.Values(obj.geomweight(:,:,1)<=0)=0;
			for i=1:length(p)
			bp=p{i};
			gamma{i}(:,j)=n(bp(1:2,:)').(ur(bp(1:2,:)').Nr{i}'+uz(bp(1:2,:)').*Nz{i}');
			end
			end
			% return the boundary position p and the corresponding flux
			% gamma
			Gamma.p=p;
			Gamma.gamma=gamma;
			end

			function gamma_species = MetallicFlux_species(obj,timestep,subdiv)
			% Computes the particle flux at time obj.t2d(timestep) on the
			% metallic boundaries for species 2, supposedly the ions

			if nargin<3
			subdiv=1;
			end
			% We find the borderpoints
			p=obj.borderpoints(subdiv);
			gamma=cell(size(p));
			Nr=cell(size(p));
			Nz=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{i}=-interp2(obj.zgrid,obj.rgrid,obj.geomweight(:,:,3),bp(1,:),bp(2,:),'linear',0);
			Nz{i}=-interp2(obj.zgrid,obj.rgrid,obj.geomweight(:,:,2),bp(1,:),bp(2,:),'linear',0);
			norm=sqrt(Nr{i}.^2+Nz{i}.^2);
			norm(norm==0)=1;
			norm(isnan(norm))=1;
			Nr{i}=Nr{i}./norm;
			Nz{i}=Nz{i}./norm;
			else
			Nr{i}=bp(4,:);
			Nz{i}=bp(3,:);
			end

			gamma{i}=zeros(size(bp,2),length(timestep));
			end
			[z,r]=ndgrid(obj.zgrid,obj.rgrid);
			N=obj.species_moments.N(:,:,timestep(1));
			n=griddedInterpolant(z,r,N'); % choose as function of species # (dep espece)
			uz=griddedInterpolant(z,r,obj.species_moments.fluidUZ(:,:,timestep(1))');
			ur=griddedInterpolant(z,r,obj.species_moments.fluidUR(:,:,timestep(1))');
			% we get the density and fluid velocities at the desired time
			% steps and interpolate them at the boundary position
			for j=1:length(timestep)
			n.Values=obj.species_moments.N(:,:,timestep(j))'; % dependence esp
			uz.Values=obj.species_moments.fluidUZ(:,:,timestep(j))';
			ur.Values=obj.species_moments.fluidUR(:,:,timestep(j))';
			for i=1:length(p)
			bp=p{i};
			gamma{i}(:,j)=n(bp(1:2,:)').(ur(bp(1:2,:)').Nr{i}'+uz(bp(1:2,:)').*Nz{i}');
			end
			end
			% return the boundary position p and the corresponding flux
			% gamma
			gamma_species.p=p;
			gamma_species.gamma=gamma;


			end

			function [I, pos]=OutCurrents(obj,timestep, subdiv)
			% Computes the Outgoing currens at the simulation axial boundaries at timestep timestep
			% This is simply the surface integral of the axial flux
			if nargin<3
			subdiv=1;
			end
			flux=obj.Axialflux(timestep,[1 obj.nz+1]);
			Iz=squeeze(trapz(obj.rgrid,flux.obj.rgrid)2piobj.qsim/obj.weight);
			Iz(1,:)=-Iz(1,:);

			gamm=obj.Metallicflux(timestep, subdiv);
			Im=zeros(length(gamm.p),length(timestep));
			pos=cell(size(gamm.p));
			for i=1:length(gamm.p)
			p=gamm.p{i};
			pos{i}=p;
			flux=gamm.gamma{i}';
			for j=1:length(timestep)
			Im(i,j)=pi/2sum((p(2,1:end-1)+p(2,2:end)).(flux(j,2:end)+flux(j,1:end-1))...
			.*sqrt((p(1,2:end)-p(1,1:end-1)).^2+(p(2,2:end)-p(2,1:end-1)).^2));
			%AxialDensity = flux(j,:).* p(2,:);
			%Im(i,j) = 2pitrapz(p(1,:),AxialDensity);
			end
			end
			I=-cat(1,Iz,Im*obj.qsim/obj.weight);
			end

			function [I, pos]=OutCurrents_species(obj,timestep, subdiv)
			% Computes the Outgoing currens at the simulation axial boundaries at timestep timestep
			% This is simply the surface integral of the axial flux for
			% ions
			if nargin<3
			subdiv=1;
			end
			flux=obj.Axialflux(timestep,[1 obj.nz+1],2);
			Iz=squeeze(trapz(obj.species(2).rgrid,flux.obj.species(2).rgrid)2piobj.species(2).q);
			Iz(1,:)=-Iz(1,:);
			mflux = obj.MetallicFlux_species(timestep,subdiv);
			qe = abs(obj.species(2).q); % ions charge

			Im=zeros(length(mflux.p),length(timestep));
			pos=cell(size(mflux.p));

			for i=1:length(mflux.p)
			p=mflux.p{i};
			pos{i}=p;
			flux=mflux.gamma{i}';
			for j=1:length(timestep)
			Im(i,j)=pi/2sum((p(2,1:end-1)+p(2,2:end)).(flux(j,2:end)+flux(j,1:end-1))...
			.*sqrt((p(1,2:end)-p(1,1:end-1)).^2+(p(2,2:end)-p(2,1:end-1)).^2));
			%AxialDensity = flux(j,:).* p(2,:);
			%Im(i,j) = 2pitrapz(p(1,:),AxialDensity);
			end
			end
			I=cat(1,Iz,Im*qe);
			end

			function [pot] = PotentialWell(obj,fieldstep,fieldaligned)
			%PotentialWell Computes the potential well at the given timestep on the FEM grid points
			% interpolates the model data on rgrid and zgrid
			if nargin<3
			fieldaligned=false;
			end
			model=obj.potentialwellmodel(fieldstep);
			z=model.z;
			modpot=model.pot;

			if fieldaligned
			r=model.rathet;
			lvls=linspace(min(obj.rAthet(:)),max(obj.rAthet(:)),400);
			[Zmesh,Rmesh]=meshgrid(obj.zgrid,lvls);
			else
			r=model.r;
			[Zmesh,Rmesh]=meshgrid(obj.zgrid,obj.rgrid);
			end
			pot=zeros(size(Zmesh,2),size(Zmesh,1),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)
			%Ekin Computes the classical kinetic energy of particles indices{1} at
			% time obj.tpart(indices{2}) in Joules
			if(~iscell(varargin))
			indices=mat2cell(varargin);
			else
			indices=varargin;
			end
			if strcmp(indices{1},':')
			p=1:obj.species(1).R.nparts;
			else
			p=indices{1};
			end
			if strcmp(indices{2},':')
			t=1:length(obj.tpart);
			else
			t=indices{2};
			end
			% if track is true we look at specific particles with their
			% index and follow them in time
			% if it is false we just care about the distribution function
			% and specific particles can have different positions in the
			% resulting array for each timestep
			if size(indices,1)>2
			track=indices{3};
			else
			track=false;
			end
			Vr=obj.species(1).VR(p,t,track);
			Vthet= obj.species(1).VTHET(p,t,track);
			Vz=obj.species(1).VZ(p,t,track);
			Ekin=0.5obj.msim/obj.weight(Vr.^2+Vthet.^2+Vz.^2);
			end

			function sig=sigio(obj,E,init)
			%sigio returns the total ionisation cross-section in m^2
			% at energy E[eV]
			% init is only used during the loading of the h5 file
			if nargin <3
			init=false;
			end
			sig=zeros(size(E));

			if(~init &&( ~obj.neutcol.present \|\| isempty(obj.neutcol.io_cross_sec)))
			sig=zeros(size(E));
			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 sig=sigmio(obj,E)
			%sigmio returns the total ionisation cross-section for momentum exchange for the incoming electron in m^2
			% at energy E[eV]
			sig=zeros(size(E));
			if(~obj.neutcol.present \|\| isempty(obj.neutcol.iom_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.iom_cross_sec);
			end
			end
			end

			function sigm=sigmela(obj,E)
			%sigmela returns the elastic collision cross-section for momentum exchange for the incoming electron in m^2
			% at energy E[eV]
			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)
			%sigmela returns the elastic collision cross-section for the incoming electron in m^2
			% at energy E[eV]
			% if used this will give the frequency of elastic collisions
			E0=obj.neutcol.E0;
			chi=E./(0.25*E0+E);
			sig=(2chi.^2)./((1-chi).((1+chi).log((1+chi)./(1-chi))-2chi)).*obj.sigmela(E);
			end


			%------------------------------------------
			% Functions for plotting evolving quantities

			function line=displaysplbound(obj,ax,rescale,markers)
			%displaysplbound display on axis ax the boundary of the
			%simulation domain and the Dirichlet and Neumann walls defined
			%with spline curves
			if nargin<2
			ax=gca;
			end
			if nargin<3
			rescale=1;
			end
			if nargin<4
			markers=true;
			end
			hold on
			for i=1:obj.spl_bound.nbsplines
			if(obj.spl_bound.boundary(i).type==0)
			%continue
			end
			knots=obj.spl_bound.boundary(i).knots(1:end);
			coeffs=obj.spl_bound.boundary(i).coefs'*rescale;
			pp=spmak(knots,coeffs);
			sizec=size(coeffs,2);
			order=length(knots)-sizec;
			s=linspace(knots(order),knots(sizec+1),max(200,3*length(knots)));
			fittedpos=fnval(pp,s);
			line=plot(fittedpos(1,:),fittedpos(2,:),'-','linewidth',2);
			%ine=plot(coeffs(1,:),coeffs(2,:),'-','linewidth',2);
			if markers
			plot(coeffs(1,:),coeffs(2,:),'rx','markersize',14)
			for j=1:size(coeffs,2)
			text(coeffs(1,j),coeffs(2,j),sprintf('%i',j),'fontsize',14)
			end
			end
			end
			end

			function f=displayrprofile(obj,t,zpos,init)
			% plot the initial and final radial profile at the axial center of the simulation space
			% also plot the azimuthal fluid rotation frequency profile
			% t: time index considered
			% zpos: axial position index
			% init: initial time considered for comparison
			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(mean(obj.N(:,zpos,t),3),2);

			%inverse of radius
			Rinv=1./obj.rgrid;
			Rinv(isnan(Rinv))=0;
			%azimuthal velocity and azimuthal rotation frequency in m/s and
			%1/s
			vthet=mean(obj.fluidUTHET(:,zpos,t),3);
			omegare=(vthet.*Rinv);
			% plot the initial density
			if(init)
			plot(obj.rgrid,mean(obj.N(:,zpos,tinit),2),'bx-','DisplayName',sprintf('t=%1.2f[ns]',obj.t2d(tinit)*1e9),'linewidth',lw)
			end
			hold on
			%plot the time averaged current density
			p=plot(obj.rgrid,locdens,'r-','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','Northeast')
			% limit the axis to the simulation domain
			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();
			% plot the metallic walls for a constant radius coaxial
			% configuration
			if obj.conformgeom
			plot(obj.rgrid(1)*[1 1],ylimits,'k--')
			plot(obj.rgrid(end)*[1 1],ylimits,'k--')
			else
			geomw=obj.geomweight(:,floor(mean(zpos)),1)>0;
			changes=find(abs(diff(geomw))>0);
			pos=(obj.rgrid(changes)+obj.rgrid(changes+1))./2;
			for i=1:length(pos)
			plot(pos(i)*[1 1],ylimits,'k--')
			end
			end
			p.YData(~geomw)=NaN;
			legend(p)

			%yyaxis right
			% plot the azimuthal fluid rotation frequency profile
			%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]',mean(obj.zgrid(zpos))))
			obj.savegraph(f,sprintf('%srProf%d_%d',obj.name,mean(zpos),mean(t)),[15 10]);
			end

			function displayenergy(obj)
			%% Plot the time evolution of the system energy and number of simulated macro particles
			tmin=2;
			tmax=length(obj.ekin);
			f=figure('Name', sprintf('%s Energy',obj.name));
			subplot(1,2,1)
			semilogy(obj.t0d(tmin:tmax)*1e9,obj.ekin(tmin:tmax),'-',...
			obj.t0d(tmin:tmax)*1e9,obj.epot(tmin:tmax),':',...
			obj.t0d(tmin:tmax)*1e9,obj.etot(tmin:tmax),'-.', ... 'linewidth',4)
			obj.t0d(tmin:tmax)*1e9,obj.etot0(tmin:tmax),'--','linewidth',4)
			%obj.t0d(tmin:tmax),obj.ekin(tmin:tmax)-obj.epot(tmin:tmax),'--',
			legend('E_{kin}', 'E_{pot}', 'E_{tot}','E_{ref}','location','east')
			xlabel('Time [ns]')
			ylabel('Energies [J]')
			grid on
			xlim([obj.t0d(tmin) obj.t0d(tmax)]*1e9);
			xlimits=xlim();
			%ylim([1e-7 1e-3])
			set(gca,'fontsize',14)
			subplot(1,2,2)
			try
			semilogy(obj.t0d(tmin:tmax)*1e9,abs(obj.eerr(tmin:tmax)./obj.etot0(tmin:tmax)),'-','linewidth',1.5)
			catch
			semilogy(obj.t0d(tmin:tmax)*1e9,abs(obj.eerr(tmin:tmax)/obj.etot(2)),'-','linewidth',1.5)
			end
			hold on
			xlabel('Time [ns]')
			ylabel('Total energy relative error [-]')
			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
			ylimits=ylim;
			for i=1:length(obj.restarttimes)
			plot(obj.restarttimes(i)[1 1]1e9,ylimits,'k--')
			end
			set(gca,'fontsize',14)
			Position=f.Position;
			Position(3)=1.5*Position(3);
			f.Position=Position;
			obj.savegraph(f,sprintf('%s/%sEnergy',obj.folder,obj.name));
			end

			function out=displaycharge(obj,scalet)
			%% Plot the time evolution of the system charge of electrons
			tmin=1;
			tmax=length(obj.t0d);
			if nargin<2
			scalet =false;
			end

			f=figure('Name', sprintf('%s Charge',obj.name));
			ax=f.CurrentAxes;
			if isempty(ax)
			ax=axes(f);
			end
			if scalet
			if obj.neutcol.present
			vexb0=(obj.Ez(:,:,1).obj.Br'-obj.Er(:,:,1).obj.Bz')./(obj.B'.^2);
			vexb0(obj.geomweight(:,:,1)<=0)=0;
			E=0.5obj.msim/obj.weightmean(abs(vexb0(:)))^2/obj.qe;
			taucol=1/(obj.neutcol.neutdensmean(abs(vexb0(:)))(obj.sigio(E)+obj.sigmela(E)+obj.sigmio(E)));
			try
			Sio_S=1e17(obj.neutcol.neutdensmean(abs(vexb0(:)))obj.sigio(E))/(obj.maxwellsrce.frequencyobj.weight/(pi(obj.maxwellsrce.rlim(2)^2-obj.maxwellsrce.rlim(1)^2)diff(obj.maxwellsrce.zlim)))
			catch
			end
			tlabel='t/\tau_d [-]';
			else
			taucol=2*pi/obj.omece;
			tlabel='t/\tau_ce [-]';
			end
			else
			taucol=1e-9;
			tlabel='t [ns]';
			end
			plot(ax,obj.t0d(tmin:tmax)/taucol,abs(obj.npart(tmin:tmax)*obj.qsim),'linewidth',2)
			hold on
			ylabel(ax,'Total charge [C]')
			xlabel(ax,tlabel)
			grid on
			set(ax,'fontsize',12)

			out.time = obj.t0d(tmin:tmax)/taucol;
			out.charge = abs(obj.npart(tmin:tmax)*obj.qsim);
			end

			function displaySimParticles(obj)
			%% Plot the time evolution of the number of simulated markers in the main specie
			f=figure('Name', sprintf('%s Trapped particles',obj.name));
			plot(obj.t0d,obj.npart,'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
			% the larmor radius is calculated by considering that the
			% initial perpendicular velocity \approx the ExB velocity
			if time2d>0
			Er=obj.Er(:,:,time2d);
			Ez=obj.Ez(:,:,time2d);
			else
			Er=obj.Erxt(:,:,1);
			Ez=obj.Ezxt(:,:,1);
			end

			rl=abs(obj.me/obj.qe(-Er.obj.Bz'+Ez.*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-','linewidth',3)
			if time2d>0
			n=obj.density(:,:,time2d);
			maxN=max(n (:));
			n=n/maxN*mean(rl(:));
			contour(obj.zgrid,obj.rgrid,n,linspace(0,1,6)*mean(rl(:)),'r:','linewidth',3)
			end
			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.species(1).R.nt>=2)
			%tstart=obj.species(1).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.species(1).R.nparts);
			Hloc=H(obj,{1:obj.nbparts(tstart),tstart,false});
			h1=histogram(Hloc,20,'BinLimits',[min(Hloc(:)) max(Hloc(:))],'DisplayName',sprintf("t=%2.3d [ns]",obj.tpart(tstart)*1e9));
			hold on
			Hloc=H(obj,{1:partsmax,obj.species(1).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(tstart),tstart,false});
			h2=histogram(Ploc,50,'BinLimits',[min(Ploc(:)) max(Ploc(:))],'DisplayName',sprintf("t=%2.3d [ns]",obj.tpart(tstart)*1e9));
			hold on
			Ploc=P(obj,{1:partsmax,obj.species(1).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.95h2.BinLimits(1) 1.05h2.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));
			vr2=obj.species(1).VR(:,:,false);
			vr2=mean(vr2.^2,1)-mean(vr2,1).^2;
			vz2=obj.species(1).VZ(:,:,false);
			vz2=mean(vz2.^2,1)-mean(vz2,1).^2;
			vthet2=obj.species(1).VTHET(:,:,false);
			vthet2=mean(vthet2.^2,1)-mean(vthet2,1).^2;
			plot(obj.tpart,0.5obj.mevr2/obj.qe,'displayname','T_r')
			hold on
			plot(obj.tpart,0.5obj.mevz2/obj.qe,'displayname','T_z')
			plot(obj.tpart,0.5obj.mevthet2/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/1e3,obj.potout/1e3,obj.Rcurv))
			legend
			grid
			obj.savegraph(f,sprintf('%s/%s_partstemp',obj.folder,obj.name));

			end

			function displayCurrentsevol(obj,timesteps,species_id)
			% Computes and display the time evolution of the outgoing currents on each domain boundary
			% at timesteps timesteps
			P=obj.neutcol.neutdensobj.kb300/100;
			if nargin<2
			timesteps=1:length(obj.t2d);
			end
			if nargin<3
			species_id=1;
			end
			if species_id ==1
			currents=obj.OutCurrents(timesteps);
			currents = currents/P;
			else
			currents=obj.OutCurrents_species(timesteps);
			currents = currents/P;
			end
			f=figure('Name',sprintf('%s Currents',obj.name));
			if(obj.B(1,1)>obj.B(end,1))
			lname='HFS';
			rname='LFS';
			else
			lname='LFS';
			rname='HFS';
			end
			plot(obj.t2d(timesteps),currents(1,:),'Displayname',lname,'linewidth',1.8);
			hold on
			plot(obj.t2d(timesteps),currents(2,:),'Displayname',rname,'linewidth',1.8);
			plot(obj.t2d(timesteps),currents(3,:),'Displayname','outer cylinder','linewidth',1.8);
			plot(obj.t2d(timesteps),currents(4,:),'Displayname','inner cylinder','linewidth',1.8);
			if size(currents,1)>=5
			for ii=5:size(currents,1)
			plot(obj.t2d(timesteps),currents(ii,:),'Displayname',sprintf('metallic wall %i',ii),'linewidth',1.8);
			end
			end
			legend('location','Northeast')
			xlabel('time [s]')
			ylabel('I [A/mbar]')
			grid on
			set(gca,'fontsize',12)
			title(sprintf('\\phi_b-\\phi_a=%.2g kV, R=%.1f',(obj.potout-obj.potinn)/1e3,obj.Rcurv))
			obj.savegraph(f,sprintf('%s/%s_outCurrents',obj.folder,obj.name),[16 12]);
			end

			function displayChargeLossevol(obj,timesteps,toptitle,scalet,dens)
			% Computes and display the time evolution of the outgoing currents on each domain boundary
			% at time obj.t2d(timesteps)
			%scalet=true scales the time by the ellastic collision
			%frequency
			%dens = true plot the time evolution of the maximum electron
			%density in the simulation domain otherwise plot the total
			%number of electrons in the domain
			if nargin<2
			timesteps=1:length(obj.t2d);
			end
			if nargin<4
			scalet=true;
			end
			if nargin <5
			dens=true;
			end

			if scalet
			if obj.neutcol.present
			vexb0=(obj.Ez(:,:,1).obj.Br'-obj.Er(:,:,1).obj.Bz')./(obj.B'.^2);
			vexb0(obj.geomweight(:,:,1)<=0)=0;
			E=0.5obj.msim/obj.weightmean(abs(vexb0(:)))^2/obj.qe;
			taucol=1/(obj.neutcol.neutdensmean(abs(vexb0(:)))(obj.sigio(E)+obj.sigmela(E)+obj.sigmio(E)));
			try
			Sio_S=1e17(obj.neutcol.neutdensmean(abs(vexb0(:)))obj.sigio(E))/(obj.maxwellsrce.frequencyobj.weight/(pi(obj.maxwellsrce.rlim(2)^2-obj.maxwellsrce.rlim(1)^2)diff(obj.maxwellsrce.zlim)))
			catch
			end
			tlabel='t/\tau_d [-]';
			else
			taucol=2*pi/obj.omece;
			tlabel='t/\tau_ce [-]';
			end
			else
			taucol=1e-9;
			tlabel='t [ns]';
			end

			if dens
			N=obj.N(:,:,timesteps);
			geomw=obj.geomweight(:,:,1);
			geomw(geomw<0)=0;
			geomw(geomw>0)=1;
			N=N.*geomw;

			nmax=squeeze(max(max(N,[],1),[],2));
			tn=(obj.t2d(timesteps));
			nlabel='n_{e,max} [m^{-3}]';
			ndlabel='n_{e,max}';
			else
			t0dst=find(obj.t0d>=obj.t2d(timesteps(1)),1,'first');
			t0dlst=find(obj.t0d<=obj.t2d(timesteps(end)),1,'last');
			tn=obj.t0d(t0dst:t0dlst);
			nmax=obj.npart(t0dst:t0dlst)*obj.weight;
			nlabel='Nb e^-';
			ndlabel='Nb e^-';
			end


			[currents,pos]=obj.OutCurrents(timesteps);
			P=obj.neutcol.neutdensobj.kb300/100;% pressure at room temperature in mbar
			currents=currents/P;
			f=figure('Name',sprintf('%s Charges',obj.name));
			% Plot the evolution of nb of particles
			yyaxis right
			p=plot(tn/taucol,nmax,'b-.','linewidth',1.8,'Displayname',ndlabel);
			ylabel(nlabel)
			ax=gca;
			ax.YAxis(2).Color=p.Color;
			ylim([0 inf])

			if(obj.B(1,1)>obj.B(end,1))
			lname='HFS';
			rname='LFS';
			else
			lname='LFS';
			rname='HFS';
			end

			yyaxis left
			mincurr=max(currents(:))*5e-3;
			if (max(currents(1,:)>mincurr))
			plot(obj.t2d(timesteps)/taucol,currents(1,:),'r:','Displayname',lname,'linewidth',1.8);
			end
			hold on
			if (max(currents(2,:)>mincurr))
			plot(obj.t2d(timesteps)/taucol,currents(2,:),'r--','Displayname',rname,'linewidth',1.8);
			end
			if (max(currents(3,:)>mincurr))
			plot(obj.t2d(timesteps)/taucol,currents(3,:),'r-','Displayname','outer cylinder','linewidth',1.8);
			end
			if (max(currents(4,:)>mincurr))
			plot(obj.t2d(timesteps)/taucol,currents(4,:),'Displayname','inner cylinder','linewidth',1.8);
			end
			if (size(currents,1)>=5 && max(currents(5,:)>mincurr))
			plot(obj.t2d(timesteps)/taucol,currents(5,:),'r-','Displayname','ellipse','linewidth',1.8);
			end
			xlabel(tlabel)
			ylabel('I/p_n [A/mbar]')
			grid on
			set(gca,'fontsize',12)
			ax.YAxis(1).Color='red';

			legend('Orientation','horizontal','location','south','numcolumns',3)

			if nargin <3
			title(sprintf('\\phi_b-\\phi_a=%.2g kV, B=%f T',(obj.potout-obj.potinn)/1e3,max(obj.B(:))))
			elseif ~isempty(toptitle)
			title(toptitle)
			end
			obj.savegraph(f,sprintf('%s/%s_ChargeEvol%i%i',obj.folder,obj.name,scalet,dens),[16 14]);
			end


			function out = displaytotcurrevol(obj,timesteps,species_id,all_cur_id, toptitle,scalet,dens,subdiv,nmean)
			% Computes and display the time evolution of the outgoing
			% currents at time obj.t2d(timesteps)
			%scalet=true scales the time by the ellastic collision
			%frequency
			%dens = true plot the time evolution of the maximum electron
			%density in the simulation domain otherwise plot the total
			%number of electrons in the domain
			% also plot in a subplot the color coded boundary corresponding
			% to each current
			if nargin<2
			timesteps=1:length(obj.t2d);
			end
			if nargin<3
			species_id =1;
			end
			if nargin<4
			all_cur_id = 0;
			end
			if nargin<5
			toptitle="";
			end
			if nargin<6
			scalet=true;
			end
			if nargin <7
			dens=true;
			end
			if nargin<8
			subdiv=1;
			end
			if nargin<9
			nmean=1;
			end

			if scalet
			if obj.neutcol.present
			vexb0=(obj.Ez(:,:,1).obj.Br'-obj.Er(:,:,1).obj.Bz')./(obj.B'.^2);
			vexb0(obj.geomweight(:,:,1)<=0)=0;
			E=0.5obj.msim/obj.weightmean(abs(vexb0(:)))^2/obj.qe
			taucol=1/(obj.neutcol.neutdensmean(abs(vexb0(:)))(obj.sigio(E)+obj.sigmela(E)+obj.sigmio(E)))
			try
			Sio_S=1e17(obj.neutcol.neutdensmean(abs(vexb0(:)))obj.sigio(E))/(obj.maxwellsrce.frequencyobj.weight/(pi(obj.maxwellsrce.rlim(2)^2-obj.maxwellsrce.rlim(1)^2)diff(obj.maxwellsrce.zlim)));
			catch
			end
			tlabel='t/\tau_d [-]';
			else
			taucol=2*pi/obj.omece;
			tlabel='t/\tau_ce [-]';
			end
			else
			taucol=1e-9;
			tlabel='t [ns]';
			end

			if dens
			if species_id==1
			N=obj.N(:,:,timesteps);
			else
			N=obj.species_moments.N(:,:,timesteps);
			end
			geomw=obj.geomweight(:,:,1);
			geomw(geomw<0)=0;
			geomw(geomw>0)=1;
			N=N.*geomw;
			%[~,idl]=max(N(:,:,end),[],'all','linear');
			%[ir,iz]=ind2sub(size(geomw),idl);
			%nmax=squeeze(max(max(N,[],1),[],2));
			tn=(obj.t2d(timesteps));
			nmax=zeros(2,length(tn));
			nrhalf= floor(0.5*length(obj.rgrid));%find(obj.rgrid>0.07);

			nmax(1,:)=squeeze(max(max(N(1:nrhalf,:,:),[],1),[],2));
			nmax(2,:)=squeeze(max(max(N(nrhalf+1:end,:,:),[],1),[],2));


			nlabel='n_{e,max} [m^{-3}]';
			ndlabel='n_{e,max}';
			else
			t0dst=find(obj.t0d>=obj.t2d(timesteps(1)),1,'first');
			t0dlst=find(obj.t0d<=obj.t2d(timesteps(end)),1,'last');
			tn=obj.t0d(t0dst:t0dlst);
			nmax=obj.npart(t0dst:t0dlst)*obj.weight;
			nlabel='Nb e^-';
			ndlabel='Nb e^-';
			end

			if species_id ==1
			[currents,pos]=obj.OutCurrents(timesteps,subdiv);
			P=obj.neutcol.neutdensobj.kb300/100;% pressure at room temperature in mbar
			currents=currents/P;
			else % for ionic currents
			[currents,pos]=obj.OutCurrents_species(timesteps,subdiv);
			currents=currents/obj.species(species_id).weight*obj.species(1).weight;
			P=obj.neutcol.neutdensobj.kb300/100;% pressure at room temperature in mbar
			currents=currents/P;
			end
			% hack to sum ionic + electronic current
			if all_cur_id==1
			[el_curr,pos]=obj.OutCurrents(timesteps,subdiv);
			P=obj.neutcol.neutdensobj.kb300/100;% pressure at room temperature in mbar
			out.elec_currents=el_curr/P;

			[ion_curr,pos]=obj.OutCurrents_species(timesteps,subdiv);
			P=obj.neutcol.neutdensobj.kb300/100;% pressure at room temperature in mbar
			out.ion_currents=ion_curr/P;

			totCurr_el = sum(out.elec_currents(1:end-1,:),1);
			totCurr_ion = sum(out.ion_currents,1);

			TotalCurrent = totCurr_el+totCurr_ion;
			out.temps = tn/taucol;
			end
			% plot
			f=figure('Name',sprintf('%s Charges',obj.name));
			tiledlayout(2,1)
			nexttile
			% Plot the evolution of nb of particles
			yyaxis right
			for i=1:size(nmax,1)
			p=plot(tn/taucol,nmax(i,:),'b','linewidth',2.2,'Displayname',sprintf('%s, %d',ndlabel,i));
			hold on
			end
			ylabel(nlabel)
			axl=gca;
			axl.YAxis(2).Color=p.Color;
			ylim([0 inf])

			if(obj.B(1,1)>obj.B(end,1))
			lname='HFS';
			rname='LFS';
			else
			lname='LFS';
			rname='HFS';
			end

			yyaxis( 'left');
			map=colormap(lines);
			set(axl,'linestyleorder',{'-',':','--','*','+'},...
			'ColorOrder',map(2:7,:), 'NextPlot','replacechildren')
			p(1)=plot(axl,obj.t2d(timesteps)/taucol,movmean(currents(1,:),nmean),'Displayname',lname,'linewidth',1.8);
			hold on
			p(2)=plot(axl,obj.t2d(timesteps)/taucol,movmean(currents(2,:),nmean),'Displayname',rname,'linewidth',1.8);
			% Plot the currents
			for i=3:size(currents,1)
			p(i)=plot(axl,obj.t2d(timesteps)/taucol,movmean(currents(i,:),nmean),'Displayname',sprintf('border %i',i-2),'linewidth',1.8);
			end

			if all_cur_id ==0
			plot(axl,obj.t2d(timesteps)/taucol,movmean(sum(currents(:,:),1,'omitnan'),nmean),'k-','Displayname','total','linewidth',1.8);
			else
			% plot(axl,obj.t2d(timesteps)/taucol,movmean(sum(out.elec_currents + out.ion_currents,1,'omitnan'),nmean),'k-','Displayname','total','linewidth',1.8);
			plot(axl,obj.t2d(timesteps)/taucol,movmean(TotalCurrent,nmean),'k-','Displayname','total','linewidth',1.8);
			hold on
			p(size(currents,1))=plot(axl,obj.t2d(timesteps)/taucol,movmean(out.ion_currents(end,:),nmean),'Displayname',sprintf('border %i',i-2),'linewidth',1.8);
			end
			hold on
			xlabel(tlabel)
			ylabel('I/p_n [A/mbar]')
			grid on
			set(gca,'fontsize',12)
			ax.YAxis(1).Color='black';

			%legend('Orientation','horizontal','location','north','numcolumns',3)

			if ~isempty(toptitle)
			title(toptitle)
			end

			ax2=nexttile;
			geomw=obj.geomweight(:,:,1);
			geomw(geomw<=0)=0;
			geomw(geomw>0)=NaN;
			[c1,hContour]=contourf(ax2,obj.zgrid1000,obj.rgrid1000,geomw, [0 0]);
			hold on
			drawnow;
			grid on;

			for i=1:length(pos)
			plot(ax2,pos{i}(1,:)1000,pos{i}(2,:)1000,'linestyle',p(i+2).LineStyle,...
			'color',p(i+2).Color,'marker',p(i+2).Marker,...
			'displayname',sprintf('border %i',i),'linewidth',1.8)
			hold on
			end
			title('Domain')
			plot(ax2,ones(size(obj.rgrid))obj.zgrid(1)1000,obj.rgrid*1000,'linestyle',p(1).LineStyle,...
			'color',p(1).Color,'marker',p(1).Marker,...
			'displayname',lname,'linewidth',1.8)
			plot(ax2,ones(size(obj.rgrid))obj.zgrid(end)1000,obj.rgrid*1000,'linestyle',p(2).LineStyle,...
			'color',p(2).Color,'marker',p(2).Marker,...
			'displayname',rname,'linewidth',1.8)
			xlabel('z [mm]')
			ylabel('r [mm]')
			grid on
			set(gca,'fontsize',12)
			hFills=hContour.FacePrims;
			[hFills.ColorType] = deal('truecoloralpha'); % default = 'truecolor'
			try
			hFills(1).ColorData = uint8([150;150;150;255]);
			for idx = 2 : numel(hFills)
			hFills(idx).ColorData(4) = 0; % default=255
			end
			catch
			end
			%legend('Orientation','horizontal','location','north','numcolumns',4)


			fprintf('mean total current: %f [A/mbar]\n',mean(sum(currents(:,max(1,size(currents,2)-30):end),1,'omitnan')));

			%if nargin <3
			% sgtitle(sprintf('\\phi_b-\\phi_a=%.2g kV, B=%f T',(obj.potout-obj.potinn)*obj.phinorm/1e3,mean(obj.B(:))))
			%elseif ~isempty(toptitle)
			% sgtitle(toptitle)
			%end
			if length(subdiv)>1
			obj.savegraph(f,sprintf('%s/%s_totIEvol%i%i_div',obj.folder,obj.name,scalet,dens),[16 14]);
			else
			obj.savegraph(f,sprintf('%s/%s_totIEvol%i%i',obj.folder,obj.name,scalet,dens),[16 14]);
			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(obj.rgrid-rpos));
			end
			crpos=obj.rgrid(rpos);
			id=find(timestep==0);
			timestep(id)=[];
			n=obj.density(:,:,timestep);
			if(~isempty(timestep==0))
			N0=zeros(obj.density.nr+1,obj.density.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,1e3crpos))
			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 displaypotentialwellrprof(obj,timestep,zpos)
			% Display the potential well radial profile at the axial
			% position zpos and time timestep
			if iscell(timestep)
			timestep=cell2mat(timestep);
			end
			f=figure('Name',sprintf('%s Potential well',obj.name));

			Pot=obj.PotentialWell(timestep);
			Potrprof=Pot(zpos,:)
			z=obj.zgrid(zpos);

			plot(obj.rgrid,Potrprof)
			ylabel('Potential [eV]')
			xlabel('r [m]')
			xlim([obj.rgrid(1) obj.rgrid(end)])
			title(sprintf('Potential well t=%1.2f [ns] z=%1.2f [mm]',obj.t2d(timestep)1e9,1e3z))

			obj.savegraph(f,sprintf('%s/%s_well1Dr_%d',obj.folder,obj.name,zpos));
			end

			function displayVdistribRThetZ(obj,timestep, rpos, zpos)
			%displayVdistribRThetZ plot the velocity distribution function
			% in m/s
			%extracted from the markers at position window from rpos(1)
			%rpos(end) and zpos(1) to zpos(end)
			% and at time obj.tpart(timestep)
			%rpos and zpos are given as grid indices
			if(obj.species(1).R.nt>=2)
			if nargin<2
			timesteppart=length(obj.tpart);
			timestep=timesteppart;
			else
			timesteppart=timestep;
			end
			if nargin<3 \|\| isempty(rpos)
			rpos=1:length(obj.rgrid);
			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)
			zpos=1:length(obj.zgrid);
			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.species(1).R.nparts);
			R=obj.species(1).R(1:nbp,1,false);
			Z=obj.species(1).Z(1:nbp,1,false);
			Vr=obj.species(1).VR(1:nbp,1,false);
			Vz=obj.species(1).VZ(1:nbp,1,false);
			Vthet=obj.species(1).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.species(1).R.nparts);
			Rend=obj.species(1).R(1:nbp,timesteppart,false);
			Zend=obj.species(1).Z(1:nbp,timesteppart,false);
			Vrend=obj.species(1).VR(1:nbp,timesteppart,false);
			Vzend=obj.species(1).VZ(1:nbp,timesteppart,false);
			Vthetend=obj.species(1).VTHET(1:nbp,timesteppart,false);
			ids=Rend>=rspan(1) & Rend<=rspan(2) & Zend>=zspan(1) & Zend<=zspan(2);
			nbtot=sum(ids)
			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));

			[~,time2did]=min(abs(obj.t2d-obj.tpart(timestep)));


			subplot(1,4,1);
			obj.dispV(Vr,Vrend,'V_r [m/s]',[1,timesteppart])

			[~,time2did]=min(abs(obj.t2d-obj.tpart(timestep)));
			if length(rpos)==1
			vexb=-obj.Er(rpos,zpos,time2did)/obj.Bz(zpos,rpos)';
			vexb=mean(vexb(:));

			if ~isempty(obj.neutcol.ela_cross_sec) % plot the radial drift velocity as nu_dE_r/(B\Omega_c)
			vdr=obj.neutcol.neutdensobj.sigmela(vexb^2obj.me0.5/obj.qe)vexb*-obj.Er(rpos,zpos,time2did)...
			./(obj.B(zpos,rpos)'.obj.B(zpos,rpos)'obj.qe/obj.me);
			vdr=mean(vdr(:));
			ylimits=ylim;
			plot(vdr*[1 1],ylimits,'k--','displayname',sprintf('V_{d,pred}=%1.2g [m/s]',vdr))
			end
			end


			subplot(1,4,2);
			obj.dispV(Vthet,Vthetend,'V\theta [m/s]',[1,timesteppart])
			hold on
			drawnow
			ylimits=ylim;
			if length(rpos)==1
			if ~isempty(obj.Erxt)
			vexbext=-obj.Erxt(rpos,zpos)/obj.Bz(zpos,rpos)';
			plot(vexbext*[1 1],ylimits,'k--','displayname',sprintf('V_{ExB,ext}=%1.2g [m/s]',vexbext))
			end

			plot(vexb*[1 1],ylimits,'k-.','displayname',sprintf('V_{ExB,tot}=%1.2g [m/s]',vexb))
			end

			subplot(1,4,3);
			obj.dispV(Vz,Vzend,'Vz [m/s]',[1,timesteppart])

			subplot(1,4,4);
			obj.dispV(sqrt(Vr.^2+(Vthet).^2+Vz.^2),sqrt(Vrend.^2+(Vthetend).^2+Vzend.^2),'Vtot [m/s]',[1,timesteppart],'maxwell')


			sgtitle(sprintf('t=%1.2e[ns] r=[%2.1f, %2.1f] [mm] z=[%2.1f, %2.1f] [mm]',obj.tpart(timestep)1e9, rspan1e3, zspan*1e3))
			obj.savegraph(f,sprintf('%s/%sParts_V_RZ',obj.folder,obj.name),[25 14]);


			end

			end

			function displayEkin(obj,timestep, rpos, zpos)
			%displayEkin plot the kinetic energy distribution function in
			%eV
			%extracted from the markers at position window from rpos(1)
			%rpos(end) and zpos(1) to zpos(end)
			% and at time obj.tpart(timestep)
			%rpos and zpos are given as grid indices
			if(obj.species(1).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(1)) r(rpos(end))];
			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(1)) z(zpos(end))];
			end


			nbp=min(obj.nbparts(timesteppart(1)),obj.species(1).R.nparts);
			R=obj.species(1).R(1:nbp,timesteppart(1),false);
			Z=obj.species(1).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.species(1).R.nparts);
			Rend=obj.species(1).R(1:nbp,timesteppart(2),false);
			Zend=obj.species(1).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.dt1e9, rspan1e3, zspan*1e3))
			obj.savegraph(f,sprintf('%s/%sParts_E_kin',obj.folder,obj.name));


			end

			end

			function displayEkinavg(obj,timestep, rpos, zpos)
			%displayEkin plot the kinetic energy distribution function in
			%eV
			%extracted from the markers at position window from rpos(1)
			%rpos(end) and zpos(1) to zpos(end)
			%and averaged over the time interval timestep
			%rpos and zpos are given as grid indices
			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(1)) r(rpos(end))];
			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(1)) z(zpos(end))];
			end
			tspan=[obj.tpart(timesteppart(1)) obj.tpart(timesteppart(end))];

			Ekintot=[];
			for i=timesteppart(1):timesteppart(end)
			nbp=min(obj.nbparts(i),obj.species(1).R.nparts);
			R=obj.species(1).R(1:nbp,i,false);
			Z=obj.species(1).Z(1:nbp,i,false);
			Ekin=obj.Ekin(1:nbp,i,false);
			ids=R>=rspan(1) & R<=rspan(2) & Z>=zspan(1) & Z<=zspan(2);
			Ekintot=[Ekintot Ekin(ids)'/obj.qe];
			end
			[Counts,edges]=histcounts(Ekintot,'binmethod','sqrt');
			Counts=Counts*obj.weight;
			tot=sum(Counts)
			f=figure('Name',sprintf("%s E_k distrib",obj.file));
			plot([edges(1) 0.5*(edges(2:end)+edges(1:end-1)) edges(end)],[0 Counts 0]);
			sgtitle(sprintf('t=[%2.1f, %2.1f] [ns] r=[%2.1f, %2.1f] [mm] z=[%2.1f, %2.1f] [mm]',tspan1e9, rspan1e3, zspan*1e3))
			end

			function displayVdistribParPer(obj,timestep, rpos, zpos, gcs)
			%displayVdistribParPer plot the velocity distribution function
			% in m/s for the parallel and perpendicular velocity
			%extracted from the markers at position window from rpos(1)
			%rpos(end) and zpos(1) to zpos(end)
			% and at time obj.tpart(timestep)
			%rpos and zpos are given as grid indices
			% gcs define if you give the perpendicular velocity in the
			% guiding center frame or in the laboratory frame
			if(obj.species(1).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(1)) r(rpos(end))];
			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(1)) z(zpos(end))];
			end
			if nargin<5
			gcs=false; % define if we look in the guiding center system
			end

			nbp=min(obj.nbparts(1),obj.species(1).R.nparts);
			R=obj.species(1).R(1:nbp,1,false);
			Z=obj.species(1).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,gcs);
			Vpar=obj.Vpar(1:nbp,1,false);
			Vperp=Vperp(ids);
			Vpar=Vpar(ids);

			nbp=min(obj.nbparts(timesteppart),obj.species(1).R.nparts);
			R=obj.species(1).R(1:nbp,timesteppart,false);
			Z=obj.species(1).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,gcs);
			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)
			if gcs
			lgd='v_\perp gcs [m/s]';
			else
			lgd='v_\perp [m/s]';
			end
			obj.dispV(Vperp,Vperpend,lgd,[1,timesteppart], 'maxwell')

			subplot(1,2,2)
			obj.dispV(abs(Vpar),abs(Vparend),'v_{par} [m/s]',[1,timesteppart],'None')

			sgtitle(sprintf('t=%1.2e[ns] r=[%2.1f, %2.1f] [mm] z=[%2.1f, %2.1f] [mm]',obj.tpart(timestep)1e9, rspan1e3, zspan*1e3))
			obj.savegraph(f,sprintf('%s/%sParts_V_parper',obj.folder,obj.name));
			if gcs
			obj.savegraph(f,sprintf('%s/%sParts_V_parpergcs',obj.folder,obj.name));
			else
			obj.savegraph(f,sprintf('%s/%sParts_V_parper',obj.folder,obj.name));
			end

			end

			end

			function display2DVdistrib(obj,timestep, rpos, zpos, gcs)
			%display2DVdistrib plot the velocity distribution function
			% in m/s for the parallel and perpendicular velocity
			% and for the radial azimuthal velocity
			% as a 2D contour plot the show the velocity phase space distribution
			%extracted from the markers at position window from rpos(1)
			%rpos(end) and zpos(1) to zpos(end)
			% and at time obj.tpart(timestep)
			%rpos and zpos are given as grid indices
			% gcs define if you give the perpendicular velocity in the
			% guiding center frame or in the laboratory frame
			if(obj.species(1).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(1)) r(rpos(end)+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(1)) z(zpos(end)+1)];
			end
			if nargin<5
			gcs=false; % define if we look in the guiding center system
			end

			nbp=min(obj.nbparts(timesteppart),obj.species(1).R.nparts);
			R=obj.species(1).R(1:nbp,timesteppart,false);
			Z=obj.species(1).Z(1:nbp,timesteppart,false);
			ids=R>=rspan(1) & R<=rspan(2) & Z>=zspan(1) & Z<=zspan(2);

			Vperp=obj.Vperp(1:nbp,timesteppart,false,gcs);
			Vpar=obj.Vpar(1:nbp,timesteppart,false);
			Vr=obj.species(1).VR(1:nbp,timesteppart,false);
			Vthet=obj.species(1).VTHET(1:nbp,timesteppart,false);
			Vper=Vperp(ids);
			Vpar=Vpar(ids);

			Vr=Vr(ids);
			Vthet=Vthet(ids);
			nbp=sum(ids(:));

			f=figure('Name',sprintf("%s v parper distrib",obj.file));
			subplot(2,1,1)
			[N,Xedges,Yedges] = histcounts2(Vpar,Vper,20);
			Xedges=(Xedges(1:end-1)+Xedges(2:end))/2;
			Yedges=(Yedges(1:end-1)+Yedges(2:end))/2;
			contourf(Xedges,Yedges,N')
			xlabel('v_{par} [m/s]')
			ylabel('v_{\perp} [m/s]')
			c=colorbar;
			c.Label.String='Counts';

			subplot(2,1,2)
			[N,Xedges,Yedges] = histcounts2(Vthet,Vr,20);
			Xedges=(Xedges(1:end-1)+Xedges(2:end))/2;
			Yedges=(Yedges(1:end-1)+Yedges(2:end))/2;
			contourf(Xedges,Yedges,N')
			%histogram2(Vthet,Vr,'displaystyle','tile','binmethod','auto')
			%scatter(Vthet,Vr)
			xlabel('v_\theta [m/s]')
			ylabel('v_r [m/s]')
			c=colorbar;
			c.Label.String='Counts';

			sgtitle(sprintf('t=%1.2e[ns] r=[%2.1f, %2.1f] [mm] z=[%2.1f, %2.1f] [mm] N=%3i',mean(obj.tpart(timestep))1e9, rspan1e3, zspan*1e3,nbp))
			mkdir(sprintf('%s/vdist',obj.folder))
			if gcs
			obj.savegraph(f,sprintf('%s/vdist/%sParts_V_2dparpergcs_r%iz%it%i',obj.folder,obj.name,floor(mean(rpos)),floor(mean(zpos)),floor(mean(timestep))));
			else
			obj.savegraph(f,sprintf('%s/vdist/%sParts_V_2dparper_r%iz%it%i',obj.folder,obj.name,floor(mean(rpos)),floor(mean(zpos)),floor(mean(timestep))));
			end

			end

			end

			function displaytavgdensity(obj,fieldstart,fieldend)
			%displaytavgdensity plot the fluid density averaged in time for times t2d(fieldstart:fieldend)
			% Displays also the geometry of the metallic boundaries and the
			% magnetic field lines
			f=figure('Name', sprintf('%sfluid_dens',obj.name));
			ax1=gca;
			geomw=obj.geomweight(:,:,1);
			dens=mean(obj.density(:,:,fieldstart:fieldend),3);
			dens(geomw<=0)=0;
			maxdens=max(dens(:));
			[h,curve]=contourf(ax1,obj.zgrid1000,obj.rgrid1000,dens,30);
			hold on
			set(curve,'linecolor','none');
			%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)))]*1000)
			else
			ylim(ax1,[obj.rgrid(1) obj.rgrid(end)]*1000)
			end
			xlim(ax1,[obj.zgrid(1) obj.zgrid(end)]*1000)
			xlabel(ax1,'z [mm]')
			ylabel(ax1,'r [mm]')
			%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= 'electron density [m^{-3}]';
			view(ax1,2)
			colormap(flipud(hot));

			% Draws the magnetic field lines
			Blines=obj.rAthet;
			levels=linspace(min(Blines(obj.geomweight(:,:,1)>0)),max(Blines(obj.geomweight(:,:,1)>0)),20);
			Blines(obj.geomweight(:,:,1)<0)=NaN;
			[~,h1]=contour(ax1,obj.zgrid1000,obj.rgrid1000,Blines,real(levels),'-.','color','k','linewidth',1.5,'Displayname','Magnetic field lines');

			% Draw the metallic boundaries and the geometry itself
			[c1,hContour]=contourf(ax1,obj.zgrid1000,obj.rgrid1000,-geomw,[0,0],'linewidth',1.5);
			drawnow;
			xlim(ax1,[obj.zgrid(1)1000 obj.zgrid(end)1000])
			% Change the color of the metallic boundaries to grey
			hFills=hContour.FacePrims;
			[hFills.ColorType] = deal('truecoloralpha'); % default = 'truecolor'
			try
			hFills(end).ColorData = uint8([150;150;150;255]);
			for idx = 1 : numel(hFills)-1
			hFills(idx).ColorData(4) = 0; % default=255
			end
			catch
			end

			grid on;
			hold on;
			f.PaperOrientation='landscape';
			f.PaperUnits='centimeters';
			papsize=[16 14];
			f.PaperSize=papsize;
			set(ax1,'fontsize',14)
			%axis equal
			obj.savegraph(f,sprintf('%sfluid_dens',obj.name))
			end

			function displaymagfield(obj)
			%displaymagfield display the magnetic field lines and the
			%amplitude of the magnetic field using a contour
			% also show the domain boundaries

			B=obj.B';

			f=figure('Name', sprintf('%s B field',obj.name));
			B(obj.geomweight(:,:,1)<0)=NaN;
			ax1=gca;
			title(sprintf('Configuration'))
			Blv=linspace(min(B(:)),max(B(:)),25);
			if length(Blv)<2
			Blv=min(B(:))*[1 1];
			end
			h=contourf(ax1,obj.zgrid1000,obj.rgrid1000,B,Blv,'Displayname','B [T]', 'linestyle','none');
			hold on;

			%% Magnetic field lines
			Blines=obj.rAthet;
			levels=linspace(min(Blines(obj.geomweight(:,:,1)>0)),max(Blines(obj.geomweight(:,:,1)>0)),30);
			Blines(obj.geomweight(:,:,1)<0)=NaN;
			[~,h1]=contour(ax1,obj.zgrid1000,obj.rgrid1000,Blines,real(levels),'r-.','linewidth',1.5,'Displayname','Magnetic field lines');

			colormap(ax1,'parula')

			% Grey outline
			geomw=obj.geomweight(:,:,1);
			geomw(geomw>0)=NaN;
			geomw(geomw<0)=min(B(:));
			[c1,hContour]=contourf(ax1,obj.zgrid1000,obj.rgrid1000,geomw, [0 0]);

			drawnow;
			caxis([0 0.6])
			xlim(ax1,[obj.zgrid(1)1000 obj.zgrid(end)1000])
			if(obj.conformgeom)
			ylim([ax1 ],[obj.rgrid(1)1000 obj.rgrid(rgridend)1000])
			else
			ylim([ax1],[obj.rgrid(1)1000 obj.rgrid(end)1000])
			end
			%legend([h1],{'Magnetic field lines'},'location','northwest')
			xlabel(ax1,'z [mm]')
			ylabel(ax1,'r [mm]')
			%title(ax1,sprintf('Density t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(fieldstart),M.t2d(fieldend),double(maxdens)))
			c = colorbar(ax1);
			c.Label.String= '\|B\| [T]';
			view(ax1,2)
			%set(h,'edgecolor','none');
			grid on;
			hFills=hContour.FacePrims;
			[hFills.ColorType] = deal('truecoloralpha'); % default = 'truecolor'
			%caxis([min(B(:)) max(B(:))])
			try
			hFills(1).ColorData = uint8([150;150;150;255]);
			for idx = 2 : numel(hFills)
			hFills(idx).ColorData(4) = 0; % default=255
			end
			catch
			end
			[~, name, ~] = fileparts(obj.file);
			if( obj.walltype >=2 && obj.walltype<=4)
			rectangle('Position',[obj.zgrid(1) obj.r_b obj.zgrid(end)-obj.zgrid(1) 0.001]*1e3,'FaceColor',[150 150 150]/255,'Edgecolor','none')
			ylimits=ylim;
			ylim([ylimits(1),ylimits(2)+1])
			end
			if(isempty(obj.spl_bound))
			rectangle('Position',[obj.zgrid(1) obj.r_a-0.001 obj.zgrid(end)-obj.zgrid(1) 0.001]*1e3,'FaceColor',[150 150 150]/255,'Edgecolor','none')
			ylimits=ylim;
			ylim([ylimits(1)-1,ylimits(2)])
			end
			f.PaperUnits='centimeters';
			%axis equal

			set(gca,'fontsize',14)

			papsize=[14 5 ];
			pos=f.Position;
			pos(3)=floor(1.5*pos(3));
			f.Position=pos;
			axis equal

			obj.savegraph(f,sprintf('%s/%s_Bfield',obj.folder,obj.name),papsize);
			end

			function displaySurfFlux(obj,timestep, subdiv, species_id,scalet)
			%displaySurfFlux plot the current densities
			%on the domain boundaries for time t2d(timestep)
			%directly on the boundaries themselves
			%make it easier to see where the currents are collected

			if nargin<3
			subdiv=1;
			end
			if nargin<4
			species_id =1;
			end
			if nargin<5
			scalet = false;
			end
			if species_id ==1
			mflux= obj.Metallicflux(timestep,subdiv);
			lflux= -squeeze(obj.Axialflux(timestep,1,species_id))';
			rflux= squeeze(obj.Axialflux(timestep,length(obj.zgrid),species_id))';
			qe = obj.qsim/obj.weight;
			else
			mflux= obj.MetallicFlux_species(timestep,subdiv);
			lflux= -squeeze(obj.Axialflux(timestep,1,species_id))';
			rflux= squeeze(obj.Axialflux(timestep,length(obj.zgrid),species_id))';
			qe = obj.species(species_id).q/obj.species(species_id).weight*obj.species(1).weight;
			end
			time=obj.t2d(timestep);

			if nargin<3
			ids=1:length(mflux);
			end

			%%
			P=obj.neutcol.neutdensobj.kb300/100;% pressure at room temperature in mbar
			f=figure('name','fluxevol');
			linew=3;
			%obj.displaysplbound(gca,1e3);

			% Grey outline
			geomw=obj.geomweight(:,:,1);
			geomw(geomw>0)=NaN;
			geomw(geomw<0)=0;
			[c1,hContour]=contourf(obj.zgrid1000,obj.rgrid1000,geomw, [0 0]);
			hold on
			contour(obj.zgrid1e3,obj.rgrid1e3,obj.geomweight(:,:,1),[0 0],'b-','linewidth',1.5);
			hold on
			for i=1:length(mflux.p)
			x=mflux.p{i}(1,:)*1000;
			y=mflux.p{i}(2,:)*1000;
			y(end)=NaN;
			c=mean(mflux.gamma{i},2)'*qe/(100^2)/P;
			c(c>0)=0;
			c=abs(c);
			%c=mflux.gamma{i}'*obj.qe/(100^2)/P
			%c=mflux.gamma{i}'*qe/(100^2)/P;
			c(c==0)=NaN;
			patch(x,y,c,'EdgeColor','interp','LineWidth',linew);
			hold on
			end


			x=obj.zgrid(1)ones(size(obj.rgrid))1000;
			y=obj.rgrid*1000;
			y(end)=NaN;
			c=mean(lflux,1)*qe/(100^2)/P;
			%c=lflux*obj.qe/(100^2)/P;
			%c=lflux*qe/(100^2)/P;
			c(c==0)=NaN;
			patch(x,y,c,'EdgeColor','interp','LineWidth',linew);


			x=obj.zgrid(end)ones(size(obj.rgrid))1e3;
			y=obj.rgrid*1000;
			y(end)=NaN;
			c=mean(rflux,1)*qe/(100^2)/P;
			%c=rflux*obj.qe/(100^2)/P;
			%c=rflux*qe/(100^2)/P;
			c(c==0)=NaN;
			patch(x,y,c,'EdgeColor','interp','LineWidth',linew);

			if scalet == false

			title(sprintf('t=%4.2f [ns]',mean(time)*1e9))
			else
			if obj.neutcol.present
			vexb0=(obj.Ez(:,:,1).obj.Br'-obj.Er(:,:,1).obj.Bz')./(obj.B'.^2);
			vexb0(obj.geomweight(:,:,1)<=0)=0;
			E=0.5obj.msim/obj.weightmean(abs(vexb0(:)))^2/obj.qe;
			taucol=1/(obj.neutcol.neutdensmean(abs(vexb0(:)))(obj.sigio(E)+obj.sigmela(E)+obj.sigmio(E)));
			try
			Sio_S=1e17(obj.neutcol.neutdensmean(abs(vexb0(:)))obj.sigio(E))/(obj.maxwellsrce.frequencyobj.weight/(pi(obj.maxwellsrce.rlim(2)^2-obj.maxwellsrce.rlim(1)^2)diff(obj.maxwellsrce.zlim)))
			catch
			end
			tlabel='t/\tau_d [-]';
			else
			taucol=2*pi/obj.omece;
			tlabel='t/\tau_ce [-]';
			end
			title(sprintf('t/\\tau_d=%4.2f [-]', mean(time)/taucol))
			end
			title('')

			c=colorbar;
			c.Label.String= 'j\cdotn [A/(cm^2 mbar)]';
			xlabel('z [mm]')
			ylabel('r [mm]')
			colormap(jet)
			set(gca,'colorscale','log')
			set(gca,'fontsize',10)
			caxis([1e-3 10])
			axis equal
			f.Position=f.Position.*[1 1 1.2 1];
			%% Magnetic field lines
			Blines=obj.rAthet;
			levels=linspace(min(Blines(obj.geomweight(:,:,1)>0)),max(Blines(obj.geomweight(:,:,1)>0)),15);
			Blines(obj.geomweight(:,:,1)<0)=NaN;
			[~,h1]=contour(obj.zgrid1000,obj.rgrid1000,Blines,real(levels),'k-.','linewidth',1.5,'Displayname','Magnetic field lines');
			%axis equal
			drawnow


			% Set boundaries to gray
			hFills=hContour.FacePrims;
			[hFills.ColorType] = deal('truecoloralpha'); % default = 'truecolor'
			try
			hFills(1).ColorData = uint8([150;150;150;255]);
			for idx = 2 : numel(hFills)
			hFills(idx).ColorData(4) = 0; % default=255
			end
			catch
			end

			pos=f.Position;
			%pos(3)=floor(1.5*pos(3));
			f.Position=pos;
			%axis equal

			xlim([obj.zgrid(1) obj.zgrid(end)]1e3+0.1[-1 1])
			%xlim([201e-3 obj.zgrid(end)]1e3+0.1*[-1 1])
			ylim([obj.rgrid(1) obj.rgrid(end)]1e3+0.1[-1 1])
			drawnow
			pause(0.5)
			c.Ticks=[1e-2 1e-1 1 1e1 1e2];
			obj.savegraph(f,sprintf('%s/%s_surfFlux_it2d_%i_spec%i',obj.folder,obj.name,floor(mean(timestep)),species_id),[16 14]);
			end

			% interactive window to display the terms of the pressure tensor
			guiPressure(obj,logdensity,showgrid,fixed,temperature)

			% interactive window to display the electron density, magnetic
			% field lines, electric potential and field at given time steps
			guiFields(obj,logdensity,showgrid,fixed,parper)


			function displaycollfreq(obj)
			%displaycollfreq plot the collision frequencies in Hz/mbar for a range of
			%electron kinetic energies in eV for the different collision
			%processes considered: ionisation and elastic collisions
			E=logspace(1,4,1000);

			v=sqrt(2/obj.msimobj.weightE*obj.qe);
			P=obj.neutcol.neutdensobj.kb300/100;% pressure at room temperature in mbar
			tauio=P./(obj.neutcol.neutdensobj.sigio(E).v);
			tauiom=P./(obj.neutcol.neutdensobj.sigmio(E).v);

			tauelam=P./(obj.neutcol.neutdensobj.sigmela(E).v);

			f=figure('name','t scales coll');
			loglog(E,1./tauio,'displayname','ionisation','linewidth',1.5)
			hold on
			loglog(E,1./tauiom,'displayname','ionisation momentum','linewidth',1.5)
			loglog(E,1./tauelam,'displayname','elastic','linewidth',1.5)
			loglog(E,1./tauio+1./tauiom+1./tauelam,'displayname','total drag','linewidth',1.5)

			loglog([E(1) E(end)],1./(2pi/obj.omece) [1 1],'--','displayname','cyclotronic','linewidth',1.5)

			xlabel('Electron kinetic energy [eV]')
			ylabel('\nu [Hz/mbar]')
			legend('location','southeast')
			grid on

			obj.savegraph(f,sprintf('%s/collfreqscales',obj.folder),[14 12]);
			end

			function displaycrosssec(obj)
			%displaycrosssec plot the collision crosssections in m^2 for a range of
			%electron kinetic energies in eV for the different collision
			%processes considered: ionisation and elastic collisions
			E=logspace(1,4,1000);


			sig_io=obj.sigio(E);
			sig_iom=obj.sigmio(E);

			sig_elam=obj.sigmela(E);

			f=figure('name','t scales coll');
			loglog(E,sig_io,'displayname','ionisation','linewidth',1.5)
			hold on
			loglog(E,sig_iom,'displayname','ionisation momentum','linewidth',1.5)
			loglog(E,sig_elam,'displayname','elastic','linewidth',1.5)
			loglog(E,sig_io+sig_elam+sig_iom,'displayname','total drag','linewidth',1.5)


			xlabel('Energy [eV]')
			ylabel('\sigma [m^{2}]')
			legend('location','southeast')
			grid on

			obj.savegraph(f,sprintf('%s/coll_cross_sec_scales',obj.folder),[14 12]);
			end

			%------------------------------------------
			% Helper functions needed for other functions

			function [zpos,rpos]=getpos(obj,tstep)
			% interactive window to return an specific axial and radial
			% position picked from the cloud density
			if nargin<2
			tstep=length(obj.t2d);
			end
			n=obj.density(:,:,tstep);
			n(obj.geomweight(:,:,1)<0)=NaN;

			figure
			contourf(obj.zgrid,obj.rgrid,n);
			xlabel('z [m]')
			ylabel('r [m]')
			[x,y]=ginput(1);
			zpos=find(x>obj.zgrid,1,'last');
			rpos=find(y>obj.rgrid,1,'last');

			hold on
			plot(obj.zgrid(zpos),obj.rgrid(rpos),'rx')

			fprintf('zpos=%i rpos=%i z=%1.4f r=%1.4f\n',zpos,rpos,obj.zgrid(zpos),obj.rgrid(rpos))

			end

			function changed=ischanged(obj)
			%ischanged Check if the file has been changed since the initial loading of the file
			%and if 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, vd)
			%dispV generic functio to plot the velocity distribution and
			%comapare two timesteps V and Vend at time t(1) and t(2)
			if nargin<6
			dist='gaussian';
			end
			if nargin<7
			vd=0;
			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','sqrt');
			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','sqrt');
			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));

			if strcmp(dist,'maxwell')
			vfit=linspace(0,edges(end),300);
			a=vmeanend/sqrt(2);
			dist=sqrt(2/pi)vfit.^2.exp(-((vfit).^2-vd^2)/2/a^2)/a^3;
			dist=dist/max(dist);
			plot(vfit,max(Counts)*dist,'displayname',sprintf('Maxw mu=%2.2g sigma=%2.2g',vmeanend,vthermend))
			elseif strcmp(dist,'gaussian')
			vfit=linspace(edges(1),edges(end),300);
			dist=exp(-(vfit-vmeanend).^2/2/vthermend^2);
			plot(vfit,max(Counts)*dist,'displayname',sprintf('gauss mu=%2.2g sigma=%2.2g',vmeanend,vthermend))
			end
			ylabel('counts')
			xlabel(label)
			grid on
			legend('location','southoutside','orientation','vertical')
			end

			function cross_sec=fit_cross_sec(obj,energy,crosssec_table)
			%Interpolate the cross-section at the given energy using the
			%crosssec_table and an exponential fitting
			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=min(size(crosssec_table,1)-1,id);
			id=max(1,id);

			cross_sec=crosssec_table(id,2)*(energy/crosssec_table(id,1))^crosssec_table(id,3);
			%cross_sec=crosssec_table(id,2)+(crosssec_table(id+1,2)-crosssec_table(id,2))*(energy-crosssec_table(id,1))/((crosssec_table(id+1,1)-crosssec_table(id,1)));
			end

			function fighandle=savegraph(obj, fighandle, name, papsize)
			%% Saves the given figure as a pdf a .fig and an eps using export_fig
			% fighandle is the handle to the figure being saved
			% name is the path to the saved file without extension
			% papsize is the dimension of the pdf print in cm
			if (nargin < 4)
			papsize=[14 16];
			end
			%export_fig(fighandle,name,'-png','-r300')
			%exportgraphics(fighandle,sprintf('%s.png',name),'Resolution',300)
			%print(fighandle,name,'-dpdf','-fillpage')
			%savefig(fighandle,name)
			%fighandle.PaperUnits='centimeters';
			%set(fighandle, 'Color', 'w');
			%fighandle.PaperSize=papsize;
			%exportgraphics(fighandle,sprintf('%s.eps',name))
			%export_fig(fighandle,name,'-eps','-painters')

			end

			function sig=dsigmaio(obj,Ekin, Ebar, Ei, E0, chi, gamma)
			% calculates the integrand used for the ionisation collision
			% cross section for momentum exchange for the incoming electron
			% it is only used by obj.sigmiopre
			gamma=reshape(gamma,1,[],1);
			chi=reshape(chi,1,1,[]);

			siggamma=sin(gamma).(E0^2+8(1-chi)(Ekin-Ei)E0)./(E0+4(1-chi)(Ekin-Ei)-4(1-chi)(Ekin-Ei).*cos(gamma)).^2/2;

			sigchi=(Ekin-Ei)./(Ebaratan((Ekin-Ei)/(2Ebar)).(1+(chi(Ekin-Ei)/Ebar).^2));

			dp=1- trapz(gamma,sqrt((1-chi).(1-Ei/Ekin)).cos(gamma).siggamma,2);%- trapz(gamma,sqrt((1-chi).(1-Ei/Ekin)).cos(gamma).siggamma,2);
			sig=sigchi.*dp;
			end

			function sigm=sigmiopre(obj,E, init)
			% returns the precalculated values used for the interpolation
			% of the ionisation collision cross-section for momentum
			% exchange for the incoming electron
			if nargin <3
			init=false;
			end
			if(~init &&( ~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=numel(E);

			nchi=300;
			ngamma=300;
			gamma=linspace(0,pi,ngamma);
			chi=linspace(0,0.5,nchi);
			%sigm2=zeros(nE,nchi);
			sigm=zeros(size(E));

			for i=1:nE
			if(E(i)>=Ei)
			sigm2=zeros(nchi,1);
			for j=1:nchi
			%sigm2(j)=trapz(alpha,trapz(gamma,obj.dsigmaio(E(i),Ebar,Ei,E0,chi(j),alpha,gamma),2),1);
			sigm2(j)=obj.dsigmaio(E(i),Ebar,Ei,E0,chi(j),gamma);
			end
			sigm(i)=trapz(chi,sigm2)*obj.sigio(E(i),init);
			%sigm(i)=trapz(chi,trapz(alpha,trapz(gamma,dsigmaio(obj,E(i),Ebar,Ei,E0,chi,alpha,gamma),2),1),3)*obj.sigio(E(i),init);
			end
			end
			end

			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

			function displaygeomweight(obj)
			%% Display the geometric weight
			f=figure;
			ax1=gca;

			geomw=obj.geomweight(:,:,1);
			[h,curve]=contourf(ax1,obj.zgrid1000,obj.rgrid1000,geomw);
			hold on
			set(curve,'linecolor','none');
			xlabel(ax1,'z [mm]')
			ylabel(ax1,'r [mm]')
			c = colorbar(ax1);
			c.Label.String= 'geometric weight';
			caxis([0 1]);
			view(ax1,2)

			% Draw the metallic boundaries and the geometry itself
			[c1,hContour]=contourf(ax1,obj.zgrid1000,obj.rgrid1000,-geomw,[0,0],'linewidth',1.5);
			drawnow;
			% Change the color of the metallic boundaries to grey
			hFills=hContour.FacePrims;
			[hFills.ColorType] = deal('truecoloralpha');
			hFills(end).ColorData = uint8([150;150;150;255]);
			for idx = 1 : numel(hFills)-1
			hFills(idx).ColorData(4) = 0;
			end

			grid on
			f.PaperOrientation='landscape';
			f.PaperUnits='centimeters';
			papsize=[16 14];
			f.PaperSize=papsize;
			set(ax1,'fontsize',14);
			pos=f.Position;
			pos(3)=1.5*pos(3);
			f.Position=pos;
			end

			function displaydensityprofile(obj,zpos,fieldstep)
			%% Display the density in 1D at the axial position zpos and time obj.t2d(timestep)
			f=figure;
			geomw=obj.geomweight(:,zpos,1);
			r=obj.rgrid*1e3;
			dens=obj.N(:,zpos,fieldstep);
			dens(geomw<=0)=0;
			densmax=max(dens)
			ax1=gca;
			p=plot(ax1,r,dens,'linewidth',1.5);
			xlim(ax1,[obj.rgrid(1) obj.rgrid(end)]*1e3)
			xlabel(ax1,'r [mm]')
			ylabel(ax1,'n_e[m^{-3}]');
			title(['t = ',num2str(obj.t2d(fieldstep)),' s']);
			end

			function displaymaxdensity(obj,timesteps)
			%% Displays the time evolution of the maximum density in the electron cloud
			t=obj.t2d(timesteps)*1e9;
			geomw=obj.geomweight(:,:,1);
			maxdens=zeros(length(t),1);
			for i=1:length(timesteps)
			dens=obj.N(:,:,timesteps(i));
			dens(geomw<=0)=0;
			maxdens(i)=max(dens(:));
			end
			f=figure;
			ax=axes(f);
			plot(t,maxdens,'linewidth',1.8);
			grid on
			xlabel('t [ns]')
			ylabel('n_{e,max} [m^{-3}]')
			set(ax,'fontsize',12)
			end

			function displayavgdensity(obj,timesteps)
			%% Displays the time evolution of the average density in the electron cloud
			t=obj.t2d(timesteps)*1e9;
			geomw=obj.geomweight(:,:,1);
			avgdens=zeros(length(t),1);
			for i=1:length(timesteps)
			dens=obj.N(:,:,timesteps(i));
			dens(geomw<=0)=0;
			maxdens=max(dens(:));
			threshold=maxdens/100;
			dens(dens<threshold)=0;
			a=(dens>0);
			avgdens(i)=sum(dens(:))/sum(a(:));
			end
			f=figure;
			ax=axes(f);
			plot(t,avgdens,'linewidth',1.8);
			grid on
			xlabel('t [ns]')
			ylabel('n_{e,avg} [m^{-3}]')
			set(ax,'fontsize',12)
			end

			function displaychargecloud(obj)
			%% Displays the time evolution of the total charge in the electron cloud
			t=obj.t2d(1:10:end);
			rlim=59:120;
			zlim=125:275;
			r=obj.rgrid(rlim);
			z=obj.zgrid(zlim);
			geomw=obj.geomweight(rlim,zlim,1);
			dz=z(2)-z(1);
			dV=2pidzr.[diff(r)' 0]';
			Qcloud=zeros(length(t),1);
			for i=1:length(t)
			dens=obj.density(rlim,zlim,(i-1)*10+1);
			dens(geomw<=0)=0;
			Nelec=sum(sum(diag(dV)*dens));
			Qcloud(i)=obj.qe*Nelec;
			end
			f=figure;
			plot(t,Qcloud);
			xlabel('t [s]')
			ylabel('Q_{cloud} [C]')
			end

			function displaycurrentstrex(obj,timesteps,species_id)
			%% Display the time evolution of the outgoing currents
			if nargin <3
			species_id=1;
			end
			nmean=1;
			subdiv=1;
			if species_id ==1
			[currents,pos]=obj.OutCurrents(timesteps,subdiv);
			P=obj.neutcol.neutdensobj.kb300/100;% pressure at room temperature in mbar
			currents=currents/P;
			else % for ionic currents
			[currents,pos]=obj.OutCurrents_species(timesteps,subdiv);
			currents=currents/obj.species(species_id).weight*obj.species(1).weight;
			P=obj.neutcol.neutdensobj.kb300/100;% pressure at room temperature in mbar
			currents=currents/P;
			end
			% plot
			f=figure;
			p(1)=plot(obj.t2d(timesteps)*1e9,movmean(currents(4,:),nmean),'Displayname','top flange','linewidth',1.8);
			hold on
			p(2)=plot(obj.t2d(timesteps)*1e9,movmean(currents(5,:),nmean),'Displayname','inner electrode','linewidth',1.8);
			hold on
			p(3)=plot(obj.t2d(timesteps)*1e9,movmean(currents(6,:),nmean),'Displayname','Faraday','linewidth',1.8);
			hold on
			p(4)=plot(obj.t2d(timesteps)*1e9,movmean(currents(7,:),nmean),'Displayname','outer electrode','linewidth',1.8);
			legend;
			xlabel('t [ns]')
			ylabel('I/p_n [A/mbar]')
			grid on
			set(gca,'fontsize',12)
			fprintf('mean O.E. current: %f [A/mbar]\n',mean(currents(7,:),'omitnan'));
			fprintf('mean T.F. current: %f [A/mbar]\n',mean(currents(4,:),'omitnan'));
			end

			function displaypartslosses(obj,tsteps,nmean)
			%% Display the time evolution of the number of lost particles at each boundary
			if nargin <2
			tsteps=1:length(obj.t0d);
			end
			if nargin <3
			nmean=1;
			end
			t=obj.t0d(tsteps)*1e9;
			P=obj.neutcol.neutdensobj.kb300/100; % pressure at room temperature in mbar
			%P=1;
			deltat=obj.dt(1)*double(obj.it0);
			partsloss=obj.losses(:,tsteps)abs(obj.qsim)/(Pdeltat);

			f=figure;
			ax1=gca;
			map=colormap(lines);
			set(ax1,'linestyleorder',{'-',':','--','*','+'},'ColorOrder',map(2:7,:), 'NextPlot','replacechildren')
			p(1)=plot(ax1,t,movmean(partsloss(5,:),nmean),'Displayname','top flange','linewidth',1.8);
			hold on
			p(2)=plot(ax1,t,movmean(partsloss(6,:),nmean),'Displayname','inner electrode','linewidth',1.8);
			hold on
			p(3)=plot(ax1,t,movmean(partsloss(7,:),nmean),'Displayname','Faraday probe','linewidth',1.8);
			hold on
			p(4)=plot(ax1,t,movmean(partsloss(9,:),nmean),'Displayname','outer electrode','linewidth',1.8);
			legend;
			xlabel('t [ns]')
			ylabel('I/p_n [A/mbar]')
			grid on
			set(ax1,'fontsize',12)
			fprintf('mean O.E. current: %f [A/mbar]\n',mean(partsloss(9,:),'omitnan'));
			fprintf('mean T.F. current: %f [A/mbar]\n',mean(partsloss(5,:),'omitnan'));
			end

			function displaydensityrz(obj,thetindex,tstep,zlims,rlims)
			%% Plot the density for time t2d(tstep) in the r-z plane
			if nargin<2
			thetindex=1;
			end
			if nargin<3
			tstep=1;
			end
			if nargin<4
			zlims=[obj.zgrid(1),obj.zgrid(end)]*1e3;
			end
			if nargin<5
			rlims=[obj.rgrid(1),obj.rgrid(end)]*1e3;
			end

			geomw=obj.geomweight(:,:,1);
			dens=squeeze(obj.N(:,thetindex,:,tstep));
			dens(geomw<=0)=0;
			maxdens=max(dens(:));

			f=figure('Name', 'Density','Color',[1 1 1],'Position', [200, 200, 700, 400]);
			ax1=gca;
			[h,curve]=contourf(ax1,obj.zgrid1000,obj.rgrid1000,dens,30);
			hold on
			set(curve,'linecolor','none');
			xlabel(ax1,'z [mm]')
			ylabel(ax1,'r [mm]')
			set(ax1,'fontsize',14)
			c = colorbar(ax1);
			c.Label.String= 'n_e [m^{-3}]';
			c.FontSize=14;
			colormap(flipud(hot));
			caxis([0 maxdens]);

			% Draw the metallic boundaries and the geometry itself
			[c1,hContour]=contourf(ax1,obj.zgrid1000,obj.rgrid1000,-geomw,[0,0],'linewidth',1.5);
			drawnow;
			xlim(ax1,zlims)
			ylim(ax1,rlims)
			% Change the color of the metallic boundaries to grey
			hFills=hContour.FacePrims;
			[hFills.ColorType] = deal('truecoloralpha');
			hFills(end).ColorData = uint8([150;150;150;255]);
			for idx = 1 : numel(hFills)-1
			hFills(idx).ColorData(4) = 0;
			end

			% % Magnetic field lines
			% Blines=obj.rAthet;
			% levels=linspace(min(Blines(obj.geomweight(:,:,1)>0)),max(Blines(obj.geomweight(:,:,1)>0)),35);
			% [~,h1]=contour(ax1,obj.zgrid1000,obj.rgrid1000,Blines,real(levels),'-.','color','k','linewidth',1.5,'Displayname','Magnetic field lines');

			% Equipotential lines
			Pot=squeeze(obj.pot(:,thetindex,:,tstep));
			Pot(obj.geomweight(:,:,1)<0)=NaN;
			potcolor='b';
			[c1,h2]=contour(ax1,obj.zgrid1000,obj.rgrid1000,Pot,'--','color',potcolor,'ShowText','on','linewidth',1.5,'Displayname','Equipotentials [kV]');
			clabel(c1,h2,'Color',potcolor)
			clabel(c1,h2, 'labelspacing', 100);
			clabel(c1,h2, 'FontSize', 11);
			%legend([h1,h2],'Magnetic field lines','Equipotentials [V]'},'location','northeast')
			legend(h2,'Equipotentials [V]','location','northeast')

			title(['t = ',num2str(obj.t2d(tstep)),' s']);
			end

			function displaydensityrthet(obj,zindex,tstep,perturb)
			%% Plot the density for time t2d(tstep) in the r-theta plane
			if nargin<2
			zindex=obj.nz/2;
			end
			if nargin<3
			tstep=1;
			end
			if nargin<4
			perturb=false;
			end
			if nargin<5
			interp=false;
			end

			t=obj.t2d(tstep);
			r=obj.rgrid;
			thet=[obj.thetgrid' 2*pi];

			dens=[squeeze(obj.N(:,:,zindex,tstep)) squeeze(obj.N(:,1,zindex,tstep))];

			if perturb
			dens0=[squeeze(obj.N(:,:,zindex,1)) squeeze(obj.N(:,1,zindex,1))];
			dens=dens-dens0;
			end

			geomw=zeros(size(dens));
			for i=1:size(dens,2)
			geomw(:,i)=squeeze(obj.geomweight(:,zindex,1));
			end
			dens(geomw<=0)=NaN;

			maxdens=max(dens(:));

			[R,THET]=meshgrid(r,thet);
			X=R.*cos(THET);
			Y=R.*sin(THET);

			f=figure('Name', 'Density','Color',[1 1 1],'Position', [200, 200, 700, 600]);
			ax1=gca;
			[h,curve]=contourf(ax1,X1000,Y1000,dens',30);
			hold on
			set(curve,'linecolor','none');
			xlabel(ax1,'x [mm]')
			ylabel(ax1,'y [mm]')
			set(ax1,'fontsize',14);
			c = colorbar(ax1);
			c.Label.String= 'n_e [m^{-3}]';
			c.FontSize=14;
			colormap(flipud(hot));
			caxis([0 maxdens]);

			% Draw the metallic boundaries and the geometry itself
			[c1,hContour]=contourf(ax1,X1000,Y1000,-geomw',[0,0],'linewidth',1.5);
			drawnow;
			% Change the color of the metallic boundaries to grey
			hFills=hContour.FacePrims;
			[hFills.ColorType] = deal('truecoloralpha');
			hFills(end).ColorData = uint8([150;150;150;255]);
			for idx = 1 : numel(hFills)-1
			hFills(idx).ColorData(4) = 0;
			end

			% Equipotential lines
			% Pot=[squeeze(obj.pot(:,:,zindex,tstep)) squeeze(obj.pot(:,1,zindex,tstep))];
			% Pot(geomw<=0)=NaN;
			% potcolor='b';
			% [c1,h2]=contour(ax1,X1000,Y1000,Pot','--','color',potcolor,'ShowText','on','linewidth',1.5,'Displayname','Equipotentials [kV]');
			% clabel(c1,h2,'Color',potcolor)
			% clabel(c1,h2, 'labelspacing', 100);
			% clabel(c1,h2, 'FontSize', 11);
			% legend(h2,'Equipotentials [V]','location','northeast')

			axis equal
			title(['t = ',num2str(obj.t2d(tstep)),' s']);
			end

			function displaydensityrthet2(obj,zindex,tstep)
			%% Plot the density for time t2d(tstep) in the r-theta plane
			if nargin<2
			zindex=obj.nz/2;
			end
			if nargin<3
			tstep=1;
			end

			t=obj.t2d(tstep);
			r=obj.rgrid;
			thet=linspace(0,2*pi,1000);
			[R,THET]=meshgrid(r,thet);
			X=R.*cos(THET);
			Y=R.*sin(THET);

			densmode=squeeze(obj.rhsreal(:,zindex,:,tstep)+ 1i*obj.rhsimag(:,zindex,:,tstep));
			dens=zeros(length(r),length(thet));
			for k=1:length(thet)
			for m=1:obj.nthet/2+1
			arg=1idouble(m-1)thet(k);
			dens(:,k)=dens(:,k)+densmode(:,m)*exp(arg)/(double(obj.nthet)/2+1);
			end
			end
			dens=-dens4piobj.eps_0obj.phinorm*obj.rnorm/obj.qe;
			maxdens=max(dens(:));
			geomw=zeros(size(dens));
			for i=1:size(dens,2)
			geomw(:,i)=squeeze(obj.geomweight(:,zindex,1));
			end
			dens(geomw<=0)=NaN;

			f=figure('Name', 'Density','Color',[1 1 1],'Position', [200, 200, 700, 600]);
			ax1=gca;
			[h,curve]=contourf(ax1,X1000,Y1000,dens',30);
			hold on
			set(curve,'linecolor','none');
			xlabel(ax1,'x [mm]')
			ylabel(ax1,'y [mm]')
			set(ax1,'fontsize',14);
			c = colorbar(ax1);
			c.Label.String= 'n_e [m^{-3}]';
			c.FontSize=14;
			colormap(flipud(hot));
			caxis([0 maxdens]);

			% Draw the metallic boundaries and the geometry itself
			[c1,hContour]=contourf(ax1,X1000,Y1000,-geomw',[0,0],'linewidth',1.5);
			drawnow;
			% Change the color of the metallic boundaries to grey
			hFills=hContour.FacePrims;
			[hFills.ColorType] = deal('truecoloralpha');
			hFills(end).ColorData = uint8([150;150;150;255]);
			for idx = 1 : numel(hFills)-1
			hFills(idx).ColorData(4) = 0;
			end
			axis equal
			title(['t = ',num2str(obj.t2d(tstep)),' s']);
			end

			function displayextfield(obj,zlims,rlims)
			%% Plot the external electric field for time t2d(tstep)
			if nargin<2
			zlims=[obj.zgrid(1),obj.zgrid(end)]*1e3;
			end
			if nargin<3
			rlims=[obj.rgrid(1),obj.rgrid(end)]*1e3;
			end

			geomw=obj.geomweight(:,:,1);
			disppot=obj.potxt;
			dispEr=obj.Erxt;
			dispEthet=zeros(size(disppot));
			dispEz=obj.Ezxt;
			disppot(geomw<=0)=0;
			dispEr(geomw<=0)=0;
			dispEz(geomw<=0)=0;

			f=figure('Name', 'External field','Color',[1 1 1],'Position', [100, 100, 1200, 600]);

			%% Electrostatic potential
			subplot(2,2,1)
			ax0=gca;

			[h,curve]=contourf(ax0,obj.zgrid1000,obj.rgrid1000,disppot,30);
			hold on
			set(curve,'linecolor','none');
			xlabel(ax0,'z [mm]')
			ylabel(ax0,'r [mm]')
			set(ax0,'fontsize',14);
			c = colorbar(ax0);
			c.Label.String= '\phi_{ext} [V]';
			c.FontSize=14;

			% Draw the metallic boundaries and the geometry itself
			[c1,hContour]=contourf(ax0,obj.zgrid1000,obj.rgrid1000,-geomw,[0,0],'linewidth',1.5);
			drawnow;
			xlim(ax0,zlims)
			ylim(ax0,rlims)
			% Change the color of the metallic boundaries to grey
			hFills=hContour.FacePrims;
			[hFills.ColorType] = deal('truecoloralpha');
			hFills(end).ColorData = uint8([150;150;150;255]);
			for idx = 1 : numel(hFills)-1
			hFills(idx).ColorData(4) = 0;
			end

			%% Radial electric field
			subplot(2,2,2)
			ax1=gca;

			[h,curve]=contourf(ax1,obj.zgrid1000,obj.rgrid1000,dispEr,30);
			hold on
			set(curve,'linecolor','none');
			xlabel(ax1,'z [mm]')
			ylabel(ax1,'r [mm]')
			set(ax1,'fontsize',14);
			c = colorbar(ax1);
			c.Label.String= 'E_{r,ext} [V/m]';
			c.FontSize=14;

			% Draw the metallic boundaries and the geometry itself
			[c1,hContour]=contourf(ax1,obj.zgrid1000,obj.rgrid1000,-geomw,[0,0],'linewidth',1.5);
			drawnow;
			xlim(ax1,zlims)
			ylim(ax1,rlims)
			% Change the color of the metallic boundaries to grey
			hFills=hContour.FacePrims;
			[hFills.ColorType] = deal('truecoloralpha');
			hFills(end).ColorData = uint8([150;150;150;255]);
			for idx = 1 : numel(hFills)-1
			hFills(idx).ColorData(4) = 0;
			end

			%% Azimuthal electric field
			subplot(2,2,3)
			ax2=gca;

			[h,curve]=contourf(ax2,obj.zgrid1000,obj.rgrid1000,dispEthet,30);
			hold on
			set(curve,'linecolor','none');
			xlabel(ax2,'z [mm]')
			ylabel(ax2,'r [mm]')
			set(ax2,'fontsize',14);
			c = colorbar(ax2);
			c.Label.String= 'E_{\theta,ext} [V/m]';
			c.FontSize=14;

			% Draw the metallic boundaries and the geometry itself
			[c1,hContour]=contourf(ax2,obj.zgrid1000,obj.rgrid1000,-geomw,[0,0],'linewidth',1.5);
			drawnow;
			xlim(ax2,zlims)
			ylim(ax2,rlims)
			% Change the color of the metallic boundaries to grey
			hFills=hContour.FacePrims;
			[hFills.ColorType] = deal('truecoloralpha');
			hFills(end).ColorData = uint8([150;150;150;255]);
			for idx = 1 : numel(hFills)-1
			hFills(idx).ColorData(4) = 0;
			end

			%% Axial electric field
			subplot(2,2,4)
			ax3=gca;

			[h,curve]=contourf(ax3,obj.zgrid1000,obj.rgrid1000,dispEz,30);
			hold on
			set(curve,'linecolor','none');
			xlabel(ax3,'z [mm]')
			ylabel(ax3,'r [mm]')
			set(ax3,'fontsize',14);
			c = colorbar(ax3);
			c.Label.String= 'E_{z,ext} [V/m]';
			c.FontSize=14;

			% Draw the metallic boundaries and the geometry itself
			[c1,hContour]=contourf(ax3,obj.zgrid1000,obj.rgrid1000,-geomw,[0,0],'linewidth',1.5);
			drawnow;
			xlim(ax3,zlims)
			ylim(ax3,rlims)
			% Change the color of the metallic boundaries to grey
			hFills=hContour.FacePrims;
			[hFills.ColorType] = deal('truecoloralpha');
			hFills(end).ColorData = uint8([150;150;150;255]);
			for idx = 1 : numel(hFills)-1
			hFills(idx).ColorData(4) = 0;
			end
			end

			function displayfieldrz(obj,thetindex,tstep,zlims,rlims)
			%% Plot the electric field for time t2d(tstep) in the r-z plane
			if nargin<2
			thetindex=1;
			end
			if nargin<3
			tstep=1;
			end
			if nargin<4
			zlims=[obj.zgrid(1),obj.zgrid(end)]*1e3;
			end
			if nargin<5
			rlims=[obj.rgrid(1),obj.rgrid(end)]*1e3;
			end

			geomw=obj.geomweight(:,:,1);
			disppot=squeeze(obj.pot(:,thetindex,:,tstep));
			dispEr=squeeze(obj.Er(:,thetindex,:,tstep));
			dispEthet=squeeze(obj.Ethet(:,thetindex,:,tstep));
			dispEz=squeeze(obj.Ez(:,thetindex,:,tstep));
			disppot(geomw<=0)=0;
			dispEr(geomw<=0)=0;
			dispEthet(geomw<=0)=0;
			dispEz(geomw<=0)=0;

			f=figure('Name', 'Self-consistent field','Color',[1 1 1],'Position', [100, 100, 1200, 600]);
			suptitle(['t = ',num2str(obj.t2d(tstep)),' s']);

			%% Electrostatic potential
			subplot(2,2,1)
			ax0=gca;

			[h,curve]=contourf(ax0,obj.zgrid1000,obj.rgrid1000,disppot,30);
			hold on
			set(curve,'linecolor','none');
			xlabel(ax0,'z [mm]')
			ylabel(ax0,'r [mm]')
			set(ax0,'fontsize',14);
			c = colorbar(ax0);
			c.Label.String= '\phi [V]';
			c.FontSize=14;

			% Draw the metallic boundaries and the geometry itself
			[c1,hContour]=contourf(ax0,obj.zgrid1000,obj.rgrid1000,-geomw,[0,0],'linewidth',1.5);
			drawnow;
			xlim(ax0,zlims)
			ylim(ax0,rlims)
			% Change the color of the metallic boundaries to grey
			hFills=hContour.FacePrims;
			[hFills.ColorType] = deal('truecoloralpha');
			hFills(end).ColorData = uint8([150;150;150;255]);
			for idx = 1 : numel(hFills)-1
			hFills(idx).ColorData(4) = 0;
			end

			%% Radial electric field
			subplot(2,2,2)
			ax1=gca;

			[h,curve]=contourf(ax1,obj.zgrid1000,obj.rgrid1000,dispEr,30);
			hold on
			set(curve,'linecolor','none');
			xlabel(ax1,'z [mm]')
			ylabel(ax1,'r [mm]')
			set(ax1,'fontsize',14);
			c = colorbar(ax1);
			c.Label.String= 'E_r [V/m]';
			c.FontSize=14;

			% Draw the metallic boundaries and the geometry itself
			[c1,hContour]=contourf(ax1,obj.zgrid1000,obj.rgrid1000,-geomw,[0,0],'linewidth',1.5);
			drawnow;
			xlim(ax1,zlims)
			ylim(ax1,rlims)
			% Change the color of the metallic boundaries to grey
			hFills=hContour.FacePrims;
			[hFills.ColorType] = deal('truecoloralpha');
			hFills(end).ColorData = uint8([150;150;150;255]);
			for idx = 1 : numel(hFills)-1
			hFills(idx).ColorData(4) = 0;
			end

			%% Azimuthal electric field
			subplot(2,2,3)
			ax2=gca;

			[h,curve]=contourf(ax2,obj.zgrid1000,obj.rgrid1000,dispEthet,30);
			hold on
			set(curve,'linecolor','none');
			xlabel(ax2,'z [mm]')
			ylabel(ax2,'r [mm]')
			set(ax2,'fontsize',14);
			c = colorbar(ax2);
			c.Label.String= 'E_\theta [V/m]';
			c.FontSize=14;

			% Draw the metallic boundaries and the geometry itself
			[c1,hContour]=contourf(ax2,obj.zgrid1000,obj.rgrid1000,-geomw,[0,0],'linewidth',1.5);
			drawnow;
			xlim(ax2,zlims)
			ylim(ax2,rlims)
			% Change the color of the metallic boundaries to grey
			hFills=hContour.FacePrims;
			[hFills.ColorType] = deal('truecoloralpha');
			hFills(end).ColorData = uint8([150;150;150;255]);
			for idx = 1 : numel(hFills)-1
			hFills(idx).ColorData(4) = 0;
			end

			%% Axial electric field
			subplot(2,2,4)
			ax3=gca;

			[h,curve]=contourf(ax3,obj.zgrid1000,obj.rgrid1000,dispEz,30);
			hold on
			set(curve,'linecolor','none');
			xlabel(ax3,'z [mm]')
			ylabel(ax3,'r [mm]')
			set(ax3,'fontsize',14);
			c = colorbar(ax3);
			c.Label.String= 'E_z [V/m]';
			c.FontSize=14;

			% Draw the metallic boundaries and the geometry itself
			[c1,hContour]=contourf(ax3,obj.zgrid1000,obj.rgrid1000,-geomw,[0,0],'linewidth',1.5);
			drawnow;
			xlim(ax3,zlims)
			ylim(ax3,rlims)
			% Change the color of the metallic boundaries to grey
			hFills=hContour.FacePrims;
			[hFills.ColorType] = deal('truecoloralpha');
			hFills(end).ColorData = uint8([150;150;150;255]);
			for idx = 1 : numel(hFills)-1
			hFills(idx).ColorData(4) = 0;
			end
			end

			function displayfieldrthet(obj,zindex,tstep,perturb)
			%% Plot the electric field for time t2d(tstep) in the r-theta plane
			if nargin<2
			zindex=obj.nz/2;
			end
			if nargin<3
			tstep=1;
			end
			if nargin<4
			perturb=false;
			end

			r=obj.rgrid;
			thet=[obj.thetgrid' 2*pi];
			[R,THET]=meshgrid(r,thet);
			X=R.*cos(THET);
			Y=R.*sin(THET);

			disppot=[squeeze(obj.pot(:,:,zindex,tstep)) squeeze(obj.pot(:,1,zindex,tstep))]';
			dispEr=[squeeze(obj.Er(:,:,zindex,tstep)) squeeze(obj.Er(:,1,zindex,tstep))]';
			dispEthet=[squeeze(obj.Ethet(:,:,zindex,tstep)) squeeze(obj.Ethet(:,1,zindex,tstep))]';
			dispEz=[squeeze(obj.Ez(:,:,zindex,tstep)) squeeze(obj.Ez(:,1,zindex,tstep))]';

			if perturb
			pot0=[squeeze(obj.pot(:,:,zindex,1)) squeeze(obj.pot(:,1,zindex,1))]';
			Er0=[squeeze(obj.Er(:,:,zindex,1)) squeeze(obj.Er(:,1,zindex,1))]';
			Ethet0=[squeeze(obj.Ethet(:,:,zindex,1)) squeeze(obj.Ethet(:,1,zindex,1))]';
			Ez0=[squeeze(obj.Ez(:,:,zindex,1)) squeeze(obj.Ez(:,1,zindex,1))]';
			disppot=disppot-pot0;
			dispEr=dispEr-Er0;
			dispEthet=dispEthet-Ethet0;
			dispEz=dispEz-Ez0;
			end

			geomw=zeros(length(thet),length(r));
			for k=1:length(thet)
			geomw(k,:)=squeeze(obj.geomweight(:,zindex,1));
			end
			disppot(geomw<=0)=NaN;
			dispEr(geomw<=0)=NaN;
			dispEthet(geomw<=0)=NaN;
			dispEz(geomw<=0)=NaN;

			f=figure('Name', 'Self-consistent field','Color',[1 1 1],'Position', [100, 100, 1100, 800]);
			suptitle(['t = ',num2str(obj.t2d(tstep)),' s'])

			%% Electrostatic potential
			subplot(2,2,1);
			ax1=gca;
			[h1,curve1]=contourf(ax1,X1000,Y1000,disppot,30);
			hold on
			set(curve1,'linecolor','none');
			c1=colorbar(ax1);
			c1.Label.String= '\phi [V]';
			c1.FontSize=14;
			axis square
			xlabel('x [mm]')
			ylabel('y [mm]')
			set(ax1,'fontsize',14)

			% Draw the metallic boundaries and the geometry itself
			[c,hContour]=contourf(ax1,X1000,Y1000,-geomw,[0,0],'linewidth',1.5);
			drawnow;
			% Change the color of the metallic boundaries to grey
			hFills=hContour.FacePrims;
			[hFills.ColorType] = deal('truecoloralpha');
			hFills(end).ColorData = uint8([150;150;150;255]);
			for idx = 1 : numel(hFills)-1
			hFills(idx).ColorData(4) = 0;
			end

			%% Radial electric field
			subplot(2,2,2);
			ax2=gca;
			[h2,curve2]=contourf(ax2,X1000,Y1000,dispEr,30);
			hold on
			set(curve2,'linecolor','none');
			c2=colorbar(ax2);
			c2.Label.String= 'E_r [V/m]';
			c2.FontSize=14;
			axis square
			xlabel('x [mm]')
			ylabel('y [mm]')
			set(ax2,'fontsize',14)

			% Draw the metallic boundaries and the geometry itself
			[c,hContour]=contourf(ax2,X1000,Y1000,-geomw,[0,0],'linewidth',1.5);
			drawnow;
			% Change the color of the metallic boundaries to grey
			hFills=hContour.FacePrims;
			[hFills.ColorType] = deal('truecoloralpha');
			hFills(end).ColorData = uint8([150;150;150;255]);
			for idx = 1 : numel(hFills)-1
			hFills(idx).ColorData(4) = 0;
			end

			%% Azimuthal electric field
			subplot(2,2,3);
			ax3=gca;
			[h3,curve3]=contourf(ax3,X1000,Y1000,dispEthet,30);
			hold on
			set(curve3,'linecolor','none');
			c3=colorbar(ax3);
			c3.Label.String= 'E_\theta [V/m]';
			c3.FontSize=14;
			axis square
			xlabel('x [mm]')
			ylabel('y [mm]')
			set(ax3,'fontsize',14)

			% Draw the metallic boundaries and the geometry itself
			[c,hContour]=contourf(ax3,X1000,Y1000,-geomw,[0,0],'linewidth',1.5);
			drawnow;
			% Change the color of the metallic boundaries to grey
			hFills=hContour.FacePrims;
			[hFills.ColorType] = deal('truecoloralpha');
			hFills(end).ColorData = uint8([150;150;150;255]);
			for idx = 1 : numel(hFills)-1
			hFills(idx).ColorData(4) = 0;
			end

			%% Axial electric field
			subplot(2,2,4);
			ax4=gca;
			[h4,curve4]=contourf(ax4,X1000,Y1000,dispEz,30);
			hold on
			set(curve4,'linecolor','none');
			c4=colorbar(ax4);
			c4.Label.String= 'E_z [V/m]';
			c4.FontSize=14;
			axis square
			xlabel('x [mm]')
			ylabel('y [mm]')
			set(ax4,'fontsize',14)

			% Draw the metallic boundaries and the geometry itself
			[c,hContour]=contourf(ax4,X1000,Y1000,-geomw,[0,0],'linewidth',1.5);
			drawnow;
			% Change the color of the metallic boundaries to grey
			hFills=hContour.FacePrims;
			[hFills.ColorType] = deal('truecoloralpha');
			hFills(end).ColorData = uint8([150;150;150;255]);
			for idx = 1 : numel(hFills)-1
			hFills(idx).ColorData(4) = 0;
			end
			end

			function displaypotmode(obj,modenumber,tstep,zlims,rlims)
			%% Plot the Fourier mode of the electrostatic potential for time t2d(tstep)
			if nargin<2
			modenumber=0;
			end
			if nargin<3
			tstep=1;
			end
			if nargin<4
			zlims=[obj.zgrid(1),obj.zgrid(end)]*1e3;
			end
			if nargin<5
			rlims=[obj.rgrid(1),obj.rgrid(end)]*1e3;
			end

			geomw=obj.geomweight(:,:,1);
			real=squeeze(obj.phireal(:,:,modenumber+1,tstep));
			imag=squeeze(obj.phiimag(:,:,modenumber+1,tstep));
			phim=real+1i*imag;
			mod=abs(phim);
			arg=angle(phim);
			mod(geomw<=0)=NaN;
			arg(geomw<=0)=NaN;

			f=figure('Name', 'Electrostatic potential','Color',[1 1 1],'Position', [300, 300, 800, 400]);
			ax1=gca;
			[h,curve]=contourf(ax1,obj.zgrid1000,obj.rgrid1000,mod,30);
			hold on
			set(curve,'linecolor','none');
			xlabel(ax1,'z [mm]')
			ylabel(ax1,'r [mm]')
			set(ax1,'fontsize',14)
			c = colorbar(ax1);
			c.Label.String= ['\Phi_',num2str(modenumber),' [V]'];
			c.FontSize=14;

			% Draw the metallic boundaries and the geometry itself
			[c1,hContour]=contourf(ax1,obj.zgrid1000,obj.rgrid1000,-geomw,[0,0],'linewidth',1.5);
			drawnow;
			xlim(ax1,zlims)
			ylim(ax1,rlims)
			% Change the color of the metallic boundaries to grey
			hFills=hContour.FacePrims;
			[hFills.ColorType] = deal('truecoloralpha');
			hFills(end).ColorData = uint8([150;150;150;255]);
			for idx = 1 : numel(hFills)-1
			hFills(idx).ColorData(4) = 0;
			end

			title(['t = ',num2str(obj.t2d(tstep)),' s']);
			end

			function displaypartspos(obj,tstep)
			if nargin<2
			tstep=1;
			end

			rpos=obj.species(1).R(1:obj.npart(1),tstep)*1000;
			thetpos=obj.species(1).THET(1:obj.npart(1),tstep);
			zpos=obj.species(1).Z(1:obj.npart(1),tstep)*1000;
			xpos=rpos.*cos(thetpos);
			ypos=rpos.*sin(thetpos);

			f=figure('Position', [100, 100, 1000, 600],'Color',[1 1 1]);
			subplot(1,2,1)
			ax1=gca;
			plot(ax1,xpos,ypos,'rx','Displayname','particle positions')
			axis square
			xlabel(ax1,'x [mm]')
			ylabel(ax1,'y [mm]')
			set(ax1,'fontsize',14)

			subplot(1,2,2)
			ax2=gca;
			plot(ax2,zpos,rpos,'rx','Displayname','particle positions')
			xlabel(ax2,'z [mm]')
			ylabel(ax2,'r [mm]')
			set(ax2,'fontsize',14)
			end

			function createvideorz(obj,thetindex,tsteps)
			%% Create a video of the time evolution of the density in the r-z plane
			if nargin<2
			thetindex=1;
			end
			if nargin<3
			tsteps=1:length(obj.t2d);
			end

			t=obj.t2d(tsteps);
			geomw=obj.geomweight(:,:,1);
			dens=squeeze(obj.N(:,thetindex,:,tsteps));
			pot=squeeze(obj.pot(:,thetindex,:,tsteps));
			%nmax=max(dens(:));
			%nmax=1e15;
			%nlvls=linspace(0,nmax,100);
			minN=1e14;
			maxN=1e17;
			nlvls=logspace(log10(minN),log10(maxN),50);

			video=VideoWriter('/home/pgiroud/fennecslite/figures/videorz');
			video.FrameRate=10;
			video.Quality=70;
			open(video)

			f=figure('Position', [100, 100, 1000, 600],'Color',[1 1 1]);
			ax1=gca;
			dispdens=dens(:,:,1);
			dispdens(geomw<=0)=NaN;
			[h,curve]=contourf(ax1,obj.zgrid1000,obj.rgrid1000,dispdens,nlvls,'Displayname','n_e [m^{-3}]');
			set(curve,'linecolor','none');
			colormap(flipud(hot))
			hold on
			c=colorbar(ax1);
			c.FontSize=14;
			c.Label.String= 'electron density [m^{-3}]';
			c.Label.FontSize=14;
			set(ax1,'ColorScale','log')
			caxis(ax1,[nlvls(1) nlvls(end)]);
			xlabel(ax1,'z [mm]')
			ylabel(ax1,'r [mm]')
			set(ax1,'fontsize',14)
			title(['t = ', num2str(t(1)),' s'])

			% Equipotential lines
			disppot=pot(:,:,1); %+obj.potxt;
			disppot(geomw<=0)=NaN;
			potcolor='b';
			[c1,h2]=contour(ax1,obj.zgrid1000,obj.rgrid1000,disppot/1e3,'--','color',potcolor,'ShowText','on','linewidth',1.2,'Displayname','Equipotentials [kV]');
			clabel(c1,h2,'Color',potcolor,'fontsize',14,'LabelSpacing',800)
			legend(h2,'Equipotentials [kV]','location','southwest')

			% Draw the metallic boundaries and the geometry itself
			geom=obj.geomweight(:,:,1);
			geom(geom<1e-9)=0;
			geom(geom>1e-9)=NaN;
			[c1,hContour]=contourf(ax1,obj.zgrid1000,obj.rgrid1000,geom, [0 0]);
			hold on
			drawnow;
			grid on;
			view(ax1,2)

			% Change the color of the metallic boundaries to grey
			hFills=hContour.FacePrims;
			[hFills.ColorType] = deal('truecoloralpha'); % default = 'truecolor'
			try
			hFills(1).ColorData = uint8([150;150;150;255]);
			for idx = 2 : numel(hFills)
			hFills(idx).ColorData(4) = 0; % default=255
			end
			catch
			end

			% Run the loop to update the electric potential and electron density
			for n=1:length(tsteps)
			dispdens=dens(:,:,n);
			dispdens(geomw<=0)=NaN;
			disppot=pot(:,:,n); %+obj.potxt;
			disppot(geomw<=0)=NaN;
			curve.ZData=dispdens;
			h2.ZData=disppot/1e3;
			caxis(ax1,[nlvls(1) nlvls(end)]);
			title(['t = ', num2str(t(n)),' s'])
			drawnow;
			currframe=getframe(gcf);
			writeVideo(video,currframe)
			end
			close(video)
			end

			function createvideorthet(obj,zindex,tsteps,interp)
			%% Create a video of the time evolution of the density in the r-theta plane
			if nargin<2
			zindex=obj.nz/2;
			end
			if nargin<3
			tsteps=1:length(obj.t2d);
			end
			if nargin<4
			interp=false;
			end

			t=obj.t2d(tsteps);
			r=obj.rgrid;
			thet=[obj.thetgrid' 2*pi];

			if interp
			density=zeros(length(r),length(thet),length(t));
			%potential=zeros(length(r),length(thet),length(t));
			for i=1:length(tsteps)
			density(:,:,i)=[squeeze(obj.N(:,:,zindex,tsteps(i))) squeeze(obj.N(:,1,zindex,tsteps(i)))];
			%potential(:,:,i)=[squeeze(obj.pot(:,:,zindex,tsteps(i))) squeeze(obj.pot(:,1,zindex,tsteps(i)))];
			end
			thetinterp=linspace(0,2*pi,100);
			dens=zeros(length(r),length(thetinterp),length(t));
			%pot=zeros(length(r),length(thetinterp),length(t));
			for i=1:size(dens,1)
			for j=1:size(dens,3)
			dens(i,:,j)=interp1(thet,density(i,:,j),thetinterp);
			%pot(i,:,j)=interp1(thet,potential(i,:,j),thetinterp);
			end
			end
			thet=thetinterp;
			else
			dens=zeros(length(r),length(thet),length(t));
			for i=1:length(tsteps)
			dens(:,:,i)=[squeeze(obj.N(:,:,zindex,tsteps(i))) squeeze(obj.N(:,1,zindex,tsteps(i)))];
			%pot(:,:,i)=[squeeze(obj.pot(:,:,zindex,tsteps(i))) squeeze(obj.pot(:,1,zindex,tsteps(i)))];
			end
			end

			dens0=dens(:,:,1);
			%nmax=max(dens0(:));
			%nlvls=linspace(0,nmax,100);
			minN=1e14;
			maxN=1e17;
			nlvls=logspace(log10(minN),log10(maxN),50);
			geomw=zeros(size(dens0));
			for i=1:size(dens,2)
			geomw(:,i)=squeeze(obj.geomweight(:,zindex,1));
			end

			[R,THET]=meshgrid(r,thet);
			X=R.*cos(THET);
			Y=R.*sin(THET);

			video=VideoWriter('/home/pgiroud/fennecslite/figures/videorthet');
			video.FrameRate=10;
			video.Quality=70;
			open(video)

			f=figure('Position', [100, 100, 700, 600],'Color',[1 1 1]);
			ax1=gca;
			dispdens=dens(:,:,1);
			dispdens(geomw<=0)=NaN;
			[h,curve]=contourf(ax1,X1000,Y1000,dispdens',nlvls,'Displayname','n_e [m^{-3}]');
			hold on
			set(curve,'linecolor','none');
			colormap(flipud(hot))
			c=colorbar(ax1);
			c.FontSize=14;
			c.Label.String= 'electron density [m^{-3}]';
			c.Label.FontSize=14;
			set(ax1,'ColorScale','log')
			caxis(ax1,[nlvls(1) nlvls(end)]);
			axis equal
			xlabel(ax1,'x [mm]')
			ylabel(ax1,'y [mm]')
			set(ax1,'fontsize',14)
			title(['t = ', num2str(t(1)),' s'])

			% Equipotential lines
			% disppot=pot(:,:,1)';
			% disppot(geomw<=0)=NaN;
			% potcolor='b';
			% [c1,h2]=contour(ax1,X1000,Y1000,disppot/1e3,'--','color',potcolor,'ShowText','on','linewidth',1.5,'Displayname','Equipotentials [kV]');
			% clabel(c1,h2,'Color',potcolor,'fontsize',14,'LabelSpacing',800)
			% legend(h2,'Equipotentials [kV]','location','southwest')

			% Draw the metallic boundaries and the geometry itself
			[c1,hContour]=contourf(ax1,X1000,Y1000,-geomw',[0,0],'linewidth',1.5);
			drawnow;

			% Change the color of the metallic boundaries to grey
			hFills=hContour.FacePrims;
			[hFills.ColorType] = deal('truecoloralpha');
			hFills(end).ColorData = uint8([150;150;150;255]);
			for idx = 1 : numel(hFills)-1
			hFills(idx).ColorData(4) = 0;
			end

			% Run the loop to update the electric potential and electron density
			for n=1:length(tsteps)
			dispdens=dens(:,:,n);
			dispdens(geomw<=0)=NaN;
			%disppot=pot(:,:,n);
			%disppot(geomw<=0)=NaN;
			curve.ZData=dispdens';
			%h2.ZData=disppot';
			caxis(ax1,[nlvls(1) nlvls(end)]);
			title(['t = ', num2str(t(n)),' s'])
			drawnow
			currframe=getframe(gcf);
			writeVideo(video,currframe)
			end
			close(video)
			end

			function displayfluidvel(obj,thetindex,tstep,zlims,rlims)
			%% Plot the fluid velocity for time t2d(tstep) in the r-z plane
			if nargin<2
			thetindex=1;
			end
			if nargin<3
			tstep=1;
			end
			if nargin<4
			zlims=[obj.zgrid(1),obj.zgrid(end)]*1e3;
			end
			if nargin<5
			rlims=[obj.rgrid(1),obj.rgrid(end)]*1e3;
			end

			geomw=obj.geomweight(:,:,1);
			dens=squeeze(obj.N(:,thetindex,:,tstep));
			Ur=squeeze(obj.fluidUR(:,thetindex,:,tstep));
			Uthet=squeeze(obj.fluidUTHET(:,thetindex,:,tstep));
			Uz=squeeze(obj.fluidUZ(:,thetindex,:,tstep));
			dens(geomw<=0)=0;
			Ur(geomw<=0)=0;
			Uthet(geomw<=0)=0;
			Uz(geomw<=0)=0;

			f=figure('Name', 'Self-consistent field','Color',[1 1 1],'Position', [100, 100, 1200, 600]);

			%% Density
			subplot(2,2,1)
			ax0=gca;

			[h,curve]=contourf(ax0,obj.zgrid1000,obj.rgrid1000,dens,30);
			hold on
			set(curve,'linecolor','none');
			xlabel(ax0,'z [mm]')
			ylabel(ax0,'r [mm]')
			set(ax0,'fontsize',14);
			c = colorbar(ax0);
			c.Label.String= 'n_e [m^{-3}]';
			c.FontSize=14;
			title(['t = ',num2str(obj.t2d(tstep)),' s']);

			% Draw the metallic boundaries and the geometry itself
			[c1,hContour]=contourf(ax0,obj.zgrid1000,obj.rgrid1000,-geomw,[0,0],'linewidth',1.5);
			drawnow;
			xlim(ax0,zlims)
			ylim(ax0,rlims)
			% Change the color of the metallic boundaries to grey
			hFills=hContour.FacePrims;
			[hFills.ColorType] = deal('truecoloralpha');
			hFills(end).ColorData = uint8([150;150;150;255]);
			for idx = 1 : numel(hFills)-1
			hFills(idx).ColorData(4) = 0;
			end

			%% Radial velocity
			subplot(2,2,2)
			ax1=gca;

			[h,curve]=contourf(ax1,obj.zgrid1000,obj.rgrid1000,Ur,30);
			hold on
			set(curve,'linecolor','none');
			xlabel(ax1,'z [mm]')
			ylabel(ax1,'r [mm]')
			set(ax1,'fontsize',14);
			c = colorbar(ax1);
			c.Label.String= 'U_r [m/s]';
			c.FontSize=14;

			% Draw the metallic boundaries and the geometry itself
			[c1,hContour]=contourf(ax1,obj.zgrid1000,obj.rgrid1000,-geomw,[0,0],'linewidth',1.5);
			drawnow;
			xlim(ax1,zlims)
			ylim(ax1,rlims)
			% Change the color of the metallic boundaries to grey
			hFills=hContour.FacePrims;
			[hFills.ColorType] = deal('truecoloralpha');
			hFills(end).ColorData = uint8([150;150;150;255]);
			for idx = 1 : numel(hFills)-1
			hFills(idx).ColorData(4) = 0;
			end

			%% Azimuthal velocity
			subplot(2,2,3)
			ax2=gca;

			[h,curve]=contourf(ax2,obj.zgrid1000,obj.rgrid1000,Uthet,30);
			hold on
			set(curve,'linecolor','none');
			xlabel(ax2,'z [mm]')
			ylabel(ax2,'r [mm]')
			set(ax2,'fontsize',14);
			c = colorbar(ax2);
			c.Label.String= 'U_\theta [m/s]';
			c.FontSize=14;

			% Draw the metallic boundaries and the geometry itself
			[c1,hContour]=contourf(ax2,obj.zgrid1000,obj.rgrid1000,-geomw,[0,0],'linewidth',1.5);
			drawnow;
			xlim(ax2,zlims)
			ylim(ax2,rlims)
			% Change the color of the metallic boundaries to grey
			hFills=hContour.FacePrims;
			[hFills.ColorType] = deal('truecoloralpha');
			hFills(end).ColorData = uint8([150;150;150;255]);
			for idx = 1 : numel(hFills)-1
			hFills(idx).ColorData(4) = 0;
			end

			%% Axial velocity
			subplot(2,2,4)
			ax3=gca;

			[h,curve]=contourf(ax3,obj.zgrid1000,obj.rgrid1000,Uz,30);
			hold on
			set(curve,'linecolor','none');
			xlabel(ax3,'z [mm]')
			ylabel(ax3,'r [mm]')
			set(ax3,'fontsize',14);
			c = colorbar(ax3);
			c.Label.String= 'U_z [m/s]';
			c.FontSize=14;

			% Draw the metallic boundaries and the geometry itself
			[c1,hContour]=contourf(ax3,obj.zgrid1000,obj.rgrid1000,-geomw,[0,0],'linewidth',1.5);
			drawnow;
			xlim(ax3,zlims)
			ylim(ax3,rlims)
			% Change the color of the metallic boundaries to grey
			hFills=hContour.FacePrims;
			[hFills.ColorType] = deal('truecoloralpha');
			hFills(end).ColorData = uint8([150;150;150;255]);
			for idx = 1 : numel(hFills)-1
			hFills(idx).ColorData(4) = 0;
			end
			end

			function displaymodesprofile(obj,modes,zindex,tstep)
			if nargin<2
			modes=0:obj.nthet/2;
			end
			if nargin<3
			zindex=obj.nz/2;
			end
			if nargin<4
			tstep=1;
			end

			t=obj.t2d(tstep);
			r=obj.rgrid*1000;

			f1=figure;
			ax1=gca;

			for m=1:length(modes)
			real=squeeze(obj.phireal(:,zindex,modes(m)+1,tstep));
			imag=squeeze(obj.phiimag(:,zindex,modes(m)+1,tstep));
			phim=abs(real+1i*imag);
			plot(ax1,r,phim,'linewidth',2,'Displayname',['m = ',num2str(modes(m))]);
			hold on
			end
			legend('Location','Southeast')
			xlabel(ax1,'r [mm]');
			ylabel(ax1,'\|\phi_m\| [a.u.]');
			set(ax1,'fontsize',14);
			end

			function displaymodekz(obj,m,tstep)
			if nargin<2
			m=1;
			end
			if nargin<3
			tstep=1;
			end

			t=obj.t2d(tstep);
			Lz=obj.zgrid(end)-obj.zgrid(1);
			kzmin=pi/Lz;
			kz=kzmin*[1:obj.nz];

			phim=squeeze(obj.phireal(:,:,m+1,tstep))+1i*squeeze(obj.phiimag(:,:,m+1,tstep));
			phim=fft(phim,[],2);
			maxphim=zeros(length(kz),1);
			for i=1:length(kz)
			maxphim(i)=max(squeeze(abs(phim(:,i))));
			end

			f1=figure;
			ax1=gca;
			plot(ax1,kz,maxphim,'linewidth',2);
			xlabel(ax1,'k_z [m^{-1}]');
			ylabel(ax1,'log \|\phi_m\| [a.u.]');
			set(ax1,'fontsize',14);
			title(['m = ',num2str(m),', t = ',num2str(t)]);
			end

			function displaymodesmevolution(obj,modes,ikz,tsteps)
			% Display the growth of the Fourier modes of the potential for
			% an axial mode number kz and different azimuthal mode numbers
			if nargin<2
			modes=0:obj.nthet/2;
			end
			if nargin<3
			ikz=1;
			end
			if nargin<4
			tsteps=1:length(obj.t2d);
			end

			t=obj.t2d(tsteps);
			Lz=obj.zgrid(end)-obj.zgrid(1);
			kzmin=2*pi/Lz;
			kz=ikz*kzmin;
			gamma=zeros(length(modes),1);

			f1=figure;
			ax1=gca;
			for m=1:length(modes)
			phim=squeeze(obj.phireal(:,:,modes(m)+1,tsteps))+1i*squeeze(obj.phiimag(:,:,modes(m)+1,tsteps));
			phim=fft(phim,[],2);
			phimkz=abs(squeeze(phim(:,ikz+1,:)));
			maxphimkz=zeros(length(t),1);
			for n=1:length(t)
			maxphimkz(n)=max(squeeze(phimkz(:,n)));
			end
			logphi=log(maxphimkz);
			p=polyfit(t,logphi,1);
			gamma(m)=p(1)/2;
			plot(ax1,t,logphi,'linewidth',2,'Displayname',['m = ',num2str(modes(m))]);
			hold on
			end
			legend('Location','Southeast')
			xlabel(ax1,'t [s]');
			ylabel(ax1,'log \|\phi_{m,k_z}\| [a.u.]');
			set(ax1,'fontsize',14);
			title(['k_z = ',num2str(kz)])

			f2=figure;
			ax2=gca;
			plot(modes,gamma,'--o')
			xlabel(ax2,'m');
			ylabel(ax2,'\gamma [s^{-1}]');
			set(ax2,'fontsize',14);
			title(['k_z = ',num2str(kz)])
			end

			function displaymodeskzevolution(obj,modes,m,tsteps)
			% Display the growth of the Fourier modes of the potential for
			% an azimuthal mode number m and different axial mode numbers
			if nargin<2
			modes=0:obj.nz;
			end
			if nargin<3
			m=1;
			end
			if nargin<4
			tsteps=1:length(obj.t2d);
			end

			t=obj.t2d(tsteps);
			Lz=obj.zgrid(end)-obj.zgrid(1);
			kzmin=2*pi/Lz;
			kz=kzmin*modes;
			gamma=zeros(length(modes),1);

			f1=figure;
			ax1=gca;
			for i=1:length(modes)
			phim=squeeze(obj.phireal(:,:,m+1,tsteps))+1i*squeeze(obj.phiimag(:,:,m+1,tsteps));
			phim=fft(phim,[],2);
			phimkz=abs(squeeze(phim(:,modes(i)+1,:)));
			maxphimkz=zeros(length(t),1);
			for n=1:length(t)
			maxphimkz(n)=max(squeeze(phimkz(:,n)));
			end
			logphi=log(maxphimkz);
			p=polyfit(t,logphi,1);
			gamma(i)=p(1)/2;
			plot(ax1,t,logphi,'linewidth',2,'Displayname',['k_z = ',num2str(kz(i))]);
			hold on
			end
			legend('Location','Southeast')
			xlabel(ax1,'t [s]');
			ylabel(ax1,'log \|\phi_{m,k_z}\| [a.u.]');
			set(ax1,'fontsize',14);
			title(['m = ',num2str(m)])

			f2=figure;
			ax2=gca;
			plot(kz,gamma,'--o')
			xlabel(ax2,'k_z');
			ylabel(ax2,'\gamma [s^{-1}]');
			set(ax2,'fontsize',14);
			title(['m = ',num2str(m)])
			end

			function displaymodesevolution(obj,m,kz,tsteps)
			if nargin<2
			m=0:obj.nthet/2;
			end
			if nargin<3
			kz=0:obj.nz;
			end
			if nargin<4
			tsteps=1:length(obj.t2d);
			end

			t=obj.t2d(tsteps);
			Lz=obj.zgrid(end)-obj.zgrid(1);
			kzmin=2*pi/Lz;
			ikz=kz;
			kz=kzmin*kz;

			f1=figure;
			ax1=gca;

			gamma=zeros(length(m),length(kz));
			maxphimkz=zeros(length(t),1);
			for i=1:length(m)
			phim=squeeze(obj.phireal(:,:,m(i)+1,tsteps))+1i*squeeze(obj.phiimag(:,:,m(i)+1,tsteps));
			phim=fft(phim,[],2);
			for j=1:length(kz)
			phimkz=abs(squeeze(phim(:,ikz(j)+1,:)));
			for n=1:length(t)
			maxphimkz(n)=max(squeeze(phimkz(:,n)));
			end
			logphi=log(maxphimkz);
			p=polyfit(t,logphi,1);
			gamma(i,j)=p(1)/2;
			plot(ax1,t,logphi,'linewidth',2,'Displayname',['m = ',num2str(m(i)),', k_z = ',num2str(kz(j))]);
			hold on
			end
			end
			legend('Location','Southeast')
			xlabel(ax1,'t [s]');
			ylabel(ax1,'log \|\phi_{m,k_z}\| [a.u.]');
			set(ax1,'fontsize',14);

			f2=figure;
			ax2=gca;
			pclr=pcolor(m,kz,gamma');
			set(pclr, 'edgecolor','none');
			%shading interp
			c=colorbar(ax2);
			c.FontSize=14;
			c.Label.String= '\gamma [s^{-1}]';
			c.Label.FontSize=14;
			set(ax2,'ColorScale','log')
			xlabel(ax2,'m');
			ylabel(ax2,'k_z [m^{-1}]');
			set(ax2,'fontsize',14);
			end
			end
			end

fennecslite/matlab/@fennecslitehdf5/fennecslitehdf5.mcdf25dab51fdmaster

fennecslitehdf5.mView Options

fennecslite/matlab/@fennecslitehdf5/fennecslitehdf5.m
cdf25dab51fdmaster

fennecslitehdf5.m
View Options