85 lines
5.1 KiB
Matlab
85 lines
5.1 KiB
Matlab
function VDk = calculateVDk(this,Params,Transf,TransfRad,IncludeDDICutOff)
|
|
% makes the dipolar interaction matrix, size numel(Params.kr) * numel(Params.kz)
|
|
% Rmax and Zmax are the interaction cutoffs
|
|
% VDk needs to be multiplied by Cdd
|
|
% approach is that of Lu, PRA 82, 023622 (2010)
|
|
|
|
% == Calculating the DDI potential in Fourier space with appropriate cutoff == %
|
|
% Cylindrical (semianalytic)
|
|
% Cylindrical infinite Z, polarized along x (analytic)
|
|
% Spherical
|
|
% Custom Cylindrical
|
|
|
|
if IncludeDDICutOff
|
|
switch this.CutOffType
|
|
case 'Cylindrical' % Cylindrical (semianalytic)
|
|
Zcutoff = Params.Lz/2;
|
|
alph = acos((Transf.KX*sin(Params.theta)*cos(Params.phi)+Transf.KY*sin(Params.theta)*sin(Params.phi)+Transf.KZ*cos(Params.theta))./sqrt(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2));
|
|
alph(1) = pi/2;
|
|
|
|
% Analytic part of cutoff for slice 0<z<Zmax, 0<r<Inf Ronen, PRL 98, 030406 (2007)
|
|
cossq = cos(alph).^2;
|
|
VDk = cossq-1/3;
|
|
sinsq = 1 - cossq;
|
|
VDk = VDk + exp(-Zcutoff*sqrt(Transf.KX.^2+Transf.KY.^2)).*( sinsq .* cos(Zcutoff * Transf.KZ) - sqrt(sinsq.*cossq).*sin(Zcutoff * Transf.KZ) );
|
|
|
|
% Nonanalytic part
|
|
% For a cylindrical cutoff, we need to construct a kr grid based on the 3D parameters using Bessel quadrature
|
|
VDkNon = this.calculateNumericalHankelTransform(TransfRad.kr, TransfRad.kz, TransfRad.Rmax, Zcutoff);
|
|
|
|
% Interpolating the nonanalytic part onto 3D grid
|
|
fullkr = [-flip(TransfRad.kr)',TransfRad.kr'];
|
|
[KR,KZ] = ndgrid(fullkr,TransfRad.kz);
|
|
[KX3D,KY3D,KZ3D] = ndgrid(ifftshift(Transf.kx),ifftshift(Transf.ky),ifftshift(Transf.kz));
|
|
KR3D = sqrt(KX3D.^2 + KY3D.^2);
|
|
fullVDK = [flip(VDkNon',2),VDkNon']';
|
|
VDkNon = interpn(KR,KZ,fullVDK,KR3D,KZ3D,'spline',0); %Last argument is -1/3 for full VDk. 0 for nonanalytic piece?
|
|
VDkNon = fftshift(VDkNon);
|
|
|
|
VDk = VDk + VDkNon;
|
|
|
|
case 'CylindricalInfiniteZ' % Cylindrical infinite Z, polarized along x -- PRA 107, 033301 (2023)
|
|
alph = acos((Transf.KX*sin(Params.theta)*cos(Params.phi)+Transf.KY*sin(Params.theta)*sin(Params.phi)+Transf.KZ*cos(Params.theta))./sqrt(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2));
|
|
alph(1) = pi/2;
|
|
rhoc = max([abs(Transf.x),abs(Transf.y)]);
|
|
KR = sqrt(Transf.KX.^2+Transf.KY.^2);
|
|
func = @(n,u,v) v.^2./(u.^2+v.^2).*(v.*besselj(n,u).*besselk(n+1,v) - u.*besselj(n+1,u).*besselk(n,v));
|
|
VDk = -0.5*func(0,KR*rhoc,abs(Transf.KZ)*rhoc) + (Transf.KX.^2./KR.^2 - 0.5).*func(2,KR*rhoc,abs(Transf.KZ)*rhoc);
|
|
VDk = (1/3)*(3*(cos(alph).^2)-1) - VDk;
|
|
|
|
VDk(KR==0) = -1/3 + 1/2*abs(Transf.KZ(KR==0))*rhoc.*besselk(1,abs(Transf.KZ(KR==0))*rhoc);
|
|
VDk(Transf.KZ==0) = 1/6 + (Transf.KX(Transf.KZ==0).^2-Transf.KY(Transf.KZ==0).^2)./(KR(Transf.KZ==0).^2).*(1/2 - besselj(1,KR(Transf.KZ==0)*rhoc)./(KR(Transf.KZ==0)*rhoc));
|
|
VDk(1,1,1) = 1/6;
|
|
|
|
case 'Spherical' % Spherical
|
|
Rcut = min(Params.Lx/2,Params.Ly/2,Params.Lz/2);
|
|
alph = acos((Transf.KX*sin(Params.theta)*cos(Params.phi)+Transf.KY*sin(Params.theta)*sin(Params.phi)+Transf.KZ*cos(Params.theta))./sqrt(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2));
|
|
alph(1) = pi/2;
|
|
|
|
K = sqrt(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
|
|
VDk = (cos(alph).^2-1/3).*(1 + 3*cos(Rcut*K)./(Rcut^2.*K.^2) - 3*sin(Rcut*K)./(Rcut^3.*K.^3));
|
|
|
|
case 'CustomCylindrical'
|
|
alph = acos((Transf.KX*sin(Params.theta)*cos(Params.phi)+Transf.KY*sin(Params.theta)*sin(Params.phi)+Transf.KZ*cos(Params.theta))./sqrt(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2));
|
|
alph(1) = pi/2;
|
|
VDk = cos(alph).^2-1/3;
|
|
% Implementing a custom cylindrical cutoff:
|
|
VDr = ifftn(VDk);
|
|
VDr = fftshift(VDr);
|
|
% Define the cylindrical mask
|
|
cylinder_mask = (Transf.X.^2 + Transf.Y.^2 <= this.CustomCylindricalCutOffRadius^2) & (abs(Transf.Z) <= this.CustomCylindricalCutOffHeight / 2);
|
|
% Multiply the mask to apply the cut-off
|
|
VDr = VDr.*double(cylinder_mask);
|
|
VDr = ifftshift(VDr);
|
|
VDk = fftn(VDr);
|
|
|
|
otherwise
|
|
disp('Choose a valid DDI cutoff type!')
|
|
return
|
|
end
|
|
else
|
|
alph = acos((Transf.KX*sin(Params.theta)*cos(Params.phi)+Transf.KY*sin(Params.theta)*sin(Params.phi)+Transf.KZ*cos(Params.theta))./sqrt(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2));
|
|
alph(1) = pi/2;
|
|
VDk = cos(alph).^2-1/3;
|
|
end
|
|
end |