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

fennecs/matlab/h5parts.m
a2847b9464fageom_testpardiso

Last Change
Raw File

h5parts.m
View Options

classdef h5parts
properties
fullpath
R
Z
THET
VR
VZ
VTHET
Rindex
Zindex
partindex
partepot
nlclassical
q
m
weight
tpart
it2
nbparts
rgrid
zgrid
vnorm
parent
end
methods
function obj=h5parts(filename,hdf5group,parent)
obj.fullpath=filename;
obj.nlclassical=parent.nlclassical;
obj.it2=parent.it2;
obj.rgrid=parent.rgrid;
obj.zgrid=parent.zgrid;
obj.vnorm=parent.vnorm;
obj.tpart = h5read(obj.fullpath,sprintf('%s/time',hdf5group));
obj.nbparts = h5read(obj.fullpath,sprintf('%s/Nparts',hdf5group));
obj.parent=parent;
obj.R = h5partsquantity(obj.fullpath,hdf5group,'R');
obj.Z = h5partsquantity(obj.fullpath,hdf5group,'Z');
obj.THET = h5partsquantity(obj.fullpath,hdf5group,'THET');
obj.VR = h5partsquantity(obj.fullpath,hdf5group,'UR',obj.vnorm);
obj.VZ = h5partsquantity(obj.fullpath,hdf5group,'UZ',obj.vnorm);
obj.VTHET= h5partsquantity(obj.fullpath,hdf5group,'UTHET',obj.vnorm);
obj.q= h5readatt(obj.fullpath,hdf5group,'q');
obj.m= h5readatt(obj.fullpath,hdf5group,'m');
obj.weight= h5readatt(obj.fullpath,hdf5group,'weight');
obj.partepot = h5partsquantity(filename,hdf5group,'pot',obj.q);
try
obj.Rindex=h5partsquantity(obj.fullpath,hdf5group,'Rindex');
obj.Zindex=h5partsquantity(obj.fullpath,hdf5group,'Zindex');
obj.partindex=h5partsquantity(obj.fullpath,hdf5group,'partindex');
catch
end
% try
% obj.partindex = h5read(obj.fullpath,sprintf('%s/partindex',hdf5group));
% partindex=obj.partindex;
% partindex(partindex==-1)=NaN;
% %[~,Indices]=sort(partindex,'ascend');
% Indices=obj.partindex;
% for i=1:size(Indices,2)
% obj.H(Indices(:,i),i)=obj.H(:,i);
% obj.P(Indices(:,i),i)=obj.P(:,i);
%
% obj.R(Indices(:,i),i)=obj.R(:,i);
% obj.Z(Indices(:,i),i)=obj.Z(:,i);
% obj.Rindex(Indices(:,i),i)=obj.Rindex(:,i);
% obj.Zindex(Indices(:,i),i)=obj.Zindex(:,i);
% obj.VR(Indices(:,i),i)=obj.VR(:,i);
% obj.VZ(Indices(:,i),i)=obj.VZ(:,i);
% obj.VTHET(Indices(:,i),i)=obj.VTHET(:,i);
% obj.THET(Indices(:,i),i)=obj.THET(:,i);
% end
% catch
% end
% clear partindex;
end
function quantity=H(obj,indices)
if indices{1}==':'
p=1:obj.VR.nparts;
else
p=indices{1};
end
if indices{2}==':'
t=1:length(obj.tpart);
else
t=indices{2};
end
if size(indices,1)>2
track=indices{3};
else
track=false;
end
quantity=0.5*obj.m*(obj.VR(p,t,track).^2+obj.VTHET(p,t,track).^2+obj.VZ(p,t,track).^2)+obj.partepot(p,t,track);
end
function quantity=P(obj,indices)
if indices{1}==':'
p=1:obj.R.nparts;
else
p=indices{1};
end
if indices{2}==':'
t=1:length(obj.tpart);
else
t=indices{2};
end
if size(indices,1)>2
track=indices{3};
else
track=false;
end
quantity=obj.R(p,t,track).*(obj.VTHET(p,t,track)*obj.m+obj.q*obj.parent.Atheta(obj.R(p,t,track),obj.Z(p,t,track)));
end
function quantity=Vpar(obj,indices)
%Vpar Computes the parallel velocity for the particle indices{1} at time indices{2}
if indices{1}==':'
p=1:obj.R.nparts;
else
p=indices{1};
end
if indices{2}==':'
t=1:length(obj.tpart);
else
t=indices{2};
end
if size(indices,1)>2
track=indices{3};
else
track=false;
end
Zind=obj.Zindex(p,t,track);
Rind=obj.Rindex(p,t,track);
posind=sub2ind(size(obj.parent.B),Zind,Rind);
costhet=obj.parent.Bz(posind)./obj.parent.B(posind);
sinthet=obj.parent.Br(posind)./obj.parent.B(posind);
quantity=obj.VR(p,t,track).*sinthet+obj.VZ(p,t,track).*costhet;
end
function quantity=Vperp(obj,indices)
%Vperp Computes the perpendicular velocity in the guidind center reference frame,
% for the main specie particle indices{1} at time indices{2}
if indices{1}==':'
p=1:obj.R.nparts;
else
p=indices{1};
end
if indices{2}==':'
t=1:length(obj.tpart);
else
t=indices{2};
end
if size(indices,1)>2
track=indices{3};
else
track=false;
end
Zind=obj.Zindex(p,t,track);
Rind=obj.Rindex(p,t,track);
posind=sub2ind(size(obj.parent.B),Zind,Rind);
costhet=obj.parent.Bz(posind)./obj.parent.B(posind);
sinthet=obj.parent.Br(posind)./obj.parent.B(posind);
Vdrift=zeros(size(Zind));
for j=1:length(t)
tfield=find(obj.parent.t2d==obj.tpart(t(j)));
timeEr=obj.parent.Er(:,:,tfield);
timeEz=obj.parent.Ez(:,:,tfield);
posindE=sub2ind(size(timeEr),Rind(:,j),Zind(:,j));
Vdrift(:,j)=(timeEz(posindE).*obj.Br(posind(:,j))-timeEr(posindE).*obj.Bz(posind(:,j)))./obj.B(posind(:,j)).^2;
end
quantity=sqrt((obj.VTHET(p,t,track)-Vdrift).^2+(obj.VR(p,t,track).*costhet-obj.VZ(p,t,track).*sinthet).^2);
end
function sref = subsref(obj,s)
% obj(i) is equivalent to obj.Data(i)
switch s(1).type
case '.'
if(strcmp(s(1).subs,'H'))
sref=H(obj,s(2).subs);
elseif(strcmp(s(1).subs,'P'))
sref=P(obj,s(2).subs);
else
sref=builtin('subsref',obj,s);
end
case '()'
sref=builtin('subsref',obj,s);
case '{}'
error('MYDataClass:subsref',...
'Not a supported subscripted reference')
end
end
end
end
c4science · Help