The sun's position¶

Its apparent right ascension and declination.

Requires

the time, as a Julian day (UT).

Provides

RightAscension and Declination (rad), the ecliptic longitude (rad) and the equation of time (s: the angle over a turn, of the Earth's mean solar day).

Law

the mean anomaly and mean longitude advancing daily, the equation of the centre to two terms, the obliquity of the ecliptic slowly falling; to about a hundredth of a degree over two centuries of J2000.

Source

the Astronomical Almanac's low-precision formulae (section C); Meeus, Astronomical Algorithms, ch. 25.

Proof

curve, reference, limit.

No Trait form: not written yet: its mean longitude is lusine.sky.mean_longitude; the equation of the centre, the obliquity and the turn from the ecliptic to the equator are not written, and a Julian day does not fit binary32 (sidereal's note); until then a reference.

sky.sun(jd)

Uses, without proving it here: sky.horizontal, sky.sidereal.

The kernel follows the Almanac's low-precision coordinates. Every number it runs on is a row's: the Sun's elements lusine.bodies.sun and the Earth's obliquity and day lusine.bodies.earth, bound as a game binds them: the kernel declares no constant.

In [1]:

import os, sys
root = os.getcwd()
while not os.path.isfile(os.path.join(root, 'docsite.json')):
    root = os.path.dirname(root)
sys.path.insert(0, os.path.join(root, 'tools'))
import numpy as np
import matplotlib.pyplot as plt
from lusine import dmath as D, library
from lusine.kernels.sky._astro import deg, julian_day, to_deg
from lusine.kernels.sky.horizontal import horizontal
from lusine.kernels.sky.sidereal import sidereal
from lusine.kernels.sky.sun import sun

import os, sys root = os.getcwd() while not os.path.isfile(os.path.join(root, 'docsite.json')): root = os.path.dirname(root) sys.path.insert(0, os.path.join(root, 'tools')) import numpy as np import matplotlib.pyplot as plt from lusine import dmath as D, library from lusine.kernels.sky._astro import deg, julian_day, to_deg from lusine.kernels.sky.horizontal import horizontal from lusine.kernels.sky.sidereal import sidereal from lusine.kernels.sky.sun import sun

The reference¶

The kernel's reference, in numpy: the body under the declaration above, and what every section below runs. It stays as the kernel's test beside its Trait form.

In [2]:

library.reference(sun)

library.reference(sun)

Out[2]:

the reference, sun.py

def sun(jd):
    """THE SUN'S POSITION: its apparent right ascension and declination.

    Requires: the time, as a Julian day (UT).
    Provides: `RightAscension` and `Declination` (rad), the ecliptic longitude
    (rad) and the equation of time (s: the angle over a turn, of the Earth's
    mean solar day).
    Law: the mean anomaly and mean longitude advancing daily, the equation of
    the centre to two terms, the obliquity of the ecliptic slowly falling; to
    about a hundredth of a degree over two centuries of J2000.
    Source: the Astronomical Almanac's low-precision formulae (section C);
    Meeus, Astronomical Algorithms, ch. 25.
    Proof: curve, reference, limit."""
    s = _SUN
    d = np.asarray(jd, dtype=np.float64) - s["epoch"]
    g = _wrap360(s["mean_anomaly_epoch"] + s["mean_anomaly_rate"] * d)
    q = _wrap360(s["mean_longitude_epoch"] + s["mean_longitude_rate"] * d)
    lam = _wrap360(q + s["centre_1"] * _s(g) + s["centre_2"] * _s(2.0 * g))
    e = _obliquity(np.asarray(jd, dtype=np.float64) - _EARTH["epoch"], _EARTH)
    ra = D.datan2(_c(e) * _s(lam), _c(lam))
    dec = D.dasin(_s(e) * _s(lam))
    eqt = q - _wrap360(to_deg(ra))
    eqt = (eqt + 180.0) - 360.0 * np.floor((eqt + 180.0) / 360.0) - 180.0
    return {"ra": ra, "dec": dec, "longitude": deg(lam), "equation_of_time": eqt / 360.0 * _EARTH["mean_solar_day"]}

Curve: the sun's path through a day, and its noon altitude¶

London's altitude through the solstices and an equinox; the noon altitudes are 90° - 51.5° ± 23.44° and 38.5°.

In [3]:

lat = deg(51.5)
hours = np.arange(0, 24, 0.05)
noon = {}
fig, ax = plt.subplots(figsize=(7, 3.5))
for label, (y, m, d) in {'21 June': (2024, 6, 20), '20 March': (2024, 3, 20), '21 December': (2024, 12, 21)}.items():
    jd = julian_day(y, m, d) + hours / 24
    s = sun(jd)
    alt, _ = horizontal(s['ra'], s['dec'], lat, sidereal(jd, 0.0))
    noon[label] = float(np.max(to_deg(alt)))
    ax.plot(hours, to_deg(alt), label=f'{label}: noon {noon[label]:.2f}°')
assert abs(noon['21 June'] - (90 - 51.5 + 23.44)) < 0.1 and abs(noon['21 December'] - (90 - 51.5 - 23.44)) < 0.1
assert abs(noon['20 March'] - (90 - 51.5)) < 0.5
ax.axhline(0, color='grey', lw=0.5)
ax.set_xlabel('hour (UT)')
ax.set_ylabel('altitude (°)')
ax.legend()
plt.show()

lat = deg(51.5) hours = np.arange(0, 24, 0.05) noon = {} fig, ax = plt.subplots(figsize=(7, 3.5)) for label, (y, m, d) in {'21 June': (2024, 6, 20), '20 March': (2024, 3, 20), '21 December': (2024, 12, 21)}.items(): jd = julian_day(y, m, d) + hours / 24 s = sun(jd) alt, _ = horizontal(s['ra'], s['dec'], lat, sidereal(jd, 0.0)) noon[label] = float(np.max(to_deg(alt))) ax.plot(hours, to_deg(alt), label=f'{label}: noon {noon[label]:.2f}°') assert abs(noon['21 June'] - (90 - 51.5 + 23.44)) < 0.1 and abs(noon['21 December'] - (90 - 51.5 - 23.44)) < 0.1 assert abs(noon['20 March'] - (90 - 51.5)) < 0.5 ax.axhline(0, color='grey', lw=0.5) ax.set_xlabel('hour (UT)') ax.set_ylabel('altitude (°)') ax.legend() plt.show()

Reference: the solstice's declination¶

At the June solstice of 2024 (20 June, 20:51 UT) the sun stands at the obliquity, 23.44° north.

In [4]:

dec = float(to_deg(sun(julian_day(2024, 6, 20, 20 + 51 / 60))['dec']))
print(f'declination {dec:.3f}°')
assert abs(dec - 23.44) < 0.01

dec = float(to_deg(sun(julian_day(2024, 6, 20, 20 + 51 / 60))['dec'])) print(f'declination {dec:.3f}°') assert abs(dec - 23.44) < 0.01

declination 23.436°

Limit: at an equinox the sun rises due east¶

With the sun on the celestial equator it rises at azimuth 90° from any latitude.

In [5]:

for latitude in (0.0, 30.0, 51.5, -40.0):
    jd = julian_day(2024, 3, 20) + np.arange(0, 1, 1 / 2880)
    s = sun(jd)
    alt, az = horizontal(s['ra'], s['dec'], deg(latitude), sidereal(jd, 0.0))
    rise = np.where((alt[:-1] < 0) & (alt[1:] >= 0))[0][0]
    print(f'latitude {latitude:5.1f}: rises at azimuth {float(to_deg(az[rise])):.2f}°')
    assert abs(float(to_deg(az[rise])) - 90.0) < 1.0

for latitude in (0.0, 30.0, 51.5, -40.0): jd = julian_day(2024, 3, 20) + np.arange(0, 1, 1 / 2880) s = sun(jd) alt, az = horizontal(s['ra'], s['dec'], deg(latitude), sidereal(jd, 0.0)) rise = np.where((alt[:-1] < 0) & (alt[1:] >= 0))[0][0] print(f'latitude {latitude:5.1f}: rises at azimuth {float(to_deg(az[rise])):.2f}°') assert abs(float(to_deg(az[rise])) - 90.0) < 1.0

latitude   0.0: rises at azimuth 89.95°
latitude  30.0: rises at azimuth 89.91°
latitude  51.5: rises at azimuth 89.89°
latitude -40.0: rises at azimuth 89.94°