Return solar time for given longitude and UTC time. This function uses the low accuracy equations presented in General Solar Position Calculations (https://www.esrl.noaa.gov/gmd/grad/solcalc/solareqns.PDF, valid 2018-11-06) published by Global Radiation Group, NOAA Global Monitoring Division.

calc_solar_time.py

"""Copyright (C) 2018 Mikko Pitkänen.

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
"""
import numpy as np
import pandas as pd
from calendar import isleap


def calc_solar_time(timeUTC, longitude):
    """Return solar time for given longitude and UTC time.

    This function is an implementation of low accuracy equations presented in
    the document General Solar Position Calculations
    (https://www.esrl.noaa.gov/gmd/grad/solcalc/solareqns.PDF, valid
    2018-11-06) published by Global Radiation Group, NOAA Global Monitoring
    Division.

    Parameters:
    timeUTC     pandas DatetimeIndex
    longitude   float, east positive

    Returns:
    np.array with dtype='datetime64[s]' and same size as timeUTC

    Example:
    true_solar_time = calc_solar_time(
        pd.DatetimeIndex([pd.Timestamp('2018-09-09 12:00:00')]), 100.0)

    See also: equation_of_time_spencer71 and equation_of_time_pvcdrom at
    https://github.com/pvlib/pvlib-python/blob/master/pvlib/solarposition.py
    (valid 2018-11-07), that provide further literature sources for these
    equations and a correction for a constant in eqtime (0.000075 should be
    0.0000075 as in this file).

    """
    # check input
    if not -180.0 <= longitude <= 180:
        raise ValueError('longitude out of range [-180, 180]')

    # define logical indexing with a True value for leap year date
    isleap_vector = np.vectorize(isleap)
    is_leap_year = isleap_vector(timeUTC.year)

    # define fraction of year y
    y = np.zeros(len(is_leap_year))
    y[is_leap_year] = \
        2 * np.pi / 366.0 * (timeUTC.dayofyear[is_leap_year] - 1 +
                             (timeUTC.hour[is_leap_year] - 12) / 24)
    y[~is_leap_year] = \
        2 * np.pi / 365.0 * (timeUTC.dayofyear[~is_leap_year] - 1 +
                             (timeUTC.hour[~is_leap_year] - 12) / 24)

    # equation of time in minutes. here 0.000075 was corrected to 0.0000075
    # as stated in the docstring
    eqtime = 229.18 * (
        0.0000075 +
        0.001868 * np.cos(y) -
        0.032077 * np.sin(y) -
        0.014615 * np.cos(2 * y) -
        0.040849 * np.sin(2 * y))

    # longitudes in degrees, time_offset in minutes
    time_offset = eqtime + 4.0 * longitude

    # define pd.DatetimeIndex for solar time
    solar_time = timeUTC + pd.TimedeltaIndex(time_offset, unit='m')

    # convert to np.array with seconds as units
    return np.array(solar_time, dtype='datetime64[s]')