-
Notifications
You must be signed in to change notification settings - Fork 1
/
xdxlm.m
152 lines (143 loc) · 4.94 KB
/
xdxlm.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
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
function [X,dX,theta,dems]=xdxlm(l,m,theta,xver,tol,blox)
% [X,dX,theta,dems]=XDXLM(l,m,theta,xver,tol,blox)
%
% Calculates normalized (associate) Legendre functions DT (B.58/B.60)
% and their derivatives. This is a modification of the xlm function.
%
% INPUT:
%
% l degree (0 <= l <= infinity) [default: random]
% m order (-l <= m <= l) [default: all orders 0<=l]
% l and m can be vectors, but not both at the same time
% theta argument (0 <= theta <= pi) [default: 181 linearly spaced; not NaN!]
% xver 1 optional normalization check by Gauss-Legendre quadrature
% 0 no normalization check [default]
% tol Tolerance for optional normalization checking
% blox 0 Standard (lm) ordering, l=0:L, m=-l:l
% 1 Block-diagonal ordering, m=-L:L, l=abs(m):L
%
% OUTPUT:
%
% X The associated normalized Legendre function at the desired argument(s):
% as a scalar or a row vector with length(theta) columns, OR
% as a matrix with length(m) rows and length(theta) columns, OR
% as a matrix with length(l) rows and length(theta) columns, OR
% as a 3D matrix of size [length(l) size(theta), OR
% (L+1)^2 x length(theta) if you put in
% a degree l=[0 L] and an order []: lists orders -l to l.
% dX The derivative of X in exactly the same format
% theta The argument(s), which you might or not have specified
% dems The orders to which the Xlms belong (to verify input or block sorting)
%
% EXAMPLES:
%
% plot(xlm([0:5],0)')
% plot(xlm(5,[])')
%
% SEE ALSO:
%
% LIBBRECHT, PLM, YLM, BLM
%
% Last modified by fjsimons-at-alum.mit.edu, 01/18/2008
% Last modified by plattner-at-alumni.ethz.ch, 07/30/2012
% For single m, this will accept 2D theta's and return 3D results X
% See Research Notebook VI p 77ff
% Default values
defval('l',round(rand*10))
defval('m',[])
defval('theta',linspace(0,pi,181))
% Never make this one or it will reach internal recursion limit
defval('xver',0)
defval('tol',1e-10)
defval('blox',0)
if blox~=0 & blox~=1
error('Specify valid block-sorting option ''blox''')
end
% Revert back to cos(theta)
mu=cos(theta);
switch xver
case 0
% If the degrees go from 0 to some L and m is empty, know what to do
if min(l)==0 & max(l)>0 & isempty(m)
X=repmat(NaN,(max(l)+1)^2,length(theta));
dX=repmat(NaN,(max(l)+1)^2,length(theta));
for thel=0:max(l)
[theX,thedX]=xlm(thel,-thel:thel,theta,xver,tol);
X(thel^2+1:(thel+1)^2,:)=theX;
dX(thel^2+1:(thel+1)^2,:)=thedX;
end
[dems,dels,mz,blkm]=addmout(max(l));
if blox==1
X= X(blkm,:);
dX=dX(blkm,:);
dems=dems(blkm);
end
return
end
dems=m;
% Error handling common to PLM, XLM, YLM - note this resets defaults
[l,m,mu,xver,tol]=pxyerh(l,m,mu,xver,tol);
% Calculation for m>0
if prod(size(l))==1 & prod(size(m))==1 % SINGLE L AND M
% Note that Matlab 'sch' has the sqrt(2-(m==0)) in there so we get
% rid of it; this option also has gotten rid of the Condon-Shortley
% phase which we now need to put back in
[theX,thedX]=libbrecht(l,mu,'sch');
X=(-1)^m*sqrt(2*l+1)/sqrt(2-(m==0))/sqrt(4*pi)*...
rindeks(theX,abs(m)+1);
dX=(-1)^m*sqrt(2*l+1)/sqrt(2-(m==0))/sqrt(4*pi)*...
rindeks(thedX,abs(m)+1);
% This straight from the rule DT B.60
if m<0
X=(-1)^m*X;
dX=(-1)^m*dX;
end
elseif prod(size(l))==1 % SINGLE L MULTIPLE M
[theX,thedX]=libbrecht(l,mu,'sch');
X=sqrt(2*l+1)/sqrt(4*pi)*...
(rindeks(theX,abs(m)+1)'*...
diag(((-1).^m)./sqrt(2-(m==0))))';
dX=sqrt(2*l+1)/sqrt(4*pi)*...
(rindeks(thedX,abs(m)+1)'*...
diag(((-1).^m)./sqrt(2-(m==0))))';
for index=find(m<0)
X(index,:)=(-1)^m(index)* X(index,:);
dX(index,:)=(-1)^m(index)*dX(index,:);
end
elseif prod(size(m))==1 % MUTIPLE L SINGLE M
if min(size(mu))==1 % VECTOR ARGUMENT
X=repmat(NaN,[length(l) length(mu) 1]);
dX=repmat(NaN,[length(l) length(mu) 1]);
ini=1;
else % MATRIX ARGUMENT
X=repmat(NaN,[length(l) size(mu)]);
dX=repmat(NaN,[length(l) size(mu)]);
% The following is to bypass a dimensionality convention
[theX,thedX]=libbrecht(l,mu,'sch');
X(1,:,:)=shiftdim((-1)^m/sqrt(2-(m==0))/sqrt(4*pi)*...
theX,-1);
dX(1,:,:)=shiftdim((-1)^m/sqrt(2-(m==0))/sqrt(4*pi)*...
thedX,-1);
ini=2;
end
for index=ini:length(l)
[theX,thedX]=libbrecht(l(index),mu,'sch');
X(index,:,:)=(-1)^m*sqrt(2*l(index)+1)/sqrt(2-(m==0))/sqrt(4*pi)*...
rindeks(theX,abs(m)+1);
dX(index,:,:)=(-1)^m*sqrt(2*l(index)+1)/sqrt(2-(m==0))/sqrt(4*pi)*...
rindeks(thedX,abs(m)+1);
end
if m<0
X=(-1)^m* X;
dX=(-1)^m*dX;
end
else
error('Specify valid option')
end
case 1
% Check normalization... only for different 'ms
pxynrm(l,unique(m),tol,'X');
% Still also need to get you the answer at the desired argument
[X,theta,dems]=xlm(l,m,theta,0);
% So you can't return norms anymore - since you force not to
end