insolation

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

This code was originally inspired by 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 Berger [1978]. Further references: Berger and Loutre [1991].

climlab.solar.insolation.annual_insolation(lat, orb={'ecc': 0.017236, 'long_peri': 281.37, 'obliquity': 23.446}, S0=1365.2, days_per_year=365.2422)

Compute annual average insolation given latitude 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

Parameters:

• lat (array) – Latitude in degrees (-90 to 90).

• orb (dict) –

  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

• S0 (float) – solar constant n - unit: \(\textrm{W}/\textrm{m}^2\) n - default value: 1365.2

• days_per_year (float) – number of days in a year (optional) (default: 365.2422) Reads the length of the year from constants if available.

Returns:

Annual average solar radiation in unit \(\textrm{W}/\textrm{m}^2\).

Dimensions of output are (lat.size, ecc.size)

Return type:

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 lookup_parameters()

For more information about computation of solar insolation see the Tutorials chapter.

Note

Annual mean insolation is computed numerically by evenly sampling the instant insolation over one calendar year.

climlab.solar.insolation.coszen_daily_insolation_weighted(phi, delta)

Cosine of solar zenith angle, insolation-weighted daily average.

Inputs:

• phi: latitude in radians

• delta: solar declination angle in radians

climlab.solar.insolation.coszen_daily_time_weighted(phi, delta)

Cosine of solar zenith angle averaged in time over 24 hours.

Inputs:

• phi: latitude in radians

• delta: solar declination angle in radians

climlab.solar.insolation.coszen_daily_time_weighted_sunlit(delta, phi)

Cosine of solar zenith angle averaged in time sunlit hours only.

Inputs:

• phi: latitude in radians

• delta: solar declination angle in radians

climlab.solar.insolation.coszen_instantaneous(phi, delta, h)

Cosine of solar zenith angle (instantaneous).

Returns zero if the sun is below the horizon.

Inputs:

• phi: latitude in radians

• delta: solar declination angle in radians

• h: hour angle in radians

climlab.solar.insolation.daily_insolation(lat, day, orb={'ecc': 0.017236, 'long_peri': 281.37, 'obliquity': 23.446}, S0=1365.2, day_type=1, days_per_year=365.2422)

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

Parameters:

• lat (array) – Latitude in degrees (-90 to 90).

• day (array) – Indicator of time of year. See argument day_type for details about format.

• orb (dict) –

  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

• S0 (float) – solar constant n - unit: \(textrm{W}/\textrm{m}^2\) n - default value: 1365.2

• day_type (int) –

  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:

ValueError if day_type is neither 1 nor 2

Returns:

Daily average solar radiation in unit \(\textrm{W}/\textrm{m}^2\).

Dimensions of output are (lat.size, day.size, ecc.size)

Return type:

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 lookup_parameters()

For more information about computation of solar insolation see the Tutorials chapter.

Note

Calling insolation = daily_insolation(lat, day, orb, S0) is equivalent to:

coszen, irradiance_factor = daily_insolation_factors(lat, day, orb)
insolation = S0 * irradiance_factor * coszen

Computing the zenith angle with daily_insolation_factors allows for optional time averaging choices which may be important for certain radiative transfer calculations that are sensitive to zenith angle.

climlab.solar.insolation.daily_insolation_factors(lat, day, orb={'ecc': 0.017236, 'long_peri': 281.37, 'obliquity': 23.446}, day_type=1, days_per_year=365.2422, weighting='time')

Compute daily average cosine of solar zenith angle and irradiance factor given latitude, time of year and orbital parameters. Multiple zenith angle averaging methods are supported.

Input arguments:

• lat: Latitude in degrees (-90 to 90)

• day: Indicator of time of year. See docs for daily_insolation().

• orb: Orbital parameter dictionary. See docs for daily_insolation().

• day_type: Convention for specifying time of year. See docs for daily_insolation().

• days_per_year: length of calendar year in days (default = 365.2422)

• weighting: flag to specify averaging method for solar zenith angle. Valid options are

  • 'time' (default): unweighted 24 hour daily average

  • 'sunlit': time average over sunlit hours

  • 'insolation': insolation-weighted average

This function retuns two values:

• coszen: the cosine of the average zenith angle (dimensionless), averaged according to the weighting argument

• irradiance_factor: ratio of current total irradiance to its annual average (dimensionless)

For weighting='sunlit' and weighting='insolation', the irradiance factor is reduced commensurate with the increased coszen value to ensure that the product coszen * irrandiance_factor always gives the 24 hour time-weighted daily average irradiance relative to the solar constant.

In all cases, daily average insolation can then be computed from S0 * coszen * irradiance_factor where S0 is the solar constant.

climlab.solar.insolation.dates_to_day_index(datetime)

Convert dates and time (assumed to be UTC) to a fractional number of days from a hypothetical January 1 00h that is exactly 80 days before the spring equinox.

Input datetime can be any of - a single datetime object - a scalar np.datetime64 object - numpy array of np.datetime64 objects - xarray.DataArray of np.datetime64 objects

Output is array of floats (either numpy or Xarray).

We actually measure number of days from a specific equinox in 2025, so the result is not bounded between 1 and 365.25. When the result is passed to the insolation calculators, it is normalized by the length of year and used as argument to trig functions, so extra multiples of 365.25 don’t matter.

climlab.solar.insolation.declination_angle(obliquity, lambda_long)

Compute solar declination angle in radians.

Inputs:

• obliquity: obliquity angle in radians

• lambda_long: solar longitude angle in radians

climlab.solar.insolation.hour_angle_at_sunset(phi, delta)

Compute the hour angle (in radians) at sunset.

Formulas based on Berger [1978] eqns (8), (9).

Inputs:

• phi: latitude in radians

• delta: solar declination angle in radians

climlab.solar.insolation.instant_insolation(lat, day, lon=0.0, orb={'ecc': 0.017236, 'long_peri': 281.37, 'obliquity': 23.446}, S0=1365.2, day_type=1, days_per_year=365.2422)

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

Parameters:

• lat (array) – Latitude in degrees (-90 to 90).

• day (array) – 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.

• lon (array) – Longitude in degrees (0 to 360), optional. Defaults to zero.

• orb (dict) –

  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

• S0 (float) – solar constant n - unit: \(\textrm{W}/\textrm{m}^2\) n - default value: 1365.2

• days_per_year (float) – number of days in a year (optional) (default: 365.2422) Reads the length of the year from constants if available.

Returns:

Daily average solar radiation in unit \(\textrm{W}/\textrm{m}^2\).

Dimensions of output are (lat.size, day.size, ecc.size)

Return type:

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 lookup_parameters()

For more information about computation of solar insolation see the Tutorials chapter.

climlab.solar.insolation.instant_insolation_factors(lat, day, lon=0.0, orb={'ecc': 0.017236, 'long_peri': 281.37, 'obliquity': 23.446}, day_type=1, days_per_year=365.2422)

Compute instantaneous cosine of solar zenith angle and irradiance factor given latitude, longitude, time of year and orbital parameters.

Input arguments:

• lat: Latitude in degrees (-90 to 90)

• day: Indicator of time of year. See docs for daily_insolation().

• lon: Longitude in degrees (0 to 360) (default = 0.)

• orb: Orbital parameter dictionary. See docs for daily_insolation().

• day_type: Convention for specifying time of year. See docs for daily_insolation().

• days_per_year: length of calendar year in days (default = 365.2422)

Same options and call signature as instant_insolation but rather than returning the insolation, this function retuns two values:

• coszen: the cosine of the instantaneous zenith angle (dimensionless)

• irradiance_factor: ratio of current total irradiance to its annual average (dimensionless)

Insolation can then be computed from S0 * coszen * irradiance_factor where S0 is the solar constant (annual average total irradiance).

climlab.solar.insolation.solar_longitude(day, orb={'ecc': 0.017236, 'long_peri': 281.37, 'obliquity': 23.446}, days_per_year=365.2422)

Estimates solar longitude from calendar day.

Method is using an approximation from Berger [1978] section 3 (lambda = 0 at spring equinox).

Function-call arguments n

Parameters:

• day (array) – Indicator of time of year.

• orb (dict) –

  a dictionary with three members (as provided by 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

• days_per_year (float) – number of days in a year (optional) (default: 365.2422) Reads the length of the year from constants if available.

Returns:

solar longitude lambda_long in degrees in dimension``( day.size, ecc.size )``

Return type:

array

Works for both scalar and vector orbital parameters.