Commit a894e973 authored by Maximilian Schanner's avatar Maximilian Schanner
Browse files

Added a convenience routine to generate input points.

parent da083c4f
......@@ -387,15 +387,8 @@ def maps(args):
# check whether the epoch is in the range of the model
epoch = args.model.valid_epoch(args.epoch)
# set up a grid of equidistant points on the sphere
tps = utils.equi_sph(args.res)
# set up the points as expected by dsh_basis
z_at = np.empty((4, tps.shape[1]), order='F')
z_at[0] = np.rad2deg(tps[0])
# BUGFIX: white edges in cartopy
z_at[1] = np.rad2deg(tps[1]) - 180
z_at[2] = utils.REARTH
z_at[3] = epoch
z_at = utils.get_z_at(args.res, t=epoch)
# use the core function to get the field
field = core.field(z_at, args.model.splines)
# convert the field if necessary
......@@ -255,7 +255,21 @@ def i2lm_m(i):
def equi_sph(n, twopi=True):
'''Regularly place points on the surface of the unit sphere [Deserno]_. the
number of output points may differ from n.
number of output points may differ from n.
n : int
The approximate number of points
twopi : bool, optional
Whether to include 2*pi as a duplicate of 0 for the longitutinal angle
(useful for plotting purposes).
array of shape (2, n')
azimutal and longitudinal angles in radians of n' points, that are
equally spaced on the sphere. n' is approximately n.
......@@ -279,6 +293,54 @@ def equi_sph(n, twopi=True):
return np.array([ts, ps], dtype=float)
def get_z_at(n, R=REARTH, t=1900., random=False):
'''Get input points on the sphere. This is a convenience routine to allow
easy construction of field maps and synthetic data from the models.
n : int
The approximate number of points. Due to the points being equally
spaced on the sphere, the actual number may be slightly higher.
R : float, optional
The radius of the sphere. By default this is the Earth's radius.
t : float, optional
The epoch at which the points are generated.
random : bool, optional
If true, exactly n random points are returned. This is useful for
generating synthetic data.
z_at : array of shape (4, n')
* z_at[0] contains co-latitudes in degrees.
* z_at[1] contains longitudes in degrees.
* z_at[2] contains radii in km.
* z_at[3] contains dates in years.
n' is approximately n.
if random:
n_ret = n
n = 100*n
angles = equi_sph(n)
z_at = np.zeros((4, angles.shape[1]), order='F')
z_at[0] = np.rad2deg(angles[0])
z_at[1] = np.rad2deg(angles[1])
z_at[2] = R
z_at[3] = t
if random:
rng = np.random.default_rng()
inds = rng.choice(z_at.shape[1], size=n_ret, replace=False)
return z_at[:, inds]
# BUGFIX: white edges in cartopy
z_at[1] -= 180
return z_at
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
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