diff --git a/.gitignore b/.gitignore
index 87620ac7e74efee566c6ee9d2ed7281ebafb4788..4f69e2555c6ab11f274d76580a8cc5c320806856 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1 +1,2 @@
.ipynb_checkpoints/
+solar/__pycache__/
diff --git a/solar/README.md b/solar/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..45baeaac280d6a4b47b3e08e94ccce27e4ab7845
--- /dev/null
+++ b/solar/README.md
@@ -0,0 +1,7 @@
+Taken from climlab package by Brian Rose, SUNY Albany
+
+https://github.com/climlab/climlab/blob/main/climlab/utils/constants.py
+https://github.com/climlab/climlab/blob/main/climlab/solar/insolation.py
+commit 5168ea08d294750bd09399b523b98bb914196a66
+
+A use example is given by the notebook dailmean_insolation.ipynb
diff --git a/solar/constants.py b/solar/constants.py
new file mode 100644
index 0000000000000000000000000000000000000000..bea748404e36c64c44ee05f194e6c75a799f8edb
--- /dev/null
+++ b/solar/constants.py
@@ -0,0 +1,91 @@
+# constants.py
+
+#Part of the climlab package
+#Brian Rose, University at Albany
+#brose@albany.edu
+
+"""Contains a collection of physical constants for the atmosphere and ocean.
+
+ .. literalinclude:: ../code_input_manual/constants.py
+
+"""
+
+from __future__ import division
+from numpy import pi
+
+a = 6.373E6 # Radius of Earth (m)
+Lhvap = 2.5E6 # Latent heat of vaporization (J / kg)
+Lhsub = 2.834E6 # Latent heat of sublimation (J / kg)
+Lhfus = Lhsub - Lhvap # Latent heat of fusion (J / kg)
+cp = 1004. # specific heat at constant pressure for dry air (J / kg / K)
+Rd = 287. # gas constant for dry air (J / kg / K)
+kappa = Rd / cp
+Rv = 461.5 # gas constant for water vapor (J / kg / K)
+cpv = 1875. # specific heat at constant pressure for water vapor (J / kg / K)
+eps = Rd / Rv
+Omega = 2 * pi / 24. / 3600. # Earth's rotation rate, (s**(-1))
+g = 9.8 # gravitational acceleration (m / s**2)
+kBoltzmann = 1.3806488E-23 # the Boltzmann constant (J / K)
+c_light = 2.99792458E8 # speed of light (m/s)
+hPlanck = 6.62606957E-34 # Planck's constant (J s)
+# Stefan-Boltzmann constant (W / m**2 / K**4) derived from fundamental constants
+sigma = (2*pi**5 * kBoltzmann**4) / (15 * c_light**2 * hPlanck**3)
+
+S0 = 1365.2 # solar constant (W / m**2)
+# value is consistent with Trenberth and Fasullo, Surveys of Geophysics 2012
+
+ps = 1000. # approximate surface pressure (mb or hPa)
+
+rho_w = 1000. # density of water (kg / m**3)
+cw = 4181.3 # specific heat of liquid water (J / kg / K)
+
+tempCtoK = 273.15 # 0degC in Kelvin
+tempKtoC = -tempCtoK # 0 K in degC
+mb_to_Pa = 100. # conversion factor from mb to Pa
+
+# Some useful time conversion factors
+seconds_per_minute = 60.
+minutes_per_hour = 60.
+hours_per_day = 24.
+
+# the length of the "tropical year" -- time between vernal equinoxes
+# This value is consistent with Berger (1978)
+# "Long-Term Variations of Daily Insolation and Quaternary Climatic Changes"
+days_per_year = 365.2422
+seconds_per_hour = minutes_per_hour * seconds_per_minute
+minutes_per_day = hours_per_day * minutes_per_hour
+seconds_per_day = hours_per_day * seconds_per_hour
+seconds_per_year = seconds_per_day * days_per_year
+minutes_per_year = seconds_per_year / seconds_per_minute
+hours_per_year = seconds_per_year / seconds_per_hour
+# average lenghts of months based on dividing the year into 12 equal parts
+months_per_year = 12.
+seconds_per_month = seconds_per_year / months_per_year
+minutes_per_month = minutes_per_year / months_per_year
+hours_per_month = hours_per_year / months_per_year
+days_per_month = days_per_year / months_per_year
+
+area_earth = 4 * pi * a**2
+
+# present-day orbital parameters in dictionary form
+orb_present = {'ecc': 0.017236, 'long_peri': 281.37, 'obliquity': 23.446}
+# usage should be indentical to
+ #from climlab.solar.orbital import OrbitalTable
+ #orb_present = OrbitalTable.sel(kyear=-0)
+# But avoids reading the data file unnecessarily
+
+# Molecular weights in g/mol
+molecular_weight = {'N2': 28.,
+ 'O2': 32.,
+ 'H2O': 18.01528,
+ 'CO2': 44.,
+ 'N2O': 44.,
+ 'CH4': 16.,
+ 'CFC11': 136.,
+ 'CFC12': 120.,
+ 'O3': 48.,
+ 'SO2': 64.,
+ 'SO4': 96.,
+ 'H2O2': 34.,
+ 'dry air': 28.97,
+ }
diff --git a/solar/dailymean_insolation.ipynb b/solar/dailymean_insolation.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..641db6b2e50e20a4c6216d1d9f647b3d23f18606
--- /dev/null
+++ b/solar/dailymean_insolation.ipynb
@@ -0,0 +1,159 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "id": "fb125d4b-1488-460c-88a6-d03f73f5d534",
+ "metadata": {},
+ "source": [
+ "# Example notebook to illustrate how to compute daily-mean insolation"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "d7abe6d4-7386-49b2-b8ff-2f6b87934e75",
+ "metadata": {},
+ "source": [
+ "Load routines borrowed from climlab."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "id": "200c8a05-d4d2-4226-9374-5bcd966f10e0",
+ "metadata": {
+ "tags": []
+ },
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "\n",
+ "import sys\n",
+ "sys.path.append(\"/jetfs/home/avoigt/cloudsat-calipso-heating-rates/solar/\")\n",
+ "import insolation"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "4d3ce854-2f9b-4a04-a571-6e56ab2d2da8",
+ "metadata": {
+ "tags": []
+ },
+ "source": [
+ "Daily-mean insolation for latitude 30 deg and day 10 of the year. Note that the orbital values are those of Earth's todays orbit, which are used as the default values in \"insolation.daily_insolation\"."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "id": "7a13d583-daad-4951-9f83-c5b6eed1bdad",
+ "metadata": {
+ "tags": []
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Daily-mean insolation for 30 deg N and day 10: 238.329137700421\n"
+ ]
+ }
+ ],
+ "source": [
+ "print(\"Daily-mean insolation for 30 deg N and day 10:\", insolation.daily_insolation(lat=30, day=10))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "b274e026-320c-4aef-a126-5cfd993775b0",
+ "metadata": {},
+ "source": [
+ "One can also input an array of latitudes, which is helpful for the granules of CloudSat/Calipso."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "id": "08ea035e-0dc7-45a8-a4cd-04a890fdf972",
+ "metadata": {
+ "tags": []
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "<xarray.DataArray (lat: 180)> Size: 1kB\n",
+ "array([5.30489008e+02, 5.30407307e+02, 5.30162228e+02, 5.29753847e+02,\n",
+ " 5.29182291e+02, 5.28447734e+02, 5.27550403e+02, 5.26490575e+02,\n",
+ " 5.25268576e+02, 5.23884782e+02, 5.22339620e+02, 5.20633566e+02,\n",
+ " 5.18767145e+02, 5.16740931e+02, 5.14555550e+02, 5.12211674e+02,\n",
+ " 5.09710026e+02, 5.07051375e+02, 5.04236541e+02, 5.01266390e+02,\n",
+ " 4.98141839e+02, 4.94863848e+02, 4.91465705e+02, 4.89670063e+02,\n",
+ " 4.88918455e+02, 4.88752075e+02, 4.88984763e+02, 4.89505468e+02,\n",
+ " 4.90238471e+02, 4.91128262e+02, 4.92132230e+02, 4.93216595e+02,\n",
+ " 4.94353935e+02, 4.95521570e+02, 4.96700464e+02, 4.97874443e+02,\n",
+ " 4.99029625e+02, 5.00153992e+02, 5.01237060e+02, 5.02269627e+02,\n",
+ " 5.03243569e+02, 5.04151680e+02, 5.04987544e+02, 5.05745422e+02,\n",
+ " 5.06420171e+02, 5.07007164e+02, 5.07502236e+02, 5.07901626e+02,\n",
+ " 5.08201935e+02, 5.08400092e+02, 5.08493318e+02, 5.08479098e+02,\n",
+ " 5.08355160e+02, 5.08119454e+02, 5.07770132e+02, 5.07305533e+02,\n",
+ " 5.06724168e+02, 5.06024711e+02, 5.05205984e+02, 5.04266952e+02,\n",
+ " 5.03206709e+02, 5.02024475e+02, 5.00719588e+02, 4.99291494e+02,\n",
+ " 4.97739749e+02, 4.96064008e+02, 4.94264022e+02, 4.92339635e+02,\n",
+ " 4.90290780e+02, 4.88117474e+02, 4.85819817e+02, 4.83397991e+02,\n",
+ " 4.80852251e+02, 4.78182930e+02, 4.75390431e+02, 4.72475230e+02,\n",
+ " 4.69437871e+02, 4.66278963e+02, 4.62999185e+02, 4.59599277e+02,\n",
+ "...\n",
+ " 3.62391367e+02, 3.56653858e+02, 3.50828976e+02, 3.44918784e+02,\n",
+ " 3.38925400e+02, 3.32850993e+02, 3.26697790e+02, 3.20468069e+02,\n",
+ " 3.14164170e+02, 3.07788484e+02, 3.01343466e+02, 2.94831629e+02,\n",
+ " 2.88255547e+02, 2.81617859e+02, 2.74921270e+02, 2.68168553e+02,\n",
+ " 2.61362552e+02, 2.54506188e+02, 2.47602457e+02, 2.40654440e+02,\n",
+ " 2.33665304e+02, 2.26638309e+02, 2.19576813e+02, 2.12484281e+02,\n",
+ " 2.05364288e+02, 1.98220535e+02, 1.91056853e+02, 1.83877217e+02,\n",
+ " 1.76685759e+02, 1.69486784e+02, 1.62284785e+02, 1.55084465e+02,\n",
+ " 1.47890762e+02, 1.40708869e+02, 1.33544275e+02, 1.26402794e+02,\n",
+ " 1.19290613e+02, 1.12214338e+02, 1.05181061e+02, 9.81984275e+01,\n",
+ " 9.12747272e+01, 8.44189979e+01, 7.76411565e+01, 7.09521597e+01,\n",
+ " 6.43642048e+01, 5.78909842e+01, 5.15480122e+01, 4.53530531e+01,\n",
+ " 3.93266913e+01, 3.34931122e+01, 2.78812020e+01, 2.25261591e+01,\n",
+ " 1.74719704e+01, 1.27754758e+01, 8.51368208e+00, 4.79888133e+00,\n",
+ " 1.81842721e+00, 3.22768647e-02, 0.00000000e+00, 0.00000000e+00,\n",
+ " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
+ " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
+ " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
+ " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
+ " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00])\n",
+ "Coordinates:\n",
+ " * lat (lat) float64 1kB -90.0 -88.99 -87.99 -86.98 ... 87.99 88.99 90.0\n"
+ ]
+ }
+ ],
+ "source": [
+ "lats=np.linspace(-90,90,180)\n",
+ "print(insolation.daily_insolation(lat=lats, day=10))"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "cloudsat-calipso",
+ "language": "python",
+ "name": "cloudsat-calipso"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.12.2"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/solar/insolation.py b/solar/insolation.py
new file mode 100644
index 0000000000000000000000000000000000000000..dd22a3de816bd22e6603319b3c2b45c48c305322
--- /dev/null
+++ b/solar/insolation.py
@@ -0,0 +1,329 @@
+"""This module contains general-purpose routines for computing
+daily-average incoming solar radiation at the top of the atmosphere.
+
+ :Example:
+ Compute the timeseries of insolation at 65N at summer
+ solstice over the past 5 Myears::
+
+ import numpy as np
+ from climlab.solar.orbital import OrbitalTable
+ from climlab.solar.insolation import daily_insolation
+
+ # array with specified kyears (can be plain numpy or xarray.DataArray)
+ years = np.linspace(-5000, 0, 5001)
+
+ # subset of orbital parameters for specified time
+ orb = OrbitalTable.interp(kyear=years)
+
+ # insolation values for past 5 Myears at 65N at summer solstice (day 172)
+ S65 = daily_insolation(lat=65, day=172, orb=orb)
+ # returns an xarray.DataArray object with insolation values in W/m2
+
+.. note::
+
+ Ported and modified from MATLAB code daily_insolation.m \n
+ *Original authors:* \n
+
+ Ian Eisenman and Peter Huybers, Harvard University, August 2006
+
+ Available online at http://eisenman.ucsd.edu/code/daily_insolation.m
+
+If using calendar days, solar longitude is found using an
+approximate solution to the differential equation representing conservation
+of angular momentum (Kepler's Second Law). Given the orbital parameters
+and solar longitude, daily average insolation is calculated exactly
+following :cite:`Berger_1978`. Further references: :cite:`Berger_1991`.
+
+"""
+
+from __future__ import division
+import numpy as np
+#from climlab import constants as const
+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,
+ 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
+ ``climlab.solar.orbital.OrbitalTable`` (see example above).
+
+ Longer orbital tables are available with ``climlab.solar.orbital.LongOrbitalTable``
+
+ Inputs can be scalar, ``numpy.ndarray``, or ``xarray.DataArray``.
+
+ The return value will be ``numpy.ndarray`` if **all** the inputs are ``numpy``.
+ Otherwise ``xarray.DataArray``.
+
+ **Function-call argument** \n
+
+ :param array lat: Latitude in degrees (-90 to 90).
+ :param array day: Indicator of time of year. See argument ``day_type``
+ for details about format.
+ :param dict orb: a dictionary with three members (as provided by
+ ``climlab.solar.orbital.OrbitalTable``)
+
+ * ``'ecc'`` - eccentricity
+
+ * unit: dimensionless
+ * default value: ``0.017236``
+
+ * ``'long_peri'`` - longitude of perihelion (precession angle)
+
+ * unit: degrees
+ * default value: ``281.37``
+
+ * ``'obliquity'`` - obliquity angle
+
+ * unit: degrees
+ * default value: ``23.446``
+
+ :param float S0: solar constant \n
+ - unit: :math:`\\textrm{W}/\\textrm{m}^2` \n
+ - default value: ``1365.2``
+ :param int day_type: Convention for specifying time of year (+/- 1,2) [optional].
+
+ *day_type=1* (default):
+ day input is calendar day (1-365.24), where day 1
+ is January first. The calendar is referenced to the
+ vernal equinox which always occurs at day 80.
+
+ *day_type=2:*
+ day input is solar longitude (0-360 degrees). Solar
+ longitude is the angle of the Earth's orbit measured from spring
+ equinox (21 March). Note that calendar days and solar longitude are
+ not linearly related because, by Kepler's Second Law, Earth's
+ angular velocity varies according to its distance from the sun.
+ :raises: :exc:`ValueError`
+ if day_type is neither 1 nor 2
+ :returns: Daily average solar radiation in unit
+ :math:`\\textrm{W}/\\textrm{m}^2`.
+
+ Dimensions of output are ``(lat.size, day.size, ecc.size)``
+ :rtype: array
+
+
+ Code is fully vectorized to handle array input for all arguments. \n
+ Orbital arguments should all have the same sizes.
+ This is automatic if computed from
+ :func:`~climlab.solar.orbital.OrbitalTable.lookup_parameters`
+
+ For more information about computation of solar insolation see the
+ :ref:`Tutorial` chapter.
+
+ """
+ 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 = 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(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
+ arccos(-tan(phi)*tan(delta)),
+ # otherwise figure out if it's all night or all day
+ xr.where(phi*delta>0., pi, 0.) )
+ # 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 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
+
+
+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
+ ``climlab.solar.orbital.OrbitalTable`` (see example above).
+
+ Longer orbital tables are available with ``climlab.solar.orbital.LongOrbitalTable``
+
+ Inputs can be scalar, ``numpy.ndarray``, or ``xarray.DataArray``.
+
+ The return value will be ``numpy.ndarray`` if **all** the inputs are ``numpy``.
+ Otherwise ``xarray.DataArray``.
+
+ **Function-call argument** \n
+
+ :param array lat: Latitude in degrees (-90 to 90).
+ :param array day: Indicator of time of year. Format is calendar day (1-365.24), where day 1
+ is January first. The calendar is referenced to the
+ vernal equinox which always occurs at day 80.
+ :param array lon: Longitude in degrees (0 to 360), optional. Defaults to zero.
+ :param dict orb: a dictionary with three members (as provided by
+ ``climlab.solar.orbital.OrbitalTable``)
+
+ * ``'ecc'`` - eccentricity
+
+ * unit: dimensionless
+ * default value: ``0.017236``
+
+ * ``'long_peri'`` - longitude of perihelion (precession angle)
+
+ * unit: degrees
+ * default value: ``281.37``
+
+ * ``'obliquity'`` - obliquity angle
+
+ * unit: degrees
+ * default value: ``23.446``
+
+ :param float S0: solar constant \n
+ - unit: :math:`\\textrm{W}/\\textrm{m}^2` \n
+ - default value: ``1365.2``
+ :returns: Daily average solar radiation in unit
+ :math:`\\textrm{W}/\\textrm{m}^2`.
+
+ Dimensions of output are ``(lat.size, day.size, ecc.size)``
+ :rtype: array
+
+
+ Code is fully vectorized to handle array input for all arguments. \n
+ Orbital arguments should all have the same sizes.
+ This is automatic if computed from
+ :func:`~climlab.solar.orbital.OrbitalTable.lookup_parameters`
+
+ For more information about computation of solar insolation see the
+ :ref:`Tutorial` chapter.
+
+ """
+ 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(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
+ # to time since moving along the longitude axis produces the same effect as
+ # changing time while keeping longitude the same. the fractional time and
+ # fractional longitude are added together since they now both represent
+ # longitude/time of day. This lets us calculate the local solar time (in
+ # "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) + (lam/(2*pi))), 1.0)
+ # hour angle in rad
+ 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 input_is_xarray:
+ return Fsw
+ else:
+ # Dimensional ordering consistent with previous numpy code
+ return Fsw.transpose().values
+
+
+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
+ (lambda = 0 at spring equinox).
+
+ **Function-call arguments** \n
+
+ :param array day: Indicator of time of year.
+ :param dict orb: a dictionary with three members (as provided by
+ :class:`~climlab.solar.orbital.OrbitalTable`)
+
+ * ``'ecc'`` - eccentricity
+
+ * unit: dimensionless
+ * default value: ``0.017236``
+
+ * ``'long_peri'`` - longitude of perihelion
+ (precession angle)
+
+ * unit: degrees
+ * default value: ``281.37``
+
+ * ``'obliquity'`` - obliquity angle
+
+ * unit: degrees
+ * default value: ``23.446``
+ :param float days_per_year: number of days in a year (optional)
+ (default: 365.2422)
+ Reads the length of the year from
+ :mod:`~climlab.utils.constants` if available.
+ :returns: solar longitude ``lambda_long`` in degrees
+ in dimension``( day.size, ecc.size )``
+ :rtype: array
+
+ Works for both scalar and vector orbital parameters.
+
+
+
+ """
+ ecc = orb['ecc']
+ long_peri_rad = deg2rad( orb['long_peri'])
+ delta_lambda = (day - 80.) * 2*pi/days_per_year
+ beta = sqrt(1 - ecc**2)
+ lambda_long_m = -2*((ecc/2 + (ecc**3)/8 ) * (1+beta) * sin(-long_peri_rad) -
+ (ecc**2)/4 * (1/2 + beta) * sin(-2*long_peri_rad) + (ecc**3)/8 *
+ (1/3 + beta) * sin(-3*long_peri_rad)) + delta_lambda
+ 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 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)
+ return phi, day, ecc, long_peri, obliquity, input_is_xarray, lam
+ else:
+ return phi, day, ecc, long_peri, obliquity, input_is_xarray, lon