Homec4science
Diffusion FENNECS Lite (master) Browse

fennecslite/matlab/@fennecslitehdf5/calctrapvolume.m
7a743c8ce446master

Last Change
Raw File

calctrapvolume.m
View Options

function [Vol,VolV]= calctrapvolume(obj,time2d,mag,nblv)
% Calculates the Volume of potential well Vol in cm3
% and the total well volume ponderated by the well depth in eV cm3
if iscell(time2d)
time2d=cell2mat(time2d);
end
if nargin<3
mag=[];
end
if nargin<4
nblv=400;
end
% We prepare the geometric variables defining the boundaries
[zgrid,rgrid]=ndgrid(obj.zgrid,obj.rgrid);
geomw=obj.geomweight(:,:,1);
inside=griddedInterpolant(zgrid,rgrid,obj.geomweight(:,:,1)');
% Compute the Larmor radius assuming v_perp=v_exb
if time2d>0
Er=squeeze(obj.Er(:,1,:,time2d));
Ez=squeeze(obj.Ez(:,1,:,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)');
% Remove the points outside the geometry
rl(obj.geomweight(:,:,1)<0)=0;
% calculate the normal vector to the magnetic field lines at each position
perpB=zeros(size(obj.sinthet,1),size(obj.sinthet,2),2);
perpB(:,:,1)=sign(-Er).*obj.sinthet;
perpB(:,:,2)=-sign(-Er).*obj.costhet;
% Calculates at each position pos the upper limit of the Larmor motion posu
pos=zeros(size(obj.sinthet,1),size(obj.sinthet,2),2);
pos(:,:,1)=rgrid';
pos(:,:,2)=zgrid';
posu=pos+2*rl.*perpB;
isinside=inside(posu(:,:,2),posu(:,:,1))>0;
% Get the potential well at this time
model=obj.potentialwellmodel(time2d,true,mag,nblv);
z=model.z;
r=model.r;
Pot=model.pot;
Potdepth=scatteredInterpolant(z',r',Pot);
% calculate the effective potential well as seen by a particle along a
% larmor motion
effectdepth=0.5*(Potdepth(posu(:,:,2),posu(:,:,1))+Potdepth(pos(:,:,2),pos(:,:,1)));
effectdepth(geomw<0)=NaN;
effectdepth(~isinside)=NaN;
effectdepth(effectdepth<0)=NaN;
% Integrate the volume in cm3
rpi=pi*rgrid;
rpi(isnan(effectdepth))=0;
Vol=trapz(obj.rgrid,trapz(obj.zgrid,rpi,1),2)*1e6;
% Integrate the ponderated volume in eV cm3
effectdepth(isnan(effectdepth))=0;
VolV=trapz(obj.rgrid,trapz(obj.zgrid,effectdepth,2),1)*1e6;
end
c4science · Help