/
SHRotateCoef.f95
176 lines (149 loc) · 5.71 KB
/
SHRotateCoef.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
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
subroutine SHRotateCoef(x, cof, rcof, dj, lmax, exitstatus)
!------------------------------------------------------------------------------
!
! This subroutine will rotate a set of spherical harmonic coefficients (see
! documentation in SHRotateRealCoef for a description of how rotation angles
! are defined.) All angles are measured in radians, and spherical harmonic
! coefficients are fully normalized as in Edmonds, 1960. Note that Euler
! angles are different in Goldstein and Edmonds.
!
! This routine uses the "y-convention" where rotations are about the y axis
! instead of the x axis.
!
! djpi2 must be called before using this routine.
!
! Calling Parameters
!
! IN
! x Array of dimension 3 containing the three Euler angles
! in radians.
! dj Roation matrix with dimension (lmax+1, lmax+1, lmax+1).
! cof Indexed spherical harmonic coefficients with dimensions
! (2, (lmax+1)*(lmax+2)/2).
! lmax Maximum spherical harmonic degree.
!
! OUT
! rcof Indexed rotated spherical harmonic coefficients with
! dimensions (2, (lmax+1)*(lmax+2)/2).
!
! 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.
!
! History:
! 1. Based on routine from Guy Masters (July16, 1993)
! 2. Modified by Mark Simons (July 25, 1993)
! 3. Turned into readable f95 code by Mark Wieczorek (August, 2003)
!
! Copyright (c) 2005-2019, SHTOOLS
! All rights reserved.
!
!------------------------------------------------------------------------------
use ftypes
implicit none
real(dp), intent(in) :: cof(:,:), dj(:,:,:), x(3)
real(dp), intent(out) :: rcof(:,:)
integer(int32), intent(in) :: lmax
integer(int32), intent(out), optional :: exitstatus
real(dp) :: sum(2), temp(2,lmax+1), temp2(2,lmax+1), cgam(lmax+1), &
sgam(lmax+1), calf(lmax+1), salf(lmax+1), cbet(lmax+1), &
sbet(lmax+1), pi2, alpha, beta, gamma
integer(int32) :: ind, lp1, l, mp1, jp1, isgn, ii, indx
if (present(exitstatus)) exitstatus = 0
if (size(cof(:,1)) < 2 .or. size(cof(1,:)) < ((lmax+1)*(lmax+2))/2) then
print*, "Error --- SHRotateCoef"
print*, "COEF must be dimensioned (2, (LMAX+1)*(LMAX+2)/2) " // &
"where LMAX is ", lmax
print*, "Input array is dimensioned ", size(cof(:,1)), size(cof(1,:))
if (present(exitstatus)) then
exitstatus = 1
return
else
stop
end if
else if (size(rcof(:,1)) < 2 .or. size(rcof(1,:)) &
< ((lmax+1)*(lmax+2))/2) then
print*, "Error --- SHRotateCoef"
print*, "RCOEF must be dimensioned (2, (LMAX+1)*(LMAX+2)/2) " // &
"where LMAX is ", lmax
print*, "Input array is dimensioned ", size(rcof(:,1)), size(rcof(1,:))
if (present(exitstatus)) then
exitstatus = 1
return
else
stop
end if
else if (size(dj(:,1,1)) < lmax+1 .or. size(dj(1,:,1)) < lmax+1 &
.or. size(dj(1,1,:)) < lmax+1 ) then
print*, "Error --- SHRotateCoef"
print*, "DJ must be dimensioned (LMAX+1, LMAX+1, LMAX+1) " // &
"where LMAX is ", lmax
print*, "Input array is dimensioned ", size(dj(:,1,1)), &
size(dj(1,:,1)), size(dj(1,1,:))
if (present(exitstatus)) then
exitstatus = 1
return
else
stop
end if
end if
rcof = 0.0_dp
pi2 = 1.570796326794895_dp
alpha = x(1)
beta = x(2)
gamma = x(3)
alpha = alpha - pi2
gamma = gamma + pi2
beta = -beta
ind = 0
! Loop over degrees
do lp1 = 1, lmax + 1
l = lp1-1
cbet(lp1) = cos(l*beta)
sbet(lp1) = sin(l*beta)
cgam(lp1) = cos(l*gamma)
sgam(lp1) = sin(l*gamma)
calf(lp1) = cos(l*alpha)
salf(lp1) = sin(l*alpha)
! Alpha rotation
do mp1 = 1, lp1
indx = ind+mp1
temp(1,mp1) = cof(1,indx) * calf(mp1) - cof(2,indx) * salf(mp1)
temp(2,mp1) = cof(2,indx) * calf(mp1) + cof(1,indx) * salf(mp1)
end do
! B rotation and beta rotation
do jp1 = 1, lp1
sum(1) = dj(jp1,1,lp1) * temp(1,1)
sum(2) = 0.0_dp
isgn = 1 - 2 * mod((lp1-jp1),2)
do mp1 = 2, lp1
isgn = -isgn
ii = (3-isgn) / 2
sum(ii) = sum(ii) + 2.0_dp * dj(jp1,mp1,lp1) * temp(ii,mp1)
end do
temp2(1,jp1) = sum(1) * cbet(jp1) - sum(2) * sbet(jp1)
temp2(2,jp1) = sum(2) * cbet(jp1) + sum(1) * sbet(jp1)
end do
! Inverse B rotation and gamma rotation
do jp1 = 1, lp1
sum(1) = dj(1,jp1,lp1) * temp2(1,1)
sum(2) = 0.0_dp
isgn = 1 - 2 * mod((lp1-jp1),2)
do mp1 = 2,lp1
isgn = -isgn
ii = (3-isgn) / 2
sum(ii) = sum(ii) + 2.0_dp * dj(mp1,jp1,lp1) * temp2(ii,mp1)
end do
indx = ind + jp1
rcof(1,indx) = sum(1) * cgam(jp1) - sum(2) * sgam(jp1)
rcof(2,indx) = sum(2) * cgam(jp1) + sum(1) * sgam(jp1)
end do
ind = ind + lp1
end do
end subroutine SHRotateCoef