-
Notifications
You must be signed in to change notification settings - Fork 106
/
SphericalCapCoef.f95
71 lines (58 loc) · 1.61 KB
/
SphericalCapCoef.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
subroutine SphericalCapCoef(coef, theta, lmax)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! This routine will return the coefficients for a spherical
! cap (located at the north pole) with angular radius theta,
! using the geodesy 4-pi normalization.
!
! Calling Parameters:
! IN
! theta Angular radius of spherical cap in RADIANS.
! OUT
! coef m=0 coefficients normalized according to the
! geodesy convention, further normalized such
! that the degree-0 term is 1. (i.e., function
! has an average value of one over the entire
! sphere).
! OPTIONAL
! lmax Maximum spherical harmonic degree.
!
! Dependencies: PlBar
!
! Written by Mark Wieczorek
!
! Copyright (c) 2005, Mark A. Wieczorek
! All rights reserved.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
use SHTOOLS, only: PlBar
implicit none
real*8, intent(out) :: coef(:)
real*8, intent(in) :: theta
integer, intent(in), optional :: lmax
real*8 :: x, top, bot, pi
real*8, allocatable :: pl(:)
integer :: l, lmax2, astat
pi = acos(-1.0d0)
coef = 0.0d0
x = cos(theta)
if (present(lmax) ) then
lmax2 = lmax
else
lmax2 = size(coef) - 1
endif
allocate(pl(lmax2+3), stat = astat)
if(astat /=0) then
print*, "Error --- SphericalCapCoef"
print*, "Unable to allocate array pl", astat
stop
endif
call PlBar(pl, lmax2+2, x)
coef(1) = 1.0d0
bot = pl(1) - pl(2) / sqrt(3.0d0)
do l=1, lmax2
top = pl(l) / sqrt(dble(2*l-1)) -pl(l+2) / sqrt(dble(2*l+3))
coef(l+1) = top / (bot * sqrt(dble(2*l+1)) )
enddo
deallocate (pl)
end subroutine SphericalCapCoef