Commit 7c2f9e51 authored by Maximilian Schanner's avatar Maximilian Schanner
Browse files

Add secular variation and power spectrum to core.

parent 1db5d77f
......@@ -38,10 +38,12 @@ accessed, see also the `list of included models
# expose only utilities and core modules
__all__ = ['master_curve', 'dipole_series', 'file2splines', 'field',
'built_in_models', 'coefficients', 'utils', 'Model']
'built_in_models', 'coefficients', 'utils', 'secular_variation',
'power_spectrum', 'Model']
from pymagglobal.core import master_curve, dipole_series, file2splines, \
field, coefficients, built_in_models, Model
field, coefficients, secular_variation, power_spectrum, built_in_models, \
Model
from pymagglobal import utils
__version__ = '0.1.2'
......@@ -296,6 +296,106 @@ def coefficients(epoch, splines, cov_splines=None, check=True):
return ls, ms, coeffs, np.sqrt(covs)
def secular_variation(epoch, splines, cov_splines=None, check=True):
'''Evaluate splines at an epoch and return the secular variation
coefficients, together with the appropriate degrees and orders.
Parameters
----------
epoch : float or array
The epoch(s) at which to return the coefficients, in years.
splines : scipy.interpolate.BSpline or Model
An instance of Model or splines specifying the model.
check : bool, optional
If a Model is given, check the input validity.
Returns
-------
list
The list of coefficient degrees
list
The list of coefficient orders
ndarray
Secular variation coefficients of the model for the given epoch,
corresponding to the lists of degrees and orders.
'''
if isinstance(splines, Model):
if check:
for it in np.atleast_1d(epoch):
splines.valid_epoch(epoch)
splines = splines.splines
coeffs = splines.derivative()(epoch)
# set up degrees and orders for output purposes
ls = [i2lm_l(it) for it in range(splines.c.shape[1])]
ms = [i2lm_m(it) for it in range(splines.c.shape[1])]
if cov_splines is None:
return ls, ms, coeffs
def power_spectrum(epoch, splines, check=True):
'''Evaluate splines at an epoch and calculate the spatial power spectrum.
Parameters
----------
epoch : float or array
The epoch(s) at which to return the coefficients, in years.
splines : scipy.interpolate.BSpline or Model
An instance of Model or splines specifying the model.
check : bool, optional
If a Model is given, check the input validity.
Returns
-------
list
The list of coefficient degrees
array
The power spectrum over the given coefficients, with one power per
degree.
'''
if isinstance(splines, Model):
if check:
for it in np.atleast_1d(epoch):
splines.valid_epoch(epoch)
splines = splines.splines
ls = [i2lm_l(it) for it in range(splines.c.shape[1])]
coeffs = splines(epoch)
return np.unique(ls), _power_spectrum(ls, coeffs)
def _power_spectrum(ls, mu_gs, sd_gs=0):
'''For given coefficients, calculate the power spectrum.
Parameters
----------
ls : array-like
A list of degrees
mu_gs : array-like
Gauss coefficients for which the power spectrum will be calculated.
This can also be the mean of a multivariate Gaussian over the
coefficients.
sd_gs : array-like, optional
The standard deviations of the coefficients. If reported, the power
spectrum is considered a random variable.
Returns
-------
array
The power spectrum over the given coefficients, with one power per
degree.
'''
lmax = np.max(ls)
ps = np.full((lmax, ) + mu_gs.shape[1:], np.nan)
gs = (mu_gs**2 + sd_gs**2)
for ell in range(1, lmax + 1):
mask = np.equal(ls, ell)
ps[ell-1, ...] = (ell+1)*gs[mask].sum(axis=0)
return ps
def dipole_series(times, splines, cov_splines=None, check=True):
'''Create a dipole-moment time series from a splines object.
......
Markdown is supported
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