Request: Calculation of solar beta angle for an orbit

[scottshambaugh]

The solar beta angle is a commonly used metric in spaceflight that is necessary to calculate eclipse time fractions and I think fits right in with skyfield's focus area.

It's calculated by:

$$ \beta = \sin^{-1} \Big[ \cos(\Gamma) \sin(\Omega) \sin(i) - \sin(\Gamma) \cos(\epsilon) \cos(\Omega) \sin(i) + \sin(\Gamma) \sin(\epsilon) \cos(i) \Big]$$

• $\beta$ — Solar beta angle
• $\Gamma$ — Ecliptic true solar longitude
• $\Omega$ — Right ascension of the ascending node (RAAN)
• $i$ — Orbit inclination
• $\epsilon$ — Obliquity of the ecliptic (approximately 23.45 deg for Earth at present)

It seems to me that all of this information is available for an EarthSatellite, so it would be great if there were a function or example of how to calculate this within skyfield!

A short paper with equations: Computation of Eclipse Time for Low-Earth Orbiting Small Satellites

[brandon-rhodes]

All those angles look quite daunting. Skyfield generally operates with vector math rather than trigonometry. Here, I've tried putting together a quick solution.

from numpy import cross
from skyfield import api
from skyfield.api import load, tau
from skyfield.functions import angle_between

ts = load.timescale()
eph = load('de421.bsp')
sun = eph['sun']
earth = eph['earth']
sat = api.EarthSatellite(
    '1 25544U 98067A   20114.53491251  .00001443  00000-0  33988-4 0  9995',
    '2 25544  51.6434 258.9578 0001655 179.6828 273.6120 15.49306451223542',
)
print('Satellite epoch time:', sat.epoch.utc_strftime())

# Figure out the orbital plane with positions one minute apart.

t0 = ts.utc(2020, 4, 23)
t1 = ts.utc(2020, 4, 23, 0, 1)

p0 = sat.at(t0).xyz.km
p1 = sat.at(t1).xyz.km
plane_normal = cross(p0, p1)

# Now, compute where the Sun is, and determine the angle up or down 90°
# from the orbital plane.

s = earth.at(t0).observe(sun).xyz.km

solar_beta_degrees = 90 - angle_between(plane_normal, s) / tau * 360.0
print('solar_beta_degrees:', solar_beta_degrees)

Try it out and see if its return value is anything close to the angle you need. Feel free to play around with some variations, and let me know if the results look meaningful or not.

[scottshambaugh]

The results there look right to me, thank you! What you've provided is cleaner than the trig version IMO. Tweaking your script a bit to show a time series and dashed red lines at the max possible beta angles (inclination + earth obliquity = 51.64 + 23.45 = 75.09 deg), this plot matches what I'd expect for an ISS TLE:

Compared to eg this plot from this paper (different dates):