-
Notifications
You must be signed in to change notification settings - Fork 0
/
mclpmn.m
123 lines (121 loc) · 4.71 KB
/
mclpmn.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
function mclpmn
%This program is a direct conversion of the corresponding Fortran program in
%S. Zhang & J. Jin "Computation of Special Functions" (Wiley, 1996).
%online: http://iris-lee3.ece.uiuc.edu/~jjin/routines/routines.html
%
%Converted by f2matlab open source project:
%online: https://sourceforge.net/projects/f2matlab/
% written by Ben Barrowes (barrowes@alum.mit.edu)
%
% ============================================================
% Purpose: This program computes the associated Legendre
% functions Pmn(z)and their derivatives Pmn'(z)for
% a complex argument using subroutine CLPMN
% Input : x --- Real part of z
% y --- Imaginary part of z
% m --- Order of Pmn(z), m = 0,1,2,...,n
% n --- Degree of Pmn(z), n = 0,1,2,...,N
% Output: CPM(m,n)--- Pmn(z)
% CPD(m,n)--- Pmn'(z)
% Examples:
% n = 5, x = 0.5, y = 0.2
% m Re[Pmn(z)]Im[Pmn(z)]Re[Pmn'(z)]Im[Pmn'(z)]
% -------------------------------------------------------------
% 0 .252594D+00 -.530293D+00 -.347606D+01 -.194250D+01
% 1 .333071D+01 .135206D+01 .117643D+02 -.144329D+02
% 2 -.102769D+02 .125759D+02 .765713D+02 .598500D+02
% 3 -.572879D+02 -.522744D+02 -.343414D+03 .147389D+03
% 4 .335711D+03 -.389151D+02 -.226328D+03 -.737100D+03
% 5 -.461125D+03 .329122D+03 .187180D+04 .160494D+02
% n = 5, x = 2.5, y = 1.0
% m Re[Pmn(z)]Im[Pmn(z)]Re[Pmn'(z)]Im[Pmn'(z)]
% -------------------------------------------------------------
% 0 -.429395D+03 .900336D+03 -.350391D+02 .193594D+04
% 1 -.216303D+04 .446358D+04 -.208935D+03 .964685D+04
% 2 -.883477D+04 .174005D+05 -.123703D+04 .381938D+05
% 3 -.273211D+05 .499684D+05 -.568080D+04 .112614D+06
% 4 -.565523D+05 .938503D+05 -.167147D+05 .219713D+06
% 5 -.584268D+05 .863328D+05 -.233002D+05 .212595D+06
% ============================================================
m=[];n=[];x=[];y=[];cpm=[];cpd=[];
cpm=zeros(40+1,40+1);
cpd=zeros(40+1,40+1);
fprintf(1,'%s \n',' please enter m, n, x and y ');
% READ(*,*)M,N,X,Y
m=0;
n=5;
x=.5;
y=.2;
fprintf(1,[repmat(' ',1,1),'m =','%2g',', ','n =','%2g',', ','x =','%5.1g',', ','y =','%5.1g' ' \n'],m,n,x,y);
[dumvar1,m,n,x,y,cpm,cpd]=clpmn(40,m,n,x,y,cpm,cpd);
fprintf(1,'%s ',' m n re[pmn(z)]im[pmn(z)]');fprintf(1,'%s \n','re[pmn''(z)]im[pmn''(z)]');
fprintf(1,'%s ',' -----------------------------------');fprintf(1,'%s \n','-------------------------------');
for j=0:n;
fprintf(1,[repmat(' ',1,1),repmat('%4g',1,2),repmat(' ',1,1),repmat('%14.6g',1,2),repmat(' ',1,1),repmat('%14.6g',1,2) ' \n'],m,j,cpm(m+1,j+1),cpd(m+1,j+1));
end; j=n+1;
%format(1x,2i4,1x,2d14.6,1x,2d14.6);
%format(1x,',i2,', ',',i2,', ',',f5.1, ', ',',f5.1);
end
function [mm,m,n,x,y,cpm,cpd]=clpmn(mm,m,n,x,y,cpm,cpd,varargin);
% =========================================================
% Purpose: Compute the associated Legendre functions Pmn(z)
% and their derivatives Pmn'(z)for a complex
% argument
% Input : x --- Real part of z
% y --- Imaginary part of z
% m --- Order of Pmn(z), m = 0,1,2,...,n
% n --- Degree of Pmn(z), n = 0,1,2,...,N
% mm --- Physical dimension of CPM and CPD
% Output: CPM(m,n)--- Pmn(z)
% CPD(m,n)--- Pmn'(z)
% =========================================================
z=complex(x,y);
for i=0:n;
for j=0:m;
cpm(j+1,i+1)=complex(0.0d0,0.0d0);
cpd(j+1,i+1)=complex(0.0d0,0.0d0);
end; j=fix(m)+1;
end; i=fix(n)+1;
cpm(0+1,0+1)=complex(1.0d0,0.0d0);
if(abs(x)== 1.0d0&y == 0.0d0);
for i=1:n;
cpm(0+1,i+1)=x.^i;
cpd(0+1,i+1)=0.5d0.*i.*(i).*x.^(i);
end; i=fix(n)+1;
for j=1:n;
for i=1:m;
if(i == 1);
cpd(i+1,j+1)=complex(1.0d+300,0.0d0);
elseif(i == 2);
cpd(i+1,j+1)=-0.25d0.*(j+2).*(j).*j.*(j-1).*x.^(j);
end;
end; i=fix(m)+1;
end; j=fix(n)+1;
return;
end;
ls=1;
if(abs(z)> 1.0d0)ls=-1; end;
zq=sqrt(ls.*(1.0d0-z.*z));
zs=ls.*(1.0d0-z.*z);
for i=1:m;
cpm(i+1,i+1)=-ls.*(2.0d0.*i-1.0d0).*zq.*cpm(i-1+1,i-1+1);
end; i=fix(m)+1;
for i=0:m;
cpm(i+1,i+1+1)=(2.0d0.*i+1.0d0).*z.*cpm(i+1,i+1);
end; i=fix(m)+1;
for i=0:m;
for j=i+2:n;
cpm(i+1,j+1)=((2.0d0.*j-1.0d0).*z.*cpm(i+1,j-1+1)-(i+j- 1.0d0).*cpm(i+1,j-2+1))./(j-i);
end; j=fix(n)+1;
end; i=fix(m)+1;
cpd(0+1,0+1)=complex(0.0d0,0.0d0);
for j=1:n;
cpd(0+1,j+1)=ls.*j.*(cpm(0+1,j-1+1)-z.*cpm(0+1,j+1))./zs;
end; j=fix(n)+1;
for i=1:m;
for j=i:n;
cpd(i+1,j+1)=ls.*i.*z.*cpm(i+1,j+1)./zs+(j+i).*(j-i+1.0d0)./zq.*cpm(i-1+1,j+1);
end; j=fix(n)+1;
end; i=fix(m)+1;
return;
end