Merge pull request #129 from brian-rose/insolationwiggle

Fix insolation wiggles

climlab/solar/orbital/table.py

@@ -27,6 +27,8 @@ def _get_Berger_data(verbose=True):
                           'OBL': 'obliquity', 'PREC': 'precession'})
     # add 180 degrees to long_peri (see lambda definition, Berger 1978 Appendix)
     orbit['long_peri'] += 180.
+    # apply np.unwrap to remove discontinuities in the longitude of perihelion
+    orbit['long_peri'] = np.rad2deg(xr.apply_ufunc(np.unwrap, np.deg2rad(orbit.long_peri)))
     orbit['precession'] *= -1.
 
     orbit.attrs['Description'] = 'The Berger and Loutre (1991) orbital data table'

climlab/tests/test_insolation.py

@@ -40,6 +40,18 @@ def test_orbital_parameters():
         assert orb[k].min() > orb_expected[k][0]
         assert orb[k].max() < orb_expected[k][1]
 
+@pytest.mark.fast
+def test_orbital_interpolation():
+    # check to see if we get smooth results when we interpolate
+    # orbital parameters at high temporal frequency
+    kyears = np.linspace(-11, 0, 1001)
+    orb = OrbitalTable.interp(kyear=kyears)
+    S30 = daily_insolation(lat=30, day=172, orb=orb)
+    # there should be no abrupt changes
+    # test if insolation varies by more than 0.1 W/m2 per year
+    assert S30.diff(dim='kyear').max() < 0.1
+
+
 @pytest.mark.slow
 def test_long_orbital_parameters():
     kyears = np.arange( -1000., +500.)