Repository: climlab/climlab

Issue:
Inconsistent units for solar longitude
The function `climlab.solar.insolation.solar_longitude()` returns a value in radians.

But if you want to use these values to compute insolation, you need to call the function `climlab.solar.insolation.daily_insolation(..., day_type=2)`
But the code then applies `np.deg2rad()` on the input array because it expects the solar longitude values in degrees.

Discussion:

Patch:
diff --git a/climlab/solar/insolation.py b/climlab/solar/insolation.py
index f2f06ba..77d888d 100644
--- a/climlab/solar/insolation.py
+++ b/climlab/solar/insolation.py
@@ -42,8 +42,12 @@ from climlab import constants as const
 from numpy import sqrt, deg2rad, rad2deg, sin, cos, tan, arcsin, arccos, pi
 import xarray as xr
 
+#  suppress warning message generated by arccos here!
+oldsettings = np.seterr(invalid='ignore')
 
-def daily_insolation(lat, day, orb=const.orb_present, S0=const.S0, day_type=1):
+
+def daily_insolation(lat, day, orb=const.orb_present, S0=const.S0, day_type=1, 
+                     days_per_year=const.days_per_year):
     """Compute daily average insolation given latitude, time of year and orbital parameters.
 
     Orbital parameters can be interpolated to any time in the last 5 Myears with
@@ -113,35 +117,18 @@ def daily_insolation(lat, day, orb=const.orb_present, S0=const.S0, day_type=1):
         :ref:`Tutorial` chapter.
 
      """
-    # Inputs can be scalar, numpy vector, or xarray.DataArray.
-    #  If numpy, convert to xarray so that it will broadcast correctly
-    lat_is_xarray = True
-    day_is_xarray = True
-    if type(lat) is np.ndarray:
-        lat_is_xarray = False
-        lat = xr.DataArray(lat, coords=[lat], dims=['lat'])
-    if type(day) is np.ndarray:
-        day_is_xarray = False
-        day = xr.DataArray(day, coords=[day], dims=['day'])
-    ecc = orb['ecc']
-    long_peri = orb['long_peri']
-    obliquity = orb['obliquity']
-    # Convert precession angle and latitude to radians
-    phi = deg2rad( lat )
-
-    # lambda_long (solar longitude) is the angular distance along Earth's orbit measured from spring equinox (21 March)
+    phi, day, ecc, long_peri, obliquity, input_is_xarray, _ignored = _standardize_inputs(lat, day, orb)
+    # lambda_long (solar longitude) is the angular distance along Earth's orbit 
+    # measured from spring equinox (21 March)
     if day_type==1: # calendar days
-        lambda_long = solar_longitude(day,orb)
+        lambda_long = deg2rad(solar_longitude(day, orb, days_per_year))
     elif day_type==2: #solar longitude (1-360) is specified in input, no need to convert days to longitude
         lambda_long = deg2rad(day)
     else:
         raise ValueError('Invalid day_type.')
 
     # Compute declination angle of the sun
-    delta = arcsin(sin(deg2rad(obliquity)) * sin(lambda_long))
-
-    #  suppress warning message generated by arccos here!
-    oldsettings = np.seterr(invalid='ignore')
+    delta = arcsin(sin(obliquity) * sin(lambda_long))
     # Compute Ho, the hour angle at sunrise / sunset
     #  Check for no sunrise or no sunset: Berger 1978 eqn (8),(9)
     Ho = xr.where( abs(delta)-pi/2+abs(phi) < 0., # there is sunset/sunrise
@@ -151,15 +138,22 @@ def daily_insolation(lat, day, orb=const.orb_present, S0=const.S0, day_type=1):
     # this is not really the daily average cosine of the zenith angle...
     #  it's the integral from sunrise to sunset of that quantity...
     coszen = Ho*sin(phi)*sin(delta) + cos(phi)*cos(delta)*sin(Ho)
-    # Compute insolation: Berger 1978 eq (10)
-    Fsw = S0/pi*( (1+ecc*cos(lambda_long -deg2rad(long_peri)))**2 / (1-ecc**2)**2 * coszen)
-    if not (lat_is_xarray or day_is_xarray):
+    # Compute daily average insolation
+    Fsw = _compute_insolation_Berger(S0, ecc, lambda_long, long_peri, coszen)
+    if input_is_xarray:
+        return Fsw
+    else:
         # Dimensional ordering consistent with previous numpy code
         return Fsw.transpose().values
-    else:
-        return Fsw
 
-def instant_insolation(lat, day, lon=0., orb=const.orb_present, S0=const.S0):
+
+def _compute_insolation_Berger(S0, ecc, lambda_long, long_peri, coszen):
+    # Compute insolation: Berger 1978 eq (10)
+    return S0/pi*coszen*(1+ecc*cos(lambda_long - long_peri))**2 / (1-ecc**2)**2
+
+
+def instant_insolation(lat, day, lon=0., orb=const.orb_present, S0=const.S0, 
+                        days_per_year=const.days_per_year):
     """Compute instantaneous insolation given latitude, longitude, time of year and orbital parameters.
 
     Orbital parameters can be interpolated to any time in the last 5 Myears with
@@ -216,31 +210,12 @@ def instant_insolation(lat, day, lon=0., orb=const.orb_present, S0=const.S0):
         :ref:`Tutorial` chapter.
 
      """
-    # Inputs can be scalar, numpy vector, or xarray.DataArray.
-    #  If numpy, convert to xarray so that it will broadcast correctly
-    lat_is_xarray = True
-    day_is_xarray = True
-    lon_is_xarray = True
-    if type(lat) is np.ndarray:
-        lat_is_xarray = False
-        lat = xr.DataArray(lat, coords=[lat], dims=['lat'])
-    if type(day) is np.ndarray:
-        day_is_xarray = False
-        day = xr.DataArray(day, coords=[day], dims=['day'])
-    if type(lon) is np.ndarray:
-        lon_is_xarray = False
-        lon = xr.DataArray(lon, coords=[lon], dims=['lon'])        
-    ecc = orb['ecc']
-    long_peri = orb['long_peri']
-    obliquity = orb['obliquity']
-    # Convert precession angle and latitude to radians
-    phi = deg2rad( lat )
-
-    # lambda_long (solar longitude) is the angular distance along Earth's orbit measured from spring equinox (21 March)
-    lambda_long = solar_longitude(day,orb)
-
+    phi, day, ecc, long_peri, obliquity, input_is_xarray, lam = _standardize_inputs(lat, day, orb, lon)
+    # lambda_long (solar longitude) is the angular distance along Earth's orbit 
+    # measured from spring equinox (21 March)
+    lambda_long = deg2rad(solar_longitude(day, orb, days_per_year))
     # Compute declination angle of the sun
-    delta = arcsin(sin(deg2rad(obliquity)) * sin(lambda_long))
+    delta = arcsin(sin(obliquity) * sin(lambda_long))
     
     # np.fmod(day, 1.0) finds the "fractional" time of day with a range of [0,1)
     # where 0 is midnight, and 0.9999... is 23:59. lon/360 converts longitude 
@@ -251,25 +226,23 @@ def instant_insolation(lat, day, lon=0., orb=const.orb_present, S0=const.S0):
     # "fractional" units) and then convert to hour angle. The -0.5 is included
     # in order to assert that noon occurs when the sun is overhead (so h=0 at
     # LST=0.5 aka time=noon).
-    LST = np.fmod((np.fmod(day, 1.0) + (lon / 360)), 1.0)
+    LST = np.fmod((np.fmod(day, 1.0) + (lam/(2*pi))), 1.0)
     # hour angle in rad
-    h = (LST - 0.5) * 2 * np.pi
-
-    #  suppress warning message generated by arccos here!
-    oldsettings = np.seterr(invalid='ignore')
-    # instantaneous "coszen"
-    coszen = np.sin(phi)*np.sin(delta) + np.cos(phi)*np.cos(delta)*np.cos(h)
-    # Compute insolation: Berger 1978 eq (10)
-    Fsw = S0*( (1+ecc*cos(lambda_long -deg2rad(long_peri)))**2 / (1-ecc**2)**2 * coszen)
-    # assert |h|<Ho (same as checking Fsw > 0)
+    h = (LST - 0.5) * 2*pi
+    # instantaneous cosine of solar zenith angle
+    coszen = (sin(phi)*sin(delta) + cos(phi)*cos(delta)*cos(h)) * pi
+    # Compute insolation
+    Fsw = _compute_insolation_Berger(S0, ecc, lambda_long, long_peri, coszen)
+    # assert |h|<=Ho, i.e. it is day time (same as checking Fsw >= 0)
     Fsw = np.maximum(Fsw, 0.0)
-    if not (lat_is_xarray or day_is_xarray or lon_is_xarray):
+    if input_is_xarray:
+        return Fsw
+    else:
         # Dimensional ordering consistent with previous numpy code
         return Fsw.transpose().values
-    else:
-        return Fsw
 
-def solar_longitude( day, orb=const.orb_present, days_per_year = None ):
+
+def solar_longitude(day, orb=const.orb_present, days_per_year=const.days_per_year):
     """Estimates solar longitude from calendar day.
 
     Method is using an approximation from :cite:`Berger_1978` section 3
@@ -300,7 +273,7 @@ def solar_longitude( day, orb=const.orb_present, days_per_year = None ):
                                 (default: 365.2422)
                                 Reads the length of the year from
                                 :mod:`~climlab.utils.constants` if available.
-    :returns:                   solar longitude ``lambda_long``
+    :returns:                   solar longitude ``lambda_long`` in degrees
                                 in dimension``( day.size, ecc.size )``
     :rtype:                     array
 
@@ -309,8 +282,6 @@ def solar_longitude( day, orb=const.orb_present, days_per_year = None ):
 
 
     """
-    if days_per_year is None:
-        days_per_year = const.days_per_year
     ecc = orb['ecc']
     long_peri_rad = deg2rad( orb['long_peri'])
     delta_lambda = (day - 80.) * 2*pi/days_per_year
@@ -321,4 +292,38 @@ def solar_longitude( day, orb=const.orb_present, days_per_year = None ):
     lambda_long = ( lambda_long_m + (2*ecc - (ecc**3)/4)*sin(lambda_long_m - long_peri_rad) +
         (5/4)*(ecc**2) * sin(2*(lambda_long_m - long_peri_rad)) + (13/12)*(ecc**3)
         * sin(3*(lambda_long_m - long_peri_rad)) )
-    return lambda_long
+    return rad2deg(lambda_long)
+
+
+def _standardize_inputs(lat, day, orb, lon=None):
+    # Inputs can be scalar, numpy vector, or xarray.DataArray.
+    #  If numpy, convert to xarray so that it will broadcast correctly
+    lat_is_xarray = True
+    day_is_xarray = True
+    lon_is_xarray = True
+
+    if type(lat) is np.ndarray:
+        lat_is_xarray = False
+        lat = xr.DataArray(lat, coords=[lat], dims=['lat'])
+    if type(day) is np.ndarray:
+        day_is_xarray = False
+        day = xr.DataArray(day, coords=[day], dims=['day'])
+
+    ecc = orb['ecc']
+    long_peri = deg2rad(orb['long_peri'])
+    obliquity = deg2rad(orb['obliquity'])
+    # Convert latitude (and all other angles) to radians
+    phi = deg2rad(lat)
+
+    input_is_xarray = lat_is_xarray or day_is_xarray
+
+    if lon is not None:
+        if type(lon) is np.ndarray:
+            lon_is_xarray = False
+            lon = xr.DataArray(lon, coords=[lon], dims=['lon'])
+        input_is_xarray = input_is_xarray or lon_is_xarray
+        lam = deg2rad(lon)
+        print(lam)
+        return phi, day, ecc, long_peri, obliquity, input_is_xarray, lam
+    else:
+        return phi, day, ecc, long_peri, obliquity, input_is_xarray, lon


Test patch:
diff --git a/climlab/tests/test_insolation.py b/climlab/tests/test_insolation.py
index 3189815..c5fdfaa 100644
--- a/climlab/tests/test_insolation.py
+++ b/climlab/tests/test_insolation.py
@@ -1,7 +1,7 @@
 from __future__ import division, print_function, absolute_import
 import numpy as np
 from climlab import constants as const
-from climlab.solar.insolation import daily_insolation, instant_insolation
+from climlab.solar.insolation import daily_insolation, instant_insolation, solar_longitude
 from climlab.solar.orbital import OrbitalTable
 from climlab.solar.orbital.long import OrbitalTable as LongOrbitalTable
 from climlab import EBM_seasonal
@@ -25,6 +25,14 @@ def test_daily_insolation():
                   np.sum(np.cos(np.deg2rad(lat))) )
     np.testing.assert_almost_equal(Q_area_int, 341.384184481)
 
+@pytest.mark.fast
+def test_solar_longitude():
+    lat = np.linspace(-90., 90., 500)
+    days = np.linspace(0, const.days_per_year, 365)
+    Q1 = daily_insolation(lat, days, day_type=1)
+    Q2 = daily_insolation(lat, solar_longitude(days), day_type=2)
+    np.testing.assert_allclose(Q1, Q2)
+
 @pytest.mark.fast    
 def test_instant_insolation():
     lats = np.linspace(-90., 90., 181)