Merge pull request #15 from AssessingSolar/Add-solar-position-models

Add solar position models

Author: AdamRJensen

.github/workflows/lint.yml

@@ -22,4 +22,4 @@ jobs:
         pip install flake8
     - name: Lint with flake8
       run: |
-        flake8 . --count --max-line-length=99 --statistics
+        flake8 . --count --max-line-length=99 --statistics --ignore=E741

docs/source/solarposition.rst

@@ -8,3 +8,9 @@ Solar position algorithms
    :toctree: generated/
 
    solarposition.iqbal
+   solarposition.michalsky
+   solarposition.noaa
+   solarposition.psa
+   solarposition.sg2
+   solarposition.usno
+   solarposition.walraven

src/solposx/refraction/spa.py

@@ -44,7 +44,6 @@ def spa(elevation, pressure=101325., temperature=12., refraction_limit=0.5667):
        Applications (Revised)." 2008. NREL Report No. TP-560-34302, pp. 55
        :doi:`10.2172/15003974`.
     """  # noqa: #501
-
     pressure = pressure / 100  # convert to hPa
     # switch sets elevation when the sun is below the horizon
     switch = elevation >= -1.0 * (0.26667 + refraction_limit)

src/solposx/solarposition/__init__.py

@@ -1 +1,7 @@
 from solposx.solarposition.iqbal import iqbal  # noqa: F401
+from solposx.solarposition.michalsky import michalsky  # noqa: F401
+from solposx.solarposition.noaa import noaa  # noqa: F401
+from solposx.solarposition.psa import psa  # noqa: F401
+from solposx.solarposition.sg2 import sg2  # noqa: F401
+from solposx.solarposition.usno import usno  # noqa: F401
+from solposx.solarposition.walraven import walraven  # noqa: F401

src/solposx/solarposition/michalsky.py

@@ -0,0 +1,174 @@
+from pvlib.tools import sind, cosd, asind
+import numpy as np
+import pandas as pd
+from solposx import refraction
+from solposx.tools import _pandas_to_utc, _fractional_hour
+
+
+def michalsky(times, latitude, longitude, spencer_correction=True,
+              julian_date='original'):
+    """
+    Calculate solar position using Michalsky's algorithm.
+
+    Michalsky's algorithm [1]_ has a stated accuracy of 0.01 degrees
+    from 1950 to 2050.
+
+    Parameters
+    ----------
+    times : pandas.DatetimeIndex
+        Must be localized or UTC will be assumed.
+    latitude : float
+        Latitude in decimal degrees. Positive north of equator, negative
+        to south. [degrees]
+    longitude : float
+        Longitude in decimal degrees. Positive east of prime meridian,
+        negative to west. [degrees]
+    spencer_correction : bool, default True
+        Applies the correction suggested by Spencer [2]_ so the algorithm
+        works for all latitudes.
+    julian_date : string, default 'original'
+        Julian date calculation. Can be one of the following:
+            * ``'original'``: calculation based on Michalsky's paper [1]_.
+            * ``'pandas'``: calculation using a pandas build-in function
+
+    Returns
+    -------
+    DataFrame with the following columns (all values in degrees):
+
+        * elevation : actual sun elevation (not accounting for refraction).
+        * apparent_elevation : sun elevation, accounting for
+          atmospheric refraction.
+        * zenith : actual sun zenith (not accounting for refraction).
+        * apparent_zenith : sun zenith, accounting for atmospheric
+          refraction.
+        * azimuth : sun azimuth, east of north.
+
+    Raises
+    ------
+    ValueError
+        An error is raised if the julian_date calculation is not `original`
+        or `pandas`.
+
+    Notes
+    -----
+    The Michalsky algorithm is based equations in the Astronomical Almanac.
+
+    As pointed out by Spencer (1989) [2]_, the original Michalsky
+    algorithm [1]_ did not work for the southern hemisphere. This
+    implementation includes by default the correction provided by Spencer such
+    that it works for all latitudes.
+
+    Minor clarifications were made to the original paper have been published
+    as errata in [3]_ and [4]_.
+
+    The Julian date calculation in the original Michalsky paper [1]_ ensures
+    the stated accuracy (0.01 degrees) only for the period 1950 - 2050. Outside
+    this period, the Julian date does not handle leap years correctly. Thus,
+    there is the option to use the
+    :py:func:`pandas.DatetimeIndex.to_julian_date` method.
+
+    References
+    ----------
+    .. [1] J. J. Michalsky, "The Astronomical Almanac’s algorithm for
+       approximate solar position (1950–2050)," Solar Energy, vol. 40, no. 3.
+       Elsevier BV, pp. 227–235, 1988. :doi:`10.1016/0038-092x(88)90045-x`.
+    .. [2] J. W. Spencer, “Comments on The Astronomical Almanac’s Algorithm for
+       Approximate Solar Position (1950–2050),” Solar Energy, vol. 42, no. 4.
+       Elsevier BV, p. 353, 1989. :doi:`10.1016/0038-092x(89)90039-x`.
+    .. [3] J. J. Michalsky, “Errata,” Solar Energy, vol. 41, no. 1. Elsevier
+       BV, p. 113, 1988. :doi:`10.1016/0038-092x(88)90122-3`.
+    .. [4] J. J. Michalsky, “Errata,” Solar Energy, vol. 43, no. 5. Elsevier
+       BV, p. 323, 1989. :doi:`10.1016/0038-092x(89)90122-9`.
+    """
+    times_utc = _pandas_to_utc(times)
+
+    hour = _fractional_hour(times_utc)
+
+    if julian_date == 'original':
+        year = times_utc.year
+        day = times_utc.dayofyear
+        delta = year - 1949
+        leap = np.floor(delta / 4)
+        jd = 2432916.5 + delta * 365 + leap + day + hour / 24
+    elif julian_date == 'pandas':
+        jd = times_utc.to_julian_date()
+    else:
+        raise ValueError(
+            "Either `original` or `pandas` has to be chosen for the Julian"
+            " date calculation.")
+
+    n = jd - 2451545.0
+
+    # L - mean longitude [degrees]
+    L = 280.460 + 0.9856474 * n
+    # L has to be between 0 and 360 deg
+    L = L % 360
+
+    # g - mean anomaly [degrees]
+    g = 357.528 + 0.9856003 * n
+    # g has to be between 0 and 360 deg
+    g = g % 360
+
+    # l - ecliptic longitude [degrees]
+    l = L + 1.915 * sind(g) + 0.02 * sind(2*g)
+    # l has to be between 0 and 360 deg
+    l = l % 360
+
+    # ep - obliquity of the ecliptic [degrees]
+    ep = 23.439 - 0.0000004 * n
+
+    # ra - right ascension [degrees]
+    ra = np.rad2deg(np.arctan2(cosd(ep) * sind(l), cosd(l)))
+    # ra has to be between 0 and 360 deg
+    ra = ra % 360
+
+    # dec - declination angle [degrees]
+    dec = asind(sind(ep) * sind(l))
+
+    # gmst - Greenwich mean sidereal time [hr]
+    gmst = 6.697375 + 0.0657098242 * n + hour
+    # gmst has to be between 0 and 24 h
+    gmst = gmst % 24
+
+    # lmst - local mean sideral time [hr]
+    lmst = gmst + longitude / 15  # to convert deg to h, divide with 15
+    # lmst has to be between 0 and 24 h
+    lmst = lmst % 24
+
+    # ha - hour angle [hr]
+    ha = lmst - ra / 15
+    # ha has to be between -12 and 12 h
+    ha = (ha + 12) % 24 - 12
+
+    # el - solar elevation angle [degrees]
+    el = asind(sind(dec) * sind(latitude) + cosd(dec) * cosd(latitude)
+               * cosd(15 * ha))  # to convert h to deg, multiply with 15
+
+    # az - azimuth [degrees]
+    az = asind(-cosd(dec) * sind(15 * ha) / cosd(el))
+
+    if spencer_correction:
+        # Spencer correction for the azimuth quadrant assignment
+        cos_az = sind(dec) - sind(el) * sind(latitude)
+        az = np.where((cos_az >= 0) & (sind(az) < 0), 360 + az, az)
+        az = np.where(cos_az < 0, 180 - az, az)
+    else:
+        # Original Michalsky implementation does not work for all latitudes
+        # calcualte critical elevation
+        elc = asind(sind(dec) / sind(latitude))
+        # correct azimuth using critical elevation
+        az = np.where(el >= elc, 180 - az, az)
+        az = np.where((el <= elc) & (ha > 0), az + 360, az)
+        az = az % 360
+
+    # refraction correction
+    r = refraction.michalsky(el)
+
+    result = pd.DataFrame({
+        'elevation': el,
+        'apparent_elevation': el + r,
+        'zenith': 90 - el,
+        'apparent_zenith': 90 - (el + r),
+        'azimuth': az,
+    }, index=times)
+    return result

src/solposx/solarposition/noaa.py

@@ -0,0 +1,137 @@
+import pandas as pd
+from pvlib.tools import sind, cosd, asind, acosd, tand
+import numpy as np
+from solposx import refraction
+from solposx.tools import _pandas_to_utc, _fractional_hour
+
+
+def noaa(times, latitude, longitude, delta_t=67.0):
+    """
+    Calculate solar position using NOAA's algorithm.
+
+    NOAA's algorithm [1]_ has a stated accuracy of 0.0167 degrees
+    from years -2000 to +3000 for latitudes within +/- 72 degrees. For
+    latitudes outside this the accuracy is 0.167 degrees.
+
+    The NOAA algorithm uses by default the Hughes refraction model,
+    see :py:func:`package_name.refraction.hughes`.
+
+    Parameters
+    ----------
+    times : pandas.DatetimeIndex
+        Must be localized or UTC will be assumed.
+    latitude : float
+        Latitude in decimal degrees. Positive north of equator, negative
+        to south. [degrees]
+    longitude : float
+        Longitude in decimal degrees. Positive east of prime meridian,
+        negative to west. [degrees]
+    delta_t : float or array, optional, default 67.0
+        Difference between terrestrial time and UT1.
+        If delta_t is None, uses spa.calculate_deltat
+        using time.year and time.month from pandas.DatetimeIndex.
+        For most simulations the default delta_t is sufficient.
+        The USNO has historical and forecasted delta_t [3]_.
+        [seconds]
+
+    Returns
+    -------
+    DataFrame with the following columns (all values in degrees):
+
+        * elevation : actual sun elevation (not accounting for refraction).
+        * apparent_elevation : sun elevation, accounting for
+          atmospheric refraction.
+        * zenith : actual sun zenith (not accounting for refraction).
+        * apparent_zenith : sun zenith, accounting for atmospheric
+          refraction.
+        * azimuth : sun azimuth, east of north.
+
+    Notes
+    -----
+    The algorithm is from Astronomical Algorithms by Jean Meeus [2]_.
+
+    References
+    ----------
+    .. [1] Solar Calculation Details. Global Monitoring
+       Laboratory Earth System Research Laboratories.
+       https://gml.noaa.gov/grad/solcalc/calcdetails.html
+    .. [2] Meeus, J. "Astronomical Algorithms", 1991.
+       https://archive.org/details/astronomicalalgorithmsjeanmeeus1991
+    .. [3] USNO delta T:
+       https://maia.usno.navy.mil/products/deltaT
+    """
+    times_utc = _pandas_to_utc(times)
+    julian_date = times_utc.to_julian_date()
+    jc = (julian_date - 2451545) / 36525
+
+    # [degrees]
+    mean_long = (
+        280.46646 + jc*(36000.76983 + jc*0.0003032)
+    ) % 360
+
+    mean_anom = 357.52911 + jc * (35999.05029 - 0.0001537 * jc)
+
+    eccent_earth_orbit = (
+        0.016708634
+        - jc * (0.000042037 + 0.0000001267 * jc))
+
+    sun_eq_ctr = (
+        sind(mean_anom)*(1.914602 - jc*(0.004817 + 0.000014*jc))
+        + sind(2*mean_anom)*(0.019993 - 0.000101*jc)
+        + sind(3*mean_anom)*0.000289
+    )
+
+    sun_true_long = mean_long + sun_eq_ctr
+
+    sun_app_long = sun_true_long - 0.00569 - 0.00478*sind(125.04 - 1934.136*jc)
+
+    mean_obliq_ecliptic = 23 + (26 + (
+        21.448 - jc*(46.815 + jc*(0.00059 - jc*0.001813)))/60)/60
+
+    obliq_corr = mean_obliq_ecliptic + 0.00256*cosd(125.04 - 1934.136*jc)
+
+    sun_declin = asind(sind(obliq_corr) * sind(sun_app_long))
+
+    var_y = tand(obliq_corr / 2)**2
+
+    eot = 4*np.degrees(
+        + var_y*sind(2*mean_long)
+        - 2*eccent_earth_orbit*sind(mean_anom)
+        + 4*eccent_earth_orbit*var_y*sind(mean_anom)*cosd(2*mean_long)
+        - 0.5*(var_y**2)*sind(4*mean_long)
+        - 1.25*(eccent_earth_orbit**2)*sind(2*mean_anom)
+    )
+
+    minutes = _fractional_hour(times_utc)*60
+
+    true_solar_time = (minutes + eot + 4*longitude) % 1440
+
+    hour_angle = np.where(
+        true_solar_time / 4 < 0,
+        true_solar_time / 4 + 180,
+        true_solar_time / 4 - 180,
+    )
+
+    zenith = acosd(sind(latitude)*sind(sun_declin) +
+                   cosd(latitude)*cosd(sun_declin)*cosd(hour_angle))
+
+    azimuth = np.where(
+        hour_angle > 0,
+        (acosd(((sind(latitude)*cosd(zenith)) - sind(sun_declin)) /
+               (cosd(latitude)*sind(zenith)))+180 % 360),
+        (540 - acosd(((sind(latitude)*cosd(zenith)) - sind(sun_declin)) /
+                     (cosd(latitude)*sind(zenith)))) % 360,
+    )
+
+    elevation = 90 - zenith
+    refraction_correction = refraction.hughes(
+        elevation=elevation, pressure=101325, temperature=10)
+
+    result = pd.DataFrame({
+        'elevation': elevation,
+        'apparent_elevation': elevation + refraction_correction,
+        'zenith': zenith,
+        'apparent_zenith': zenith - refraction_correction,
+        'azimuth': azimuth,
+    }, index=times)
+    return result