-
Notifications
You must be signed in to change notification settings - Fork 1
/
spharm.m
131 lines (100 loc) · 3.19 KB
/
spharm.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
function [Y,Ytheta,Yphi] = spharm(n,m,theta,phi)
% spharm.m : scalar spherical harmonics and
% angular partial derivatives for given n,m (can take vector m).
%
% Usage:
% Y = spharm(n,m,theta,phi)
% or
% [Y,dY/dtheta,1/sin(theta)*dY/dphi] = spharm(n,m,theta,phi)
% or
% Y = spharm(n,theta,phi)
% or
% [Y,dtY,dpY] = spharm(n,theta,phi)
%
% Scalar n for the moment.
%
% If scalar m is used Y is a vector of length(theta,phi) and is
% completely compatible with previous versions of the toolbox. If vector m
% is present the output will be a matrix with rows of length(theta,phi) for
% m columns.
%
% "Out of range" n and m result in return of Y = 0
%
% PACKAGE INFO
%
% M.T. Edit: removed use of the custom legendre function, replaced it with
% the default Matlab legendre. This is more numerically stable.
if length(n)>1
error('n must be a scalar at present')
end
if nargin<4
phi=theta;
theta=m;
m=[-n:n];
end
%this is a cop out meant for future versions.
if nargout>1
mi=m;
m=[-n:n];
end
m=m(abs(m)<=n);
[theta,phi] = matchsize(theta,phi);
input_length = length(theta);
% if abs(m) > n | n < 0
% Y = zeros(input_length,1);
% Ytheta = zeros(input_length,1);
% Yphi = zeros(input_length,1);
% return
% end
% pnm = legendrerow(n,theta);
pnm=legendre(n,cos(theta),'norm')/sqrt(2*pi);
%pnm = pnm(abs(m)+1,:).';
% Why is this needed? Better do it, or m = 0 square integrals
% are equal to 1/2, not 1.
% This is either a bug in legendre or a mistake in the docs for it!
% Check this if MATLAB version changes! (Version 5.X)
%pnm(1,:) = pnm(1,:) * sqrt(2);
pnm = pnm(abs(m)+1,:); %pick the m's we potentially have.
[phiM,mv]=meshgrid(phi,m);
pnm = [(-1).^mv(m<0,:).*pnm(m<0,:);pnm(m>=0,:)];
expphi = exp(1i*mv.*phiM);
%N = sqrt((2*n+1)/(8*pi));
Y = pnm .* expphi;
% Do we want to calculate the derivatives?
if nargout <= 1
Y=Y.';
% Doesn't look like it
return
end
% We use recursion relations to find the derivatives, choosing
% ones that don't involve division by sin or cos, so no need to
% special cases to avoid division by zero
% theta derivative
% d/dtheta Y(n,m) = 1/2 exp(-i phi) sqrt((n-m)(n+m+1)) Y(n,m+1)
% - 1/2 exp(i phi) sqrt((n-m+1)(n+m)) Y(n,m-1)
ymplus=[pnm(2:end,:);zeros(1,length(theta))];
ymminus=[zeros(1,length(theta));pnm(1:end-1,:)];
Ytheta = (expphi).*(sqrt((n-mv+1).*(n+mv))/2 .* ymminus ...
- sqrt((n-mv).*(n+mv+1))/2 .* ymplus);
% phi derivative - actually 1/sin(theta) * d/dphi Y(n,m)
% Note that this is just i*m/sin(theta) * Y(n,m), but we use a
% recursion relation to avoid divide-by-zero trauma.
% 1/sin(theta) d/dphi Y(n,m) =
% i/2 * [ exp(-i phi) sqrt((2n+1)(n+m+1)(n+m+2)/(2n+3)) Y(n+1,m+1)
% + exp(i phi) sqrt((2n+1)(n-m+1)(n-m+2)/(2n+3)) Y(n+1,m-1) ]
Y2 = spharm(n+1,theta,phi).';
ymplus=Y2(3:end,:);
ymminus=Y2(1:end-2,:);
% exp(i*phi),exp(-i*phi) are useful for both partial derivatives
expplus = exp(1i*phiM);
expminus = exp(-1i*phiM);
% size(ymplus)
% size(mv)
% size(expminus)
Yphi = 1i/2 * sqrt((2*n+1)/(2*n+3)) * ...
( sqrt((n+mv+1).*(n+mv+2)) .* expminus .* ymplus ...
+ sqrt((n-mv+1).*(n-mv+2)) .* expplus .* ymminus );
Y=Y(n+mi+1,:).';
Yphi=Yphi(n+mi+1,:).';
Ytheta=Ytheta(n+mi+1,:).';
return