Skip to content
GitLab
Commit a2fcf03c authored by Maximilian Schanner's avatar Maximilian Schanner
Browse files

More suggestive names

parent 46b6164b
Hide whitespace changes
Inline Side-by-side
pymagglobal/utils.py
View file @ a2fcf03c
......@@ -227,9 +227,9 @@ def equi_sph(n, twopi=True):
return np.array([ts, ps], dtype=float)
def dsh_basis(lmax, at, B, R=REARTH):
'''Write the magnetic field basis functions, evaluated at points at, into
the array B. These are basically the derivatives of the spherical
def dsh_basis(lmax, z, out, R=REARTH):
'''Write the magnetic field basis functions, evaluated at points z, into
the array out. These are basically the derivatives of the spherical
harmonics in spherical coordinates with some scaling factors applied. This
implementation is based on a specific recursion of the Legendre
polynomials, to guarantee sane behavior at the poles. See [Du]_ for further
......@@ -245,17 +245,17 @@ def dsh_basis(lmax, at, B, R=REARTH):
The maximal spherical harmonics degree.
at : array of shape (3, n)
The points at which to evaluate the basis functions, given as
* at[0, :] contains the co-latitude in degrees.
* at[1, :] contains the longitude in degrees.
* at[2, :] contains the radius in kilometers.
* z[0, :] contains the co-latitude in degrees.
* z[1, :] contains the longitude in degrees.
* z[2, :] contains the radius in kilometers.
B : array of shape (lmax2N(lmax), 2*n)
out : array of shape (lmax2N(lmax), 2*n)
The output array in which the basis functions are stored, as follows:
* B[i, 0::3] contains the North component of the basis
corresponding to degree l and order m,
where i = lm2i(l, m).
* B[i, 1::3] contains the East component of the basis.
* B[i, 2::3] contains the Down component of the basis.
* out[i, 0::3] contains the North component of the basis
corresponding to degree l and order m,
where i = lm2i(l, m).
* out[i, 1::3] contains the East component of the basis.
* out[i, 2::3] contains the Down component of the basis.
References
----------
......@@ -269,21 +269,21 @@ def dsh_basis(lmax, at, B, R=REARTH):
_intcheck(lmax)
# check input consitency
N = lmax2N(lmax)
if B.shape[0] != N:
if out.shape[0] != N:
raise ValueError(f"Inconsistent input. Number of degrees is {N}, but "
f"the given array is of shape {B.shape}.")
f"the given array is of shape {out.shape}.")
try:
if B.shape[1] != 3*at.shape[1]:
raise ValueError(f"Inconsistent input. {at.shape[1]} points "
if out.shape[1] != 3*z.shape[1]:
raise ValueError(f"Inconsistent input. {z.shape[1]} points "
f"given, but the given field-array is of shape "
f"{B.shape}.")
f"{out.shape}.")
except IndexError:
raise IndexError(f"Both input arrays have to be two dimensional. "
f"Inputs are of shape {at.shape} and {B.shape}.")
f"Inputs are of shape {z.shape} and {out.shape}.")
# XXX at is in degrees!!!
cos_t = np.cos(np.deg2rad(at[0, :]))
B[:, :] = 0.
# XXX z is in degrees!!!
cos_t = np.cos(np.deg2rad(z[0, :]))
out[:, :] = 0.
Plm_all = np.array([lpmn(lmax, lmax, x)[0] for x in cos_t])
# Degree l=0 is needed by one of the recurrence formulas
for l in range(0, lmax + 1):
......@@ -300,23 +300,23 @@ def dsh_basis(lmax, at, B, R=REARTH):
# The pre-factors of the recurrence formula and the scaling for
# negative orders cancel out. The plus sign accounts for the
# Condon-Shortley phase
B[i, 0::3] += Plm/2
out[i, 0::3] += Plm/2
# P_l^{m-1} term: m -> m+1
if (0 < l and m < l): # Skip l=0; m=0, ..., l-1
i = lm2i(l, m+1)
# Zero and positive orders
B[i, 0::3] -= (l+m+1)*(l-m)*Plm/2
out[i, 0::3] -= (l+m+1)*(l-m)*Plm/2
# Negative orders
if (m+1 != 0): # Do not visit m=0 twice
B[i+1, 0::3] -= (l+m+1)*(l-m)*Plm/2
out[i+1, 0::3] -= (l+m+1)*(l-m)*Plm/2
# P_l^{m+1} term: m -> m-1
if (0 < l and 0 < m): # Skip l=0; m=1, ..., l
i = lm2i(l, m-1)
# Zero and positive orders
B[i, 0::3] += Plm/2
out[i, 0::3] += Plm/2
# Negative orders
if (m-1 != 0): # Do not visit m=0 twice
B[i+1, 0::3] += Plm/2
out[i+1, 0::3] += Plm/2
# East Component
# The following recurrence formulas is used:
......@@ -326,43 +326,43 @@ def dsh_basis(lmax, at, B, R=REARTH):
# P_{l-1}^{m+1} term: l -> l+1 and m -> m-1
if (l < lmax and 1 < m): # Skip l=lmax; m=2, ..., l
i = lm2i(l+1, m-1)
B[i, 1::3] -= Plm/2/(m-1)
B[i+1, 1::3] -= Plm/2/(m-1)
out[i, 1::3] -= Plm/2/(m-1)
out[i+1, 1::3] -= Plm/2/(m-1)
# P_{l-1}^{m-1} term: l -> l+1 and m -> m+1
if (l < lmax): # Skip l=lmax; m=0, ..., l
i = lm2i(l+1, m+1)
B[i, 1::3] -= (l+m+1)*(l+m+2)*Plm/2/(m+1)
B[i+1, 1::3] -= (l+m+1)*(l+m+2)*Plm/2/(m+1)
out[i, 1::3] -= (l+m+1)*(l+m+2)*Plm/2/(m+1)
out[i+1, 1::3] -= (l+m+1)*(l+m+2)*Plm/2/(m+1)
# Down Component
if (0 < l): # Skip l=0; m=0, ..., l
i = lm2i(l, m)
# Zero and positive orders
B[i, 2::3] = Plm
out[i, 2::3] = Plm
# Negative orders
if m != 0: # Do not visit m=0 twice
B[i+1, 2::3] = Plm
out[i+1, 2::3] = Plm
# Arrays of degrees and orders
m = np.fromiter(map(i2lm_m, range(N)), int, N).reshape(-1, 1)
ell = np.fromiter(map(i2lm_l, range(N)), int, N).reshape(-1, 1)
# sqrt 2 to account for real form
B *= np.sqrt(2 - (m == 0))
out *= np.sqrt(2 - (m == 0))
# Condon-Shortley Phase
B *= np.where(m % 2, -1, 1)
out *= np.where(m % 2, -1, 1)
# Schmidt semi-norm
B *= np.exp((gammaln(ell-np.abs(m)+1) - gammaln(ell+np.abs(m)+1))/2)
out *= np.exp((gammaln(ell-np.abs(m)+1) - gammaln(ell+np.abs(m)+1))/2)
# North component, 1:2 instead of 1 for easier broadcasting
B[:, 0::3] *= np.cos((m < 0)*np.pi/2 - np.abs(m)*np.deg2rad(at[1:2, :]))
out[:, 0::3] *= np.cos((m < 0)*np.pi/2 - np.abs(m)*np.deg2rad(z[1:2, :]))
# East Component
B[:, 1::3] *= np.sin((m < 0)*np.pi/2 - np.abs(m)*np.deg2rad(at[1:2, :]))
B[:, 1::3] *= -np.abs(m)
out[:, 1::3] *= np.sin((m < 0)*np.pi/2 - np.abs(m)*np.deg2rad(z[1:2, :]))
out[:, 1::3] *= -np.abs(m)
# Down component
B[:, 2::3] *= np.cos((m < 0)*np.pi/2 - np.abs(m)*np.deg2rad(at[1:2, :]))
B[:, 2::3] *= -(ell+1)
out[:, 2::3] *= np.cos((m < 0)*np.pi/2 - np.abs(m)*np.deg2rad(z[1:2, :]))
out[:, 2::3] *= -(ell+1)
# Scaling for sources of internal origin
B *= np.repeat((R/at[2:3, :])**(l+2), 3, axis=1)
out *= np.repeat((R/z[2:3, :])**(l+2), 3, axis=1)
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment