chaosmagpy.model_utils.design_gauss

chaosmagpy.model_utils.design_gauss(radius, theta, phi, nmax, *, nmin=None, mmax=None, source=None)

Computes matrices to connect the radial, colatitude and azimuthal field components to the magnetic potential field in terms of spherical harmonic coefficients (Schmidt quasi-normalized).

Parameters:

radiusndarray, shape (…)

Array containing the radius in kilometers.

thetandarray, shape (…)

Array containing the colatitude in degrees \([0^\circ,180^\circ]\).

phindarray, shape (…)

Array containing the longitude in degrees.

nmaxint, positive

Maximum degree of the sphercial harmonic expansion.

nminint, positive, optional

Minimum degree of the sphercial harmonic expansion (defaults to 1).

mmaxint, positive, optional

Maximum order of the spherical harmonic expansion (defaults to nmax). For mmax = 0, for example, only the zonal terms are returned.

source{‘internal’, ‘external’}, optional

Magnetic field source (default is an internal source).

Returns:

A_radius, A_theta, A_phindarray, shape (…, M)

Matrices for radial, colatitude and azimuthal field components. The second dimension M varies depending on nmax, nmin and mmax.

Warning

The function can also return the design matrix for the field components at the geographic poles, i.e. where theta == 0. or theta == 180.. However, users should be careful when doing this because the vector basis for spherical geocentric coordinates, \({{\mathbf{e}_r, \mathbf{e}_\theta, \mathbf{e}_\phi}}\), depends on longitude, which is not well defined at the poles. That is, at the poles, any value for the longitude maps to the same location in euclidean coordinates but gives a different vector basis in spherical geocentric coordinates. Nonetheless, by choosing a specific value for the longitude at the poles, users can define the vector basis, which then establishes the meaning of the spherical geocentric components. The vector basis for the horizontal components is defined as

\[\begin{split}\mathbf{e}_\theta &= \cos\theta\cos\phi\mathbf{e}_x - \cos\theta\sin\phi\mathbf{e}_y - \sin\theta\mathbf{e}_z\\ \mathbf{e}_\phi &= -\sin\phi\mathbf{e}_x + \cos\phi\mathbf{e}_y\end{split}\]

Hence, at the geographic north pole as given by theta = 0. and phi = 0. (chosen by the user), the returned design matrix A_theta will be for components along the direction \(\mathbf{e}_\theta = \mathbf{e}_x\) and A_phi for components along \(\mathbf{e}_\phi = \mathbf{e}_y\). However, if phi = 180. is chosen, A_theta will be for components along \(\mathbf{e}_\theta = -\mathbf{e}_x\) and A_phi along \(\mathbf{e}_\phi = -\mathbf{e}_y\).

Examples

Create the design matrices given 4 locations on the Earth’s surface:

>>> r = 6371.2
>>> theta = np.array([90., 100., 110., 120.])
>>> phi = 0.
>>> A_radius, A_theta, A_phi = design_gauss(r, theta, phi, nmax=1)
>>> A_radius
array([[ 1.22464680e-16, 2.00000000e+00, 0.00000000e+00],
       [-3.47296355e-01, 1.96961551e+00, 0.00000000e+00],
       [-6.84040287e-01, 1.87938524e+00, 0.00000000e+00],
       [-1.00000000e+00, 1.73205081e+00, 0.00000000e+00]])

Say, only the axial dipole coefficient is non-zero, what is the magnetic field in terms of spherical geocentric components?

>>> m = np.array([-30000, 0., 0.])  # model coefficients in nT
>>> Br = A_radius @ m; Br
array([-3.67394040e-12, 1.04188907e+04, 2.05212086e+04, 3.00000000e+04])
>>> Bt = A_theta @ m; Bt
array([-30000. , -29544.23259037, -28190.77862358, -25980.76211353])
>>> Bp = A_phi @ m; Bp
array([0., 0., 0., 0.])

A more complete example is given below:

import numpy as np
from chaosmagpy.model_utils import design_gauss, synth_values

nmax = 5  # desired truncation degree, i.e. 35 model coefficients
coeffs = np.arange(35)  # example model coefficients

# example locations
N = 10
radius = 6371.2  # Earth's surface
phi = np.linspace(-180., 180., num=N)  # azimuth in degrees
theta = np.linspace(1., 179., num=N)  # colatitude in degrees

# compute design matrices to compute predictions of the
# field components from the model coefficients, each is of shape N x 35
A_radius, A_theta, A_phi = design_gauss(r, theta, phi, nmax=nmax)

# magnetic components computed from the model
Br_pred = A_radius @ coeffs
Bt_pred = A_theta @ coeffs
Bp_pred = A_phi @ coeffs

# compute the magnetic components directly from the coefficients
Br, Bt, Bp = synth_values(coeffs, radius, theta, phi)
np.linalg.norm(Br - Br_pred)  # approx 0.
np.linalg.norm(Bt - Bt_pred)  # approx 0.
np.linalg.norm(Bp - Bp_pred)  # approx 0.