Test suite comparing IDL and Python results for each Geopack field model

pyspedas/geopack/tests/geopack_idl_validation_tests.py

@@ -0,0 +1,185 @@
+import unittest
+
+import numpy as np
+from geopack import geopack
+
+from pytplot import data_exists,del_data,cdf_to_tplot,tplot,subtract,tlimit
+import pyspedas
+from pyspedas import time_double
+from pyspedas.geopack import tt89
+from pyspedas.geopack import tt96
+from pyspedas.geopack import tt01
+from pyspedas.geopack import tts04
+from pyspedas.geopack.get_tsy_params import get_tsy_params
+from pyspedas.geopack.get_w_params import get_w
+from pyspedas import tinterpol
+from pytplot import join_vec, store_data, get_data
+
+from numpy.testing import assert_array_almost_equal_nulp, assert_array_max_ulp, assert_allclose
+
+
+trange = ['2015-10-16', '2015-10-17']
+
+
+def get_params(model, g_variables=None):
+    support_trange = [time_double(trange[0])-60*60*24, 
+                      time_double(trange[1])+60*60*24]
+    pyspedas.kyoto.dst(trange=support_trange)
+    pyspedas.omni.data(trange=trange)
+    join_vec(['BX_GSE', 'BY_GSM', 'BZ_GSM'])
+    if model == 't01' and g_variables is None:
+        g_variables = [6.0, 10.0]
+    else:
+        if g_variables is not None:
+            if not isinstance(g_variables, str) and not isinstance(g_variables, np.ndarray):
+                g_variables = None
+    return get_tsy_params('kyoto_dst',
+                    'BX_GSE-BY_GSM-BZ_GSM_joined',
+                    'proton_density',
+                    'flow_speed',
+                    model,
+                    pressure_tvar='Pressure',
+                    g_variables=g_variables,
+                    speed=True)
+
+
+class LoadGeopackIdlValidationTestCases(unittest.TestCase):
+    @classmethod
+    def setUpClass(cls):
+        """
+        IDL Data has to be downloaded to perform these tests
+        The SPEDAS script that creates the file: projects/themis/state/cotrans/thm_cotrans_validate.pro
+        """
+        from pyspedas.utilities.download import download
+        from pyspedas.themis.config import CONFIG
+
+        # Testing time range
+        cls.t = ['2008-03-23', '2008-03-28']
+
+        # Testing tolerance
+        cls.tol = 1e-10
+
+        # Download tplot files
+        remote_server = 'https://spedas.org/'
+        remote_name = 'testfiles/geopack_idl_validate.cdf'
+
+        datafile = download(remote_file=remote_name,
+                            remote_path=remote_server,
+                            local_path=CONFIG['local_data_dir'],
+                            no_download=False)
+        if not datafile:
+            # Skip tests
+            raise unittest.SkipTest("Cannot download data validation file")
+
+        del_data('*')
+        filename = datafile[0]
+        cdf_to_tplot(filename)
+        # Input parameters
+        cls.iopt=3
+        cls.kp=2.0
+        cls.pdyn = 2.0
+        cls.dsti = -30.0
+        cls.yimf = 0.0
+        cls.zimf = -5.0
+        cls.g1 = 6.0
+        cls.g2 = 10.0
+        cls.w1 = 8.0
+        cls.w2 = 5.0
+        cls.w3 = 9.5
+        cls.w4 = 30.0
+        cls.w5 = 18.5
+        cls.w6 = 60.0
+
+    def test_tilt(self):
+        tilt_times,tilt_vals = get_data('bt89_tilt')
+        py_tilt_vals=np.zeros(len(tilt_times))
+        for idx,t in enumerate(tilt_times):
+            py_tilt_vals[idx] = np.degrees(geopack.recalc(t))
+        store_data('py_tilt',data={'x':tilt_times,'y':py_tilt_vals})
+        idl_tilt = get_data('bt89_tilt')
+        py_tilt = get_data('py_tilt')
+        assert_allclose(idl_tilt.y, py_tilt.y, atol=.0003)
+
+    def test_igrf(self):
+        tt89('tha_state_pos_gsm', igrf_only=True)
+        py_b = get_data('tha_state_pos_gsm_bt89')
+        idl_b = get_data('bt89_igrf')
+        assert_allclose(py_b.y, idl_b.y, rtol = .01, atol=1.0)
+
+    def test_tt89(self):
+        tt89('tha_state_pos_gsm')
+        py_b = get_data('tha_state_pos_gsm_bt89')
+        idl_b = get_data('bt89')
+        assert_allclose(py_b.y, idl_b.y, atol=0.5)
+
+    def test_tt96(self):
+        tv_pos=get_data('tha_state_pos_gsm')
+        t1=tv_pos.times[0]
+        t2=tv_pos.times[-1]
+        params = np.zeros([2,10])
+        params[:,0] = self.pdyn
+        params[:,1] = self.dsti
+        params[:,2] = self.yimf
+        params[:,3] = self.zimf
+        store_data('parmod',data={'x':[t1,t2],'y':params})
+        tinterpol('parmod','tha_state_pos_gsm',method='nearest',newname='parmod_interp')
+        tplot(['parmod','parmod_interp'])
+        tt96('tha_state_pos_gsm', parmod='parmod_interp')
+        py_b = get_data('tha_state_pos_gsm_bt96')
+        idl_b = get_data('bt96')
+        subtract('bt96','tha_state_pos_gsm_bt96','bt96_diff')
+        tplot(['bt96','tha_state_pos_gsm_bt96','bt96_diff'])
+        tlimit(['2007-03-23/15:00','2007-03-23/17:00'])
+        tplot(['bt96','tha_state_pos_gsm_bt96','bt96_diff'])
+        tlimit(full=True)
+        assert_allclose(py_b.y, idl_b.y, atol=0.5)
+
+    def test_tt01(self):
+        tv_pos=get_data('tha_state_pos_gsm')
+        t1=tv_pos.times[0]
+        t2=tv_pos.times[-1]
+        params = np.zeros([2,10])
+        params[:,0] = self.pdyn
+        params[:,1] = self.dsti
+        params[:,2] = self.yimf
+        params[:,3] = self.zimf
+        params[:,4] = self.g1
+        params[:,5] = self.g2
+        store_data('parmod',data={'x':[t1,t2],'y':params})
+        tinterpol('parmod', 'tha_state_pos_gsm', method='nearest',newname='parmod_interp')
+        tt01('tha_state_pos_gsm', parmod='parmod_interp')
+        py_b = get_data('tha_state_pos_gsm_bt01')
+        idl_b = get_data('bt01')
+        subtract('bt01','tha_state_pos_gsm_bt01','bt01_diff')
+        tplot(['bt01','tha_state_pos_gsm_bt01','bt01_diff'])
+        assert_allclose(py_b.y, idl_b.y, atol=0.5)
+
+
+    def test_tts04(self):
+        tv_pos=get_data('tha_state_pos_gsm')
+        t1=tv_pos.times[0]
+        t2=tv_pos.times[-1]
+        params = np.zeros([2,10])
+        params[:,0] = self.pdyn
+        params[:,1] = self.dsti
+        params[:,2] = self.yimf
+        params[:,3] = self.zimf
+        params[:,4] = self.w1
+        params[:,5] = self.w2
+        params[:,6] = self.w3
+        params[:,7] = self.w4
+        params[:,8] = self.w5
+        params[:,9] = self.w6
+        store_data('parmod',data={'x':[t1,t2],'y':params})
+        tinterpol('parmod', 'tha_state_pos_gsm', method='nearest',newname='parmod_interp')
+        tts04('tha_state_pos_gsm', parmod='parmod_interp')
+        py_b = get_data('tha_state_pos_gsm_bts04')
+        idl_b = get_data('bts04')
+        subtract('bts04','tha_state_pos_gsm_bts04','bts04_diff')
+        tplot(['bts04','tha_state_pos_gsm_bts04','bts04_diff'])
+        assert_allclose(py_b.y, idl_b.y, atol=0.5)
+
+
+
+if __name__ == '__main__':
+    unittest.main()