chaosmagpy.model_utils.synth_values

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

Computes radial, colatitude and azimuthal field components from the magnetic potential field in terms of spherical harmonic coefficients.

Parameters:

coeffsndarray, shape (…, M)

Coefficients of the spherical harmonic expansion. The last dimension is equal to the number of coefficients.

radiusfloat or ndarray, shape (…)

Array containing the radius in kilometers.

thetafloat or ndarray, shape (…)

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

phifloat or ndarray, shape (…)

Array containing the longitude in degrees.

nmaxint, positive, optional

Maximum degree up to which expansion is to be used (default is given by the last dimension of coeffs, that is, M = nmax(nmax+2)). This value can also be smaller if only coefficients at low degree should contribute.

nminint, positive, optional

Minimum degree of the expansion (defaults to 1). This value must correspond to the minimum degree of the coefficients and cannot be larger.

mmaxint, positive, optional

Maximum order of the spherical harmonic expansion (defaults to nmax). This value can also be smaller. For example, if only the zonal part of the expansion should be used, set mmax = 0.

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

Magnetic field source (default is an internal source).

gridbool, optional

If True, field components are computed on a regular grid. Arrays theta and phi must have one dimension less than the output grid since the grid will be created as their outer product (defaults to False).

Returns:

B_radius, B_theta, B_phindarray, shape (…)

Radial, colatitude and azimuthal field components.

Warning

The function can also evaluate 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 returned by this function. 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 component B_theta will be along the direction \(\mathbf{e}_\theta = \mathbf{e}_x\) and B_phi along \(\mathbf{e}_\phi = \mathbf{e}_y\). However, if phi = 180. is chosen, B_theta will be along \(\mathbf{e}_\theta = -\mathbf{e}_x\) and B_phi along \(\mathbf{e}_\phi = -\mathbf{e}_y\).

Notes

The function can work with different grid shapes, but the inputs have to satisfy NumPy’s broadcasting rules (grid=False, default). This also applies to the dimension of the coefficients coeffs excluding the last dimension.

The optional parameter grid is for convenience. If set to True, a singleton dimension is appended (prepended) to theta (phi) for broadcasting to a regular grid. The other inputs radius and coeffs must then be broadcastable as before but now with the resulting regular grid.

Examples

The most straight forward computation uses a fully specified grid. For example, compute the magnetic field at \(N=50\) grid points on the surface.

import chaosmagpy.model_utils as cpm
import numpy as np

N = 50
coeffs = np.ones((3,))  # degree 1 coefficients for all points
radius = 6371.2 * np.ones((N,))  # radius of 50 points in km
phi = np.linspace(-180., 180., num=N)  # azimuth of 50 points in deg.
theta = np.linspace(0., 180., num=N)  # colatitude of 50 points in deg.

B = cpm.synth_values(coeffs, radius, theta, phi)
print([B[num].shape for num in range(3)])  # output shape (N,)

Instead of N points, compute the field on a regular \(N\times N\)-grid in azimuth and colatitude (slow).

radius_grid = 6371.2 * np.ones((N, N))
phi_grid, theta_grid = np.meshgrid(phi, theta)  # grid of shape (N, N)

B = cpm.synth_values(coeffs, radius_grid, theta_grid, phi_grid)
print([B[num].shape for num in range(3)])  # output shape (N, N)

But this is slow since some computations on the grid are repeatedly done. The preferred method is to use NumPy’s broadcasting rules (fast).

radius_grid = 6371.2  # float, () or (1,)-shaped array broadcasted to NxN
theta_grid = theta[:, None]  # append singleton: (N, 1)
phi_grid = phi[None, :]  # prepend singleton: (1, N)

B = cpm.synth_values(coeffs, radius_grid, theta_grid, phi_grid)
print([B[num].shape for num in range(3)])  # output shape (N, N)

For convenience, you can do the same by using grid=True option, which appends a singleton dimension to theta and inserts a singleton dimension at the second to last position in the shape of phi.

B = cpm.synth_values(coeffs, radius_grid, theta, phi, grid=True)
print([B[num].shape for num in range(3)])  # output shape (N, N)

Remember that grid=False (default) will result in (N,)-shaped outputs as in the first example.