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

Implement tests against new version of orthopoly.

parent 1c375732
......@@ -63,7 +63,7 @@ setuptools.setup(
'matplotlib>=2.2.5',
'cartopy>=0.17'
],
extras_require={'tests': ['orthopoly>=0.7', 'packaging']},
extras_require={'tests': ['orthopoly>=0.8', 'packaging']},
package_data={'pymagglobal': ['dat/*']},
include_package_data=True,
entry_points={'console_scripts': [f'pymagglobal = '
......
......@@ -21,18 +21,13 @@ import unittest
import numpy as np
from pymagglobal.utils import lm2i, i2lm_l, i2lm_m, lmax2N, REARTH, equi_sph, \
dsh_basis
dsh_basis, get_z_at
def get_test_points(lmax=1, n=100):
tps = equi_sph(n)
# set upt the points as expected by pyfield
z_at = np.empty((4, tps.shape[1]), order='F')
z_at[0] = np.rad2deg(tps[0])
z_at = get_z_at(n, t=2020)
# BUGFIX: white edges in cartopy
z_at[1] = np.rad2deg(tps[1]) - 180
z_at[2] = REARTH
z_at[3] = 2020
z_at[1] -= 180
spatial = np.empty((lmax2N(lmax), 3*z_at.shape[1]), order='F')
return z_at, spatial
......@@ -126,12 +121,6 @@ class TestBasisFunction(unittest.TestCase):
try:
from orthopoly import spherical_harmonic as sh_ref
def translate_to_ortho(z_at):
'''Translate input points to the orthopoly coordinates.'''
theta = z_at[0]
phi = (z_at[1] + 360 * (z_at[1] < 0))
return np.deg2rad([theta, phi])
def norm(l):
return np.sqrt((2*l+1) / 4 / np.pi)
......@@ -142,15 +131,19 @@ try:
n = 1000
lmax = 3
z, spatial = get_test_points(lmax=lmax, n=n)
# clip for orthopoly
z[1] = np.clip(z[1], -180, 180)
dsh_basis(lmax, z, spatial)
t, p = translate_to_ortho(z)
t, p = sh_ref.latlon2sph(90-z[0], z[1])
for it in range(lmax2N(lmax)):
north, east = sh_ref.grad_sph_har(t, p, i2lm_l(it),
-i2lm_m(it))
self.assertTrue(np.allclose(spatial[it, 0::3],
north / norm(i2lm_l(it))))
self.assertTrue(np.allclose(spatial[it, 0::3],
north / norm(i2lm_l(it))))
ell = i2lm_l(it)
m = i2lm_m(it)
north, east = sh_ref.grad_sph_har(t, p, ell, m)
self.assertTrue(np.allclose((-1)**m * spatial[it, 0::3],
north / norm(ell)))
# different sign due to convention in magnetics
self.assertTrue(np.allclose(-(-1)**m * spatial[it, 1::3],
east / norm(ell)))
except ImportError:
......
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