-
Notifications
You must be signed in to change notification settings - Fork 103
/
SHSCouplingMatrixCap.f95
160 lines (131 loc) · 5.23 KB
/
SHSCouplingMatrixCap.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
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
152
153
154
155
156
157
158
159
160
subroutine SHSCouplingMatrixCap(kij, galpha, galpha_order, lmax, nmax, &
exitstatus)
!------------------------------------------------------------------------------
!
! SHSCouplingMatrixCap returns the spherical harmonic coupling matrix that
! relates the expectation of the function expressed in a subset of the best-
! localized spherical-cap Slepian functions to the expectation of the global
! power spectrum (assumed to be stationary). The Slepian functions are
! determined by a call to `SHReturnTapers`, and each row of galpha contains
! the (lmax+1) spherical harmonic coefficients for the single angular order
! as given in galpha_order.
!
! The relationship between the global and localized power spectra is given by
!
! < S_{\tilde{f}}(l) > = \sum_{l'=0}^lmax K_{ll'} S_{f}(l')
!
! where S_{\tilde{f}} is the expectation of the power spectrum at degree l of
! the function expressed in Slepian functions, S_{f}(l') is the expectation
! of the global power spectrum, and < ... > is the expectation operator.
! The coupling matrix is given explicitly by
!
! K_{ll'} = \frac{1}{2l'+1} Sum_{m=-mmax}^mmax
! ( Sum_{alpha=1}^nmax g_{l'm}(alpha) g_{lm}(alpha) )**2
!
! where mmax is min(l, lp).
!
! Calling Parameters
!
! IN
! galpha Matrix of dimension (LMAX+1, nmax) containing the
! spherical harmonic coefficients of the spherical-cap
! Slepian functions, obtained from a call to
! SHReturnTapers, where the functions are ordered from
! best to worst concentrated.
! galpha_order An array of dimension kmax containing the real
! angular order of the functions in galpha, as
! obtained from SHReturnTapers.
! lmax Maximum spherical harmonic degree of the Slepian
! functions.
! nmax Maximum number of Slepian coefficients used to
! construct the function.
!
! OUT
! kij Spherical harmonic coupling matrix of dimension
! ((LMAX+1)**2, nmax).
!
! 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.
!
! Copyright (c) 2005-2019, SHTOOLS
! All rights reserved.
!
!------------------------------------------------------------------------------
use ftypes
implicit none
real(dp), intent(out) :: kij(:,:)
real(dp), intent(in) :: galpha(:,:)
integer(int32), intent(in) :: galpha_order(:), lmax, nmax
integer(int32), intent(out), optional :: exitstatus
integer(int32) :: l, lp, m, alpha
real(dp) :: temp
if (present(exitstatus)) exitstatus = 0
if (size(kij(:,1)) < (lmax+1) .or. size(kij(1,:)) < (lmax+1)) then
print*, "Error --- SHSCouplingMatrixCap"
print*, "KIJ must be dimensioned as (LMAX+1, LMAX+1)."
print*, "LMAX = ", lmax
print*, "Dimension of KIJ = ", size(kij(:,1)), size(kij(1,:))
if (present(exitstatus)) then
exitstatus = 1
return
else
stop
end if
else if (size(galpha(:,1)) < (lmax+1) .or. &
size(galpha(1,:)) < nmax) then
print*, "Error --- SHSCouplingMatrixCap"
print*, "GALPHA must be dimensioned as (LMAX+1, NMAX )."
print*, "LMAX = ", lmax
print*, "NMAX = ", nmax
print*, "Input array is dimensioned as ", size(galpha(:,1)), &
size(galpha(1,:))
if (present(exitstatus)) then
exitstatus = 1
return
else
stop
end if
else if (size(galpha_order) < nmax) then
print*, "Error --- SHSCouplingMatrixCap"
print*, "GALPHA_ORDER must be dimensioned as (KMAX) where NMAX is ", &
nmax
print*, "Input array is dimensioned ", size(galpha_order)
if (present(exitstatus)) then
exitstatus = 1
return
else
stop
end if
end if
kij = 0.0_dp
! Calculate k(l, lp) * (2lp +1 ), which is symmetric
do l = 0, lmax
do lp = l, lmax, 1
do m = -l, l, 1
temp = 0.0_dp
do alpha = 1, nmax, 1
if (galpha_order(alpha) /= m) cycle
temp = temp + galpha(lp+1, alpha) * galpha(l+1, alpha)
end do
kij(l+1, lp+1) = kij(l+1, lp+1) + temp**2
end do
end do
end do
! Fill lower half of matrix
do l = 1, lmax, 1
do lp=0, l-1, 1
kij(l+1, lp+1) = kij(lp+1, l+1)
end do
end do
! Multiply by 1/(2lp+1)
do lp = 0, lmax
kij(1:lmax+1, lp+1) = kij(1:lmax+1, lp+1) / dble(2*lp+1)
end do
end subroutine SHSCouplingMatrixCap