/
convert.py
207 lines (190 loc) · 8.78 KB
/
convert.py
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
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
"""
Functions for converting spherical harmonic coefficients to a different
normalization conventions.
"""
import numpy as _np
import warnings as _warnings
from scipy.special import factorial as _factorial
def convert(coeffs_in, normalization_in=None, normalization_out=None,
csphase_in=None, csphase_out=None, lmax=None):
"""
Convert an array of spherical harmonic coefficients to a different
normalization convention.
Usage
-----
coeffs_out = convert(coeffs_in, [normalization_in, normalization_out,
csphase_in, csphase_out, lmax])
Returns
-------
coeffs_out : ndarray, size (2, lmax+1, lmax+1)
An array of spherical harmonic coefficients with the new
normalization convention.
Parameters
----------
coeffs_in : ndarray
The array of imput spherical harmonic coefficients.
normalization_in : str, optional, default = None
Normalization of the output coefficients: '4pi', 'ortho'
'schmidt', or 'unnorm', for geodesy 4pi normalized,
orthonormalized, Schmidt semi-normalized, or unnormalized
coefficients, respectively.
normalization_out : str, optional, default = None
Normalization of the output coefficients: '4pi', 'ortho'
'schmidt', or 'unnorm', for geodesy 4pi normalized,
orthonormalized, Schmidt semi-normalized, or unnormalized
coefficients, respectively.
csphase_in : int, optional, default = None
Condon-Shortley phase convention of the input coefficients: 1 to
exclude the phase factor, or -1 to include it.
csphase_out : int, optional, default = None
Condon-Shortley phase convention of the output coefficients: 1 to
exclude the phase factor, or -1 to include it.
lmax : int, optional, default = coeffs.shape[1] - 1
Maximum spherical harmonic degree to output. If lmax is larger than
that of the input coefficients, the output array will be zero
padded.
Notes
-----
This routine will convert an array of spherical harmonic coefficients
to a different normalization convention and different Condon-Shortley
phase convention. Optionally, a different maximum spherical harmonic
degree can be specified. If this degree is smaller than that of the
input coefficients, the input coefficients will be truncated. If this
degree is larger than the input coefficients, then the output
coefficients will be zero padded.
"""
# check argument consistency
if normalization_in is not None:
if type(normalization_in) is not str:
raise ValueError('normalization_in must be a string. ' +
'Input type was {:s}'
.format(str(type(normalization_in))))
if normalization_in.lower() not in ('4pi', 'ortho', 'schmidt',
'unnorm'):
raise ValueError(
"normalization_in must be '4pi', 'ortho', 'schmidt', or " +
"'unnorm'. Provided value was {:s}"
.format(repr(normalization_in))
)
if normalization_out is None:
raise ValueError("normalization_in and normalization_out " +
"must both be specified.")
if normalization_out is not None:
if type(normalization_out) is not str:
raise ValueError('normalization_out must be a string. ' +
'Input type was {:s}'
.format(str(type(normalization_out))))
if normalization_out.lower() not in ('4pi', 'ortho', 'schmidt',
'unnorm'):
raise ValueError(
"normalization_out must be '4pi', 'ortho', 'schmidt', or" +
" 'unnorm'. Provided value was {:s}"
.format(repr(normalization_out))
)
if normalization_in is None:
raise ValueError("normalization_in and normalization_out " +
"must both be specified.")
if csphase_in is not None:
if csphase_in != 1 and csphase_in != -1:
raise ValueError(
"csphase_in must be 1 or -1. Input value was {:s}"
.format(repr(csphase_in)))
if csphase_out is None:
raise ValueError("csphase_in and csphase_out must both be " +
"specified.")
if csphase_out is not None:
if csphase_out != 1 and csphase_out != -1:
raise ValueError(
"csphase_out must be 1 or -1. Input value was {:s}"
.format(repr(csphase_out)))
if csphase_in is None:
raise ValueError("csphase_in and csphase_out must both be " +
"specified.")
lmaxin = coeffs_in.shape[1] - 1
if lmax is None:
lmaxout = lmaxin
else:
lmaxout = lmax
lconv = min(lmaxin, lmaxout)
if ((normalization_in == 'unnorm' or normalization_out ==
'unnorm') and lconv > 85):
_warnings.warn("Conversions with unnormalized coefficients are " +
"stable only for degrees less than or equal to " +
"85. lmax of the output coefficients will be " +
"truncated after degree 85. The spherical " +
"harmonic degree of the input coefficients was " +
"{:d}.".format(lmaxin), category=RuntimeWarning)
lconv = 85
degrees = _np.arange(lconv + 1)
if _np.iscomplexobj(coeffs_in):
coeffs = _np.zeros((2, lmaxout+1, lmaxout+1), dtype=_np.complex128)
else:
coeffs = _np.zeros((2, lmaxout+1, lmaxout+1))
coeffs[:, :lconv+1, :lconv+1] = coeffs_in[:, :lconv+1, :lconv+1]
if normalization_in == normalization_out:
pass
elif normalization_in == '4pi' and normalization_out == 'schmidt':
for l in degrees:
coeffs[:, l, :l+1] *= _np.sqrt(2. * l + 1.)
elif normalization_in == '4pi' and normalization_out == 'ortho':
coeffs *= _np.sqrt(4. * _np.pi)
elif normalization_in == '4pi' and normalization_out == 'unnorm':
for l in degrees:
ms = _np.arange(l+1)
conv = (2. * l + 1.) * _factorial(l-ms) / _factorial(l+ms)
if not _np.iscomplexobj(coeffs):
conv[1:] *= 2.
coeffs[:, l, :l+1] *= _np.sqrt(conv)
elif normalization_in == 'schmidt' and normalization_out == '4pi':
for l in degrees:
coeffs[:, l, :l+1] /= _np.sqrt(2. * l + 1.)
elif normalization_in == 'schmidt' and normalization_out == 'ortho':
for l in degrees:
coeffs[:, l, :l+1] *= _np.sqrt(4. * _np.pi / (2. * l + 1.))
elif normalization_in == 'schmidt' and normalization_out == 'unnorm':
for l in degrees:
ms = _np.arange(l+1)
conv = _factorial(l-ms) / _factorial(l+ms)
if not _np.iscomplexobj(coeffs):
conv[1:] *= 2.
coeffs[:, l, :l+1] *= _np.sqrt(conv)
elif normalization_in == 'ortho' and normalization_out == '4pi':
coeffs /= _np.sqrt(4. * _np.pi)
elif normalization_in == 'ortho' and normalization_out == 'schmidt':
for l in degrees:
coeffs[:, l, :l+1] *= _np.sqrt((2. * l + 1.) / (4. * _np.pi))
elif normalization_in == 'ortho' and normalization_out == 'unnorm':
for l in degrees:
ms = _np.arange(l+1)
conv = (2. * l + 1.) * _factorial(l-ms) \
/ 4. / _np.pi / _factorial(l+ms)
if not _np.iscomplexobj(coeffs):
conv[1:] *= 2.
coeffs[:, l, :l+1] *= _np.sqrt(conv)
elif normalization_in == 'unnorm' and normalization_out == '4pi':
for l in degrees:
ms = _np.arange(l+1)
conv = _factorial(l+ms) / (2. * l + 1.) / _factorial(l-ms)
if not _np.iscomplexobj(coeffs):
conv[1:] /= 2.
coeffs[:, l, :l+1] *= _np.sqrt(conv)
elif normalization_in == 'unnorm' and normalization_out == 'schmidt':
for l in degrees:
ms = _np.arange(l+1)
conv = _factorial(l+ms) / _factorial(l-ms)
if not _np.iscomplexobj(coeffs):
conv[1:] /= 2.
coeffs[:, l, :l+1] *= _np.sqrt(conv)
elif normalization_in == 'unnorm' and normalization_out == 'ortho':
for l in degrees:
ms = _np.arange(l+1)
conv = 4. * _np.pi * _factorial(l+ms) / (2. * l + 1.) / \
_factorial(l-ms)
if not _np.iscomplexobj(coeffs):
conv[1:] /= 2.
coeffs[:, l, :l+1] *= _np.sqrt(conv)
if csphase_in != csphase_out:
for m in degrees:
if m % 2 == 1:
coeffs[:, m:lconv+1, m] = - coeffs[:, m:lconv+1, m]
return coeffs