39 lines
1.2 KiB
Mathematica
39 lines
1.2 KiB
Mathematica
|
function VDkSemi = calculateNumericalHankelTransform(~,kr,kz,Rmax,Zmax,Nr)
|
||
|
|
|
||
|
|
% accuracy inputs for numerical integration
|
||
|
|
if(nargin==5)
|
||
|
|
Nr = 5e4;
|
||
|
|
end
|
||
|
|
|
||
|
|
Nz = 64;
|
||
|
|
farRmultiple = 2000;
|
||
|
|
|
||
|
|
% midpoint grids for the integration over 0<z<Zmax, Rmax<r<farRmultiple*Rmax (i.e. starts at Rmax)
|
||
|
|
dr=(farRmultiple-1)*Rmax/Nr;
|
||
|
|
r = ((1:Nr)'-0.5)*dr+Rmax;
|
||
|
|
dz=Zmax/Nz;
|
||
|
|
z = ((1:Nz)-0.5)*dz;
|
||
|
|
[R, Z] = ndgrid(r,z);
|
||
|
|
Rsq = R.^2 + Z.^2;
|
||
|
|
|
||
|
|
% real space interaction to be transformed
|
||
|
|
igrandbase = (1 - 3*Z.^2./Rsq)./Rsq.^(3/2);
|
||
|
|
|
||
|
|
% do the Hankel/Fourier-Bessel transform numerically
|
||
|
|
% prestore to ensure each besselj is only calculated once
|
||
|
|
% cell is faster than (:,:,krn) slicing
|
||
|
|
Nkr = numel(kr);
|
||
|
|
besselr = cell(Nkr,1);
|
||
|
|
for krn = 1:Nkr
|
||
|
|
besselr{krn} = repmat(r.*besselj(0,kr(krn)*r),1,Nz);
|
||
|
|
end
|
||
|
|
|
||
|
|
VDkSemi = zeros([Nkr,numel(kz)]);
|
||
|
|
for kzn = 1:numel(kz)
|
||
|
|
igrandbasez = repmat(cos(kz(kzn)*z),Nr,1) .* igrandbase;
|
||
|
|
for krn = 1:Nkr
|
||
|
|
igrand = igrandbasez.*besselr{krn};
|
||
|
|
VDkSemi(krn,kzn) = VDkSemi(krn,kzn) - sum(igrand(:))*dz*dr;
|
||
|
|
end
|
||
|
|
end
|
||
|
|
end
|