Skip to content

Commit

Permalink
Add files via upload
Browse files Browse the repository at this point in the history
  • Loading branch information
@waannuullll
waannuullll authored Jul 20, 2020
0 parents commit 3310e5c
Show file tree
Hide file tree
Showing 17 changed files with 1,303 additions and 0 deletions.
41 changes: 41 additions & 0 deletions A.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,41 @@
function At = A(m,n,mu,nu,kd,theta0,phi0)
format long

j = @(n,x) sqrt(pi/(2*x))*besselj(n+0.5,x); %spherical bessel
h = @(n,x) sqrt(pi/(2*x))*(besselj(n+0.5,x) + 1i*bessely(n+0.5,x)); %spherical Hankel
qmax = min([n,nu,0.5*(n+nu-abs(m-mu))]);
Qmax = min([n+1,nu,0.5*(n+1+nu-abs(m-mu))]);
%calculate Gaunt coefficient 0-->1
%calculate the denominator of a1q
%a1(1) = (pochhammer(n+1,n)*pochhammer(nu+1,nu)*factorial(n+nu--m-mu))...
% /(pochhammer(n+nu+1,n+nu)*factorial(n--m)*factorial(nu-mu));
a1(1) = 1;
a0 = a1(1);
%calculate the normalized a1q
for q=2:Qmax+1
p = n+nu-2*(q-1);
n4 = n+nu--m-mu;
S1 = 0;
for k=1:q
s1 = (poc(-m-n,2*k-2)*poc(mu-nu,2*q-2*k))...
/(factorial(k-1)*factorial(q-k)*poc(-n+0.5,k-1)*poc(-nu+0.5,q-k));
S1 = S1 + s1;
end
S2 = 0;
for j=1:q-1
s2 = a1(j)*poc(-p-q+j+0.5,q-j)/factorial(q-j);
S2 = S2 + s2;
end
a1(q) = poc(p+0.5,2*q-2)*S1/poc(-n4,2*q-2)-S2;
end
%calculate A coefficient
A1 = exp(1i*(mu-m)*(phi0))*((-1)^m*(1i)^(nu+n)*poc(n+2,n-1)*poc(nu+2,nu+1)*factorial(n+nu+m-mu))...
/(4*n*poc(n+nu+1,n+nu)*factorial(n-m)*factorial(nu+mu));
%A1 = (-1)^(m+n)*a0*((2*n+1)/(2*n*(n+1)))*exp(1i*(mu-m)*phi0);
A2 = 0;
for q=1:qmax+1
p = n+nu-2*(q-1);
sA = (-1)^(q-1)*(n*(n+1)+nu*(nu+1)-p*(p+1))*a1(q)*h(p,kd)*P(mu-m,p,cos(theta0));
A2 = A2 + sA;
end
At = A1*A2;
41 changes: 41 additions & 0 deletions A_.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,41 @@
function At = A_(m,n,mu,nu,kd,theta0,phi0)
format long

J = @(n,x) sqrt(pi/(2*x))*besselj(n+0.5,x); %spherical bessel
h = @(n,x) sqrt(pi/(2*x))*(besselj(n+0.5,x) + 1i*bessely(n+0.5,x)); %spherical Hankel
qmax = min([n,nu,0.5*(n+nu-abs(m-mu))]);
Qmax = min([n+1,nu,0.5*(n+1+nu-abs(m-mu))]);
%calculate Gaunt coefficient 0-->1
%calculate the denominator of a1q
%a1(1) = (pochhammer(n+1,n)*pochhammer(nu+1,nu)*factorial(n+nu--m-mu))...
% /(pochhammer(n+nu+1,n+nu)*factorial(n--m)*factorial(nu-mu));
a1(1) = 1;
a0 = a1(1);
%calculate the normalized a1q
for q=2:Qmax+1
p = n+nu-2*(q-1);
n4 = n+nu--m-mu;
S1 = 0;
for k=1:q
s1 = (poc(-m-n,2*k-2)*poc(mu-nu,2*q-2*k))...
/(factorial(k-1)*factorial(q-k)*poc(-n+0.5,k-1)*poc(-nu+0.5,q-k));
S1 = S1 + s1;
end
S2 = 0;
for j=1:q-1
s2 = a1(j)*poc(-p-q+j+0.5,q-j)/factorial(q-j);
S2 = S2 + s2;
end
a1(q) = poc(p+0.5,2*q-2)*S1/poc(-n4,2*q-2)-S2;
end
%calculate A coefficient
A1 = exp(1i*(mu-m)*(phi0))*((-1)^m*(1i)^(nu+n)*poc(n+2,n-1)*poc(nu+2,nu+1)*factorial(n+nu+m-mu))...
/(4*n*poc(n+nu+1,n+nu)*factorial(n-m)*factorial(nu+mu));
%A1 = (-1)^(m+n)*a0*((2*n+1)/(2*n*(n+1)))*exp(1i*(mu-m)*phi0);
A2 = 0;
for q=1:qmax+1
p = n+nu-2*(q-1);
sA = (-1)^(q-1)*(n*(n+1)+nu*(nu+1)-p*(p+1))*a1(q)*J(p,kd)*P(mu-m,p,cos(theta0));
A2 = A2 + sA;
end
At = A1*A2;
50 changes: 50 additions & 0 deletions At.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,50 @@
function Atot = At(m,n,mu,nu,k2d,tetalj,philj)
format long

%Bessel Functions
sphj= @(n,kr) sqrt(pi./(2*kr)).*besselj(n+0.5,kr);
sphy=@(n,kr) (sqrt(pi./(2*kr)).*(bessely(n+0.5,kr)));
sphh=@(n,kr) sphj(n,kr)+(1i.*sphy(n,kr));

%Translation coefficient function
Amn=@(m,n,mu,nu,phi)((((-1).^m).*((1i).^(nu+n)).*(gamma((2.*n)+1)./gamma(n+2)).*(gamma((2.*nu)+3)./gamma(nu+2)).*factorial(n+nu+m-mu))./(4.*n.*(gamma((2.*n)+(2.*nu)+1)./gamma(n+nu+1)).*factorial(n-m).*factorial(nu+mu)))*exp((1i).*(mu-m).*(philj+pi));
amn=@(n,nu,d,p,q,A,leg) ((-1).^q).*((n.*(n+1))+(nu.*(nu+1))-(p.*(p+1))).*A.*sphh(p,d).*leg;

%Gaunt coefficient function (Version1)
a1 = @(p,q) gamma(p+(2.*q)+0.5)./gamma(p+0.5);
a2 = @(m,n,mu,nu,q) gamma(n+nu+m-mu+1)./(((-1).^(2.*q)).*gamma(n+nu+m-mu-(2.*q)+1));
a3 = @(m,n,k) gamma(m+n+1)./(((-1).^(2.*k)).*gamma(m+n-2.*k+1));
a4 = @(mu,nu,q,k) gamma(nu-mu+1)./(((-1).^((2.*q)-(2.*k))).*gamma(nu-mu-(2.*q)+(2.*k)+1));
a5 = @(n,k) gamma(n+0.5)./(((-1).^k).*gamma(n-k+0.5));
a6 = @(nu,q,k) gamma(nu+0.5)./(((-1).^(q-k)).*gamma(nu-q+k+0.5));
a7 = @(p,q,j) gamma(p+q-j+0.5)./(((-1).^(q-j)).*gamma(p+0.5));

%Addition Coefficients
M = n+nu;
qmax = (n+nu-abs(m-mu))/2;
for i =0:1:qmax
pp(i+1) = M-(2.*i);
end
A = zeros(floor(qmax)+1,1);
A(1)=1;
for q=1:1:qmax
a12= a1(pp(q+1),q)./a2(m,n,mu,nu,q);
aa = 0;
bb = 0;
for k=0:1:q
aa = aa+((a3(m,n,k).*a4(mu,nu,q,k))./(factorial(k).*factorial(q-k).*a5(n,k).*a6(nu,q,k)));
end
for j=0:1:q-1
bb = bb +(((a7(pp(q+1),q,j))./factorial(q-j)).*A(j+1));
end
aq=(a12.*aa)-bb;
A(q+1) = aq;
end

asum = 0;
for q=0:1:qmax
leg = legendrePP(pp(q+1),mu-m,cos(tetalj));
asum = asum+ amn(n,nu,k2d,pp(q+1),q,A(q+1),leg);
end

Atot = Amn(m,n,mu,nu,philj).*asum;
69 changes: 69 additions & 0 deletions B.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,69 @@
function Bt = B(m,n,mu,nu,kd,theta0,phi0)

j = @(n,x) sqrt(pi/(2*x))*besselj(n+0.5,x); %spherical bessel
h = @(n,x) sqrt(pi/(2*x))*(besselj(n+0.5,x) + 1i*bessely(n+0.5,x)); %spherical Hankel
qmax = min([n,nu,0.5*(n+nu-abs(m-mu))]);
Qmax = min([n+1,nu,0.5*(n+1+nu-abs(m-mu))]);
% if m == -mu & n == nu
% Bt = 0;
% elseif m == -n & mu == -nu
% Bt = 0;
% elseif m == n & mu == nu
% Bt = 0;
% else
%calculate Gaunt coefficient 0-->1
%calculate the denominator of a2q and a3q
%a2(1) = (pochhammer(n+1+1,n+1)*pochhammer(nu+1,nu)*factorial(n+1+nu-(-m-1)-(mu+1)))...
% /(pochhammer(n+1+nu+1,n+1+nu)*factorial(n+1-(-m-1))*factorial(nu-(mu+1)));
%a3(1) = (pochhammer(n+1+1,n+1)*pochhammer(nu+1,nu)*factorial(n+1+nu--m-mu))...
% /(pochhammer(n+1+nu+1,n+1+nu)*factorial(n+1--m)*factorial(nu-mu));
a2(1) = 1; a3(1) = 1;
%b0 = ((2*n+1)*(n+nu+m-mu+1))/((2*n+2*nu+1)*(n+m+1));
%calculate the normalized a2q
for q=2:Qmax+1
p = n+1+nu-2*(q-1);
n4 = n+1+nu-(-m-1)-(mu+1);
S1 = 0;
for k=1:q
s1 = (poc(-m-1-(n+1),2*k-2)*poc(mu+1-nu,2*q-2*k))...
/(factorial(k-1)*factorial(q-k)*poc(-(n+1)+0.5,k-1)*poc(-nu+0.5,q-k));
S1 = S1 + s1;
end
S2 = 0;
for j=1:q-1
s2 = a2(j)*poc(-p-q+j+0.5,q-j)/factorial(q-j);
S2 = S2 + s2;
end
a2(q) = poc(p+0.5,2*q-2)*S1/poc(-n4,2*q-2)-S2;
end
%calculate the normalized a3q
for q=2:Qmax+1
p = n+1+nu-2*(q-1);
n4 = n+1+nu--m-mu;
S1 = 0;
for k=1:q
s1 = (poc(-m-(n+1),2*k-2)*poc(mu-nu,2*q-2*k))...
/(factorial(k-1)*factorial(q-k)*poc(-(n+1)+0.5,k-1)*poc(-nu+0.5,q-k));
S1 = S1 + s1;
end
S2 = 0;
for j=1:q-1
s2 = a3(j)*poc(-p-q+j+0.5,q-j)/factorial(q-j);
S2 = S2 + s2;
end
a3(q) = poc(p+0.5,2*q-2)*S1/poc(-n4,2*q-2)-S2;
end
%calculate B coefficient

B1 = exp(1i*(mu-m)*(phi0))*((-1)^m*(1i)^(nu+n+1)*poc(n+2,n+1)*poc(nu+2,nu+1)*factorial(n+nu+m-mu+1))...
/(4*n*(n+1)*(n+m+1)*poc(n+nu+2,n+nu+1)*factorial(n-m)*factorial(nu+mu));
%B1 = (-1)^(m+n)*1i*a0*b0*((2*n+1)/(2*n*(n+1)))*exp(1i*(mu-m)*phi0);
B2 = 0;
for q=1:Qmax+1
p = n+nu-2*(q-1);
sB = (-1)^(q-1)*(2*(n+1)*(nu-mu)*a2(q)-(p*(p+3)-nu*(nu+1)-n*(n+3)-2*mu*(n+1))*a3(q))...
*h(p+1,kd)*P(mu-m,p+1,cos(theta0));
B2 = B2 + sB;
end
Bt = B1*B2;
end
69 changes: 69 additions & 0 deletions B_.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,69 @@
function Bt = B(m,n,mu,nu,kd,theta0,phi0)

J = @(n,x) sqrt(pi/(2*x))*besselj(n+0.5,x); %spherical bessel
h = @(n,x) sqrt(pi/(2*x))*(besselj(n+0.5,x) + 1i*bessely(n+0.5,x));; %spherical Hankel
qmax = min([n,nu,0.5*(n+nu-abs(m-mu))]);
Qmax = min([n+1,nu,0.5*(n+1+nu-abs(m-mu))]);
% if m == -mu & n == nu
% Bt = 0;
% elseif m == -n & mu == -nu
% Bt = 0;
% elseif m == n & mu == nu
% Bt = 0;
% else
%calculate Gaunt coefficient 0-->1
%calculate the denominator of a2q and a3q
%a2(1) = (pochhammer(n+1+1,n+1)*pochhammer(nu+1,nu)*factorial(n+1+nu-(-m-1)-(mu+1)))...
% /(pochhammer(n+1+nu+1,n+1+nu)*factorial(n+1-(-m-1))*factorial(nu-(mu+1)));
%a3(1) = (pochhammer(n+1+1,n+1)*pochhammer(nu+1,nu)*factorial(n+1+nu--m-mu))...
% /(pochhammer(n+1+nu+1,n+1+nu)*factorial(n+1--m)*factorial(nu-mu));
a2(1) = 1; a3(1) = 1;
%b0 = ((2*n+1)*(n+nu+m-mu+1))/((2*n+2*nu+1)*(n+m+1));
%calculate the normalized a2q
for q=2:Qmax+1
p = n+1+nu-2*(q-1);
n4 = n+1+nu-(-m-1)-(mu+1);
S1 = 0;
for k=1:q
s1 = (poc(-m-1-(n+1),2*k-2)*poc(mu+1-nu,2*q-2*k))...
/(factorial(k-1)*factorial(q-k)*poc(-(n+1)+0.5,k-1)*poc(-nu+0.5,q-k));
S1 = S1 + s1;
end
S2 = 0;
for j=1:q-1
s2 = a2(j)*poc(-p-q+j+0.5,q-j)/factorial(q-j);
S2 = S2 + s2;
end
a2(q) = poc(p+0.5,2*q-2)*S1/poc(-n4,2*q-2)-S2;
end
%calculate the normalized a3q
for q=2:Qmax+1
p = n+1+nu-2*(q-1);
n4 = n+1+nu--m-mu;
S1 = 0;
for k=1:q
s1 = (poc(-m-(n+1),2*k-2)*poc(mu-nu,2*q-2*k))...
/(factorial(k-1)*factorial(q-k)*poc(-(n+1)+0.5,k-1)*poc(-nu+0.5,q-k));
S1 = S1 + s1;
end
S2 = 0;
for j=1:q-1
s2 = a3(j)*poc(-p-q+j+0.5,q-j)/factorial(q-j);
S2 = S2 + s2;
end
a3(q) = poc(p+0.5,2*q-2)*S1/poc(-n4,2*q-2)-S2;
end
%calculate B coefficient

B1 = exp(1i*(mu-m)*(phi0))*((-1)^m*(1i)^(nu+n+1)*poc(n+2,n+1)*poc(nu+2,nu+1)*factorial(n+nu+m-mu+1))...
/(4*n*(n+1)*(n+m+1)*poc(n+nu+2,n+nu+1)*factorial(n-m)*factorial(nu+mu));
%B1 = (-1)^(m+n)*1i*a0*b0*((2*n+1)/(2*n*(n+1)))*exp(1i*(mu-m)*phi0);
B2 = 0;
for q=1:Qmax+1
p = n+nu-2*(q-1);
sB = (-1)^(q-1)*(2*(n+1)*(nu-mu)*a2(q)-(p*(p+3)-nu*(nu+1)-n*(n+3)-2*mu*(n+1))*a3(q))...
*J(p+1,kd)*P(mu-m,p+1,cos(theta0));
B2 = B2 + sB;
end
Bt = B1*B2;
end
83 changes: 83 additions & 0 deletions Bt.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,83 @@
function Btot = Bt(m,n,mu,nu,k2d,tetalj,philj)
format long

%Bessel Functions
sphj=@(n,kr) sqrt(pi./(2*kr)).*besselj(n+0.5,kr);
sphy=@(n,kr) (sqrt(pi./(2*kr)).*(bessely(n+0.5,kr)));
sphh=@(n,kr) sphj(n,kr)+(1i.*sphy(n,kr));

%Translation coefficient function
Bmn = @(m,n,mu,nu,phi) ((((-1).^m).*((1i).^(nu+n+1)).*poc(n+2,n+1).*poc(nu+2,nu+1).*factorial(n+nu+m-mu+1))./(4.*n.*(n+1).*(n+m+1).*poc(n+nu+2,n+nu+1).*factorial(n-m).*factorial(nu+mu))).*exp((1i).*(mu-m).*(philj+pi));
bmn = @(n,mu,nu,d,p,q,E,G,leg) ((-1).^q).*((2.*(n+1).*(nu-mu).*E)-(((p.*(p+3))-(nu.*(nu+1))-(n.*(n+3))-(2.*mu.*(n+1))).*G)).*sphh(p+1,d).*leg;
if m == -mu & n == nu
Btot = 0;
elseif m == -n & mu == -nu
Btot = 0;
elseif m == n & mu == nu
Btot = 0;
else
a1 = @(p,q) poc(p+1.5,2*q);
a2 = @(m,n,mu,nu,q) poc(-n-nu-m+mu-1,2*q);
a3 = @(m,n,k) poc(-m-n-2,2*k);
a4 = @(mu,nu,q,k) poc(mu-nu+1,(2*q)-(2*k));
a5 = @(n,k) poc(-n+0.5-1,k);
a6 = @(nu,q,k) poc(-nu+0.5,q-k);
a7 = @(p,q,j) poc(-p-q+j+0.5-1,q-j);

M = n+nu;
qmax = (n+nu-abs(m-mu))/2;
Qmax=(n+nu+1-abs(m-mu))/2;
for i =0:1:floor(Qmax)
pp(i+1) = M-(2.*i);
end

E = zeros(floor(Qmax)+1,1);
E(1)=1;
G = zeros(floor(Qmax)+1,1);
G(1)=1;

for q=0:1:floor(Qmax)
a12= a1(pp(q+1),q)./a2(m,n,mu,nu,q);
aa = 0;
bb = 0;
for k=0:1:q
aa = aa+((a3(m,n,k).*a4(mu,nu,q,k))./(factorial(k).*factorial(q-k).*a5(n,k).*a6(nu,q,k)));
end
for j=0:1:q-1
bb = bb +(((a7(pp(q+1),q,j))./factorial(q-j)).*E(j+1));
end
aq=(a12.*aa)-bb;
E(q+1) = aq;
end
%Gaunt coefficient function (Version1)
a11 = @(p,q) poc(p+1.5,2*q);
a21 = @(m,n,mu,nu,q) poc(-n-nu-m+mu-1,2*q);
a31 = @(m,n,k) poc(-m-n-1,2*k);
a41 = @(mu,nu,q,k) poc(mu-nu,(2*q)-(2*k));
a51 = @(n,k) poc(-n+0.5-1,k);
a61 = @(nu,q,k) poc(-nu+0.5,q-k);
a71 = @(p,q,j) poc(-p-q+j+0.5-1,q-j);

for q=0:1:floor(Qmax)
a12= a11(pp(q+1),q)./a21(m,n,mu,nu,q);
aa = 0;
bb = 0;
for k=0:1:q
aa = aa+((a31(m,n,k).*a41(mu,nu,q,k))./(factorial(k).*factorial(q-k).*a51(n,k).*a61(nu,q,k)));
end
for j=0:1:q-1
bb = bb +(((a71(pp(q+1),q,j))./factorial(q-j)).*G(j+1));
end
aq=(a12.*aa)-bb;
G(q+1) = aq;
end

bsum = 0;
for q=0:1:floor(Qmax);
leg = legendrePP(pp(q+1)+1,mu-m,cos(tetalj));
bsum = bsum+ bmn(n,mu,nu,k2d,pp(q+1),q,E(q+1),G(q+1),leg);
end

Bmn(m,n,mu,nu,philj);
Btot =Bmn(m,n,mu,nu,philj).*bsum;
end
Binary file added J_C.mat
Binary file not shown.
Binary file added Johnson.mat
Binary file not shown.
Loading

0 comments on commit 3310e5c

Please sign in to comment.