Equatorial to horizontal¶

A body's altitude and azimuth.

Requires

the body's RightAscension and Declination (rad), the observer's Latitude (rad), the LocalSiderealTime (rad).

Provides

Altitude and Azimuth (rad; azimuth from north through east).

Law

with the hour angle H = LST - ra, sin alt = sin phi sin dec + cos phi cos dec cos H, and the azimuth from the atan2 of the same spherical triangle.

Source

Meeus, Astronomical Algorithms, ch. 13.

Proof

invariant, limit.

No Trait form: not written yet: a spherical triangle in binary32, waiting on nothing but a game that reads an altitude and an azimuth; until then a reference.

sky.horizontal(ra, dec, latitude, lst)

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, to_deg
from lusine.kernels.sky.horizontal import horizontal

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, to_deg from lusine.kernels.sky.horizontal import horizontal

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

library.reference(horizontal)

Out[2]:

the reference, horizontal.py

def horizontal(ra, dec, latitude, lst):
    """EQUATORIAL TO HORIZONTAL: a body's altitude and azimuth.

    Requires: the body's `RightAscension` and `Declination` (rad), the
    observer's `Latitude` (rad), the `LocalSiderealTime` (rad).
    Provides: `Altitude` and `Azimuth` (rad; azimuth from north through east).
    Law: with the hour angle H = LST - ra, sin alt = sin phi sin dec + cos phi
    cos dec cos H, and the azimuth from the atan2 of the same spherical
    triangle.
    Source: Meeus, Astronomical Algorithms, ch. 13.
    Proof: invariant, limit."""
    h = np.asarray(lst, dtype=np.float64) - np.asarray(ra, dtype=np.float64)
    phi = np.asarray(latitude, dtype=np.float64)
    dec = np.asarray(dec, dtype=np.float64)
    sin_alt = D.dsin(phi) * D.dsin(dec) + D.dcos(phi) * D.dcos(dec) * D.dcos(h)
    alt = D.dasin(np.clip(sin_alt, -1.0, 1.0))
    az = D.datan2(-D.dcos(dec) * D.dsin(h), D.dsin(dec) * D.dcos(phi) - D.dcos(dec) * D.dcos(h) * D.dsin(phi))
    return alt, az - 2.0 * D.PI * np.floor(az / (2.0 * D.PI))

Invariant: the pole stands at the latitude¶

The celestial pole's altitude is the observer's latitude, whatever the time and whatever the pole's right ascension.

In [3]:

for latitude in (-60.0, 0.0, 51.5, 89.0):
    for lst in np.linspace(0, 6.28, 7):
        alt, _ = horizontal(1.234, deg(90.0), deg(latitude), lst)
        assert abs(float(to_deg(alt)) - latitude) < 1e-9
print('the pole at the latitude')

for latitude in (-60.0, 0.0, 51.5, 89.0): for lst in np.linspace(0, 6.28, 7): alt, _ = horizontal(1.234, deg(90.0), deg(latitude), lst) assert abs(float(to_deg(alt)) - latitude) < 1e-9 print('the pole at the latitude')

the pole at the latitude

Limit: on the meridian, and overhead¶

A body on the meridian (hour angle 0) culminates at 90° - |latitude - declination|, due south or north; one at the latitude's declination passes overhead.

In [4]:

phi = deg(51.5)
for dec in (-20.0, 0.0, 30.0, 70.0):
    alt, az = horizontal(0.5, deg(dec), phi, 0.5)
    assert abs(float(to_deg(alt)) - (90 - abs(51.5 - dec))) < 1e-9
    assert abs(float(to_deg(az)) - (180.0 if dec < 51.5 else 0.0)) % 360 < 1e-6
# overhead the altitude is an arcsine of nearly 1, where a rounding in its sine is a microdegree of angle
alt, _ = horizontal(0.5, phi, phi, 0.5)
assert abs(float(to_deg(alt)) - 90) < 1e-4
print('culminations right')

phi = deg(51.5) for dec in (-20.0, 0.0, 30.0, 70.0): alt, az = horizontal(0.5, deg(dec), phi, 0.5) assert abs(float(to_deg(alt)) - (90 - abs(51.5 - dec))) < 1e-9 assert abs(float(to_deg(az)) - (180.0 if dec < 51.5 else 0.0)) % 360 < 1e-6 # overhead the altitude is an arcsine of nearly 1, where a rounding in its sine is a microdegree of angle alt, _ = horizontal(0.5, phi, phi, 0.5) assert abs(float(to_deg(alt)) - 90) < 1e-4 print('culminations right')

culminations right