-
Notifications
You must be signed in to change notification settings - Fork 103
/
PlSchmidt.f95
99 lines (87 loc) · 2.77 KB
/
PlSchmidt.f95
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
subroutine PlSchmidt(p, lmax, z, exitstatus)
!------------------------------------------------------------------------------
!
! This function evalutates all of the Schmidt normalized legendre
! polynomials up to degree lmax.
!
! Calling Parameters
!
! Out
! p A vector of all Schmidt normalized Legendgre polynomials
! evaluated at z up to lmax. The lenght must by greater or
! equal to (lmax+1).
!
! IN
! lmax Maximum degree to compute.
! z [-1, 1], cos(colatitude) or sin(latitude).
!
! OPTIONAL (OUT)
! exitstatus If present, instead of executing a STOP when an error
! is encountered, the variable exitstatus will be
! returned describing the error.
! 0 = No errors;
! 1 = Improper dimensions of input array;
! 2 = Improper bounds for input variable;
! 3 = Error allocating memory;
! 4 = File IO error.
!
! Notes:
!
! 1. The integral of Pl**2 over (-1,1) is 2 * / (2l+1).
! 2. The integral of Pl**2 over all space is 4 pi / (2l+1).
!
! Copyright (c) 2005-2019, SHTOOLS
! All rights reserved.
!
!------------------------------------------------------------------------------
use ftypes
implicit none
integer(int32), intent(in) :: lmax
real(dp), intent(out) :: p(:)
real(dp), intent(in) :: z
integer(int32), intent(out), optional :: exitstatus
real(dp) :: pm2, pm1, pl
integer(int32) :: l
if (present(exitstatus)) exitstatus = 0
if (size(p) < lmax+1) then
print*, "Error --- PlSchmidt"
print*, "P must be dimensioned as (LMAX+1) where LMAX is ", lmax
print*, "Input array is dimensioned ", size(p)
if (present(exitstatus)) then
exitstatus = 1
return
else
stop
end if
else if (lmax < 0) then
print*, "Error --- PlSchmidt"
print*, "LMAX must be greater than or equal to 0."
print*, "Input value is ", lmax
if (present(exitstatus)) then
exitstatus = 2
return
else
stop
end if
else if(abs(z) > 1.0_dp) then
print*, "Error --- PlSchmidt"
print*, "ABS(Z) must be less than or equal to 1."
print*, "Input value is ", z
if (present(exitstatus)) then
exitstatus = 2
return
else
stop
end if
end if
pm2 = 1.0_dp
p(1) = 1.0_dp
pm1 = z
p(2) = pm1
do l = 2, lmax, 1
pl = ( dble(2*l-1) * z * pm1 - dble(l-1) * pm2 ) / dble(l)
p(l+1) = pl
pm2 = pm1
pm1 = pl
end do
end subroutine PlSchmidt