How can I calculate if a satellite is currently in sunlight or eclipse using pyEphem or Skyfield (preferred)?

Question

I'd like to add this to some existing python scripts to model cubesat battery life and payload temperature but the math escapes me.

Answer 1

As a new user I cannot comment or tweak the original answer, so I'll try my own. Happy update per any recommendations that come up.

@uhoh's answer is close, but a few things to note.

The TLE tells us the number of orbits per day is 15.50995519 (line 2 columns 53–63)

$$ \frac{24 \frac{hours}{day} * 60 \frac{mins}{hour} }{15.50995519 \frac{orbits}{day}} = 92.8436 \frac{mins}{orbit}$$

So we expect our satellite to orbit once every ~93 minutes. Over an orbit, we expect to see the satellite transition into and out of eclipse once each.

@uhoh's answer shows two such transitions in that time. What's missing is consideration for the fact that there are two times when the satellite will be within the cone angle of Earth's limb. Both when it is in eclipse and directly in front of the Earth.

To build on @uhoh's answer, we can consider the distance Sun-Earth and Sun-Sat. The Sat will be in sun both when it is outside the Limb angle OR when it is closer to the Sun than the Earth.

Slightly modified code (thank's @uhoh for showing us the way to calculate Sun-Earth distance!)

Results:

To test the result, one can find the first occurrence of the transition True-False and False-True via:

sunlit[::].index(False)

and continue slicing:

sunlit[141+247::].index(False)

in this way we find transitions at 141 (True-False - Enter eclipse), then 247 (False-True - Exit eclipse), and 681 (True-False - Enter next Eclipse).

247 + 681 / 10 [each time step is 1/10th of a minute] = 92.8 minutes between eclipse entries.

TLE = """ISS (ZARYA)             
1 25544U 98067A   19203.81086311  .00000606  00000-0  18099-4 0  9996
2 25544  51.6423 184.5274 0006740 168.1171 264.4057 15.50995519180787"""

name, L1, L2 = TLE.splitlines()

import numpy as np
import matplotlib.pyplot as plt

from skyfield.api import Loader, EarthSatellite

halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]
degs, rads        = 180/pi, pi/180
Re                =   6378.137

load = Loader('~/Documents/fishing/SkyData')  # avoid duplicating large DE files
data = load('de421.bsp')
ts   = load.timescale(builtin=True)

minutes = np.arange(0, 200, 0.1)
times   = ts.utc(2019, 7, 23, 0, minutes)

data    = load('de421.bsp')
Earth   = data['earth']
Sun     = data['sun']
Sat     = Earth + EarthSatellite(L1, L2)

sunpos, earthpos, satpos = [thing.at(times).position.km for thing in (Sun, Earth, Sat)]


sunearth, sunsat         = earthpos-sunpos, satpos-sunpos

sunearthnorm, sunsatnorm = [vec/np.sqrt((vec**2).sum(axis=0)) for vec in (sunearth, sunsat)]

angle = np.arccos((sunearthnorm * sunsatnorm).sum(axis=0))

sunearthdistance = np.sqrt((sunearth**2).sum(axis=0))

sunsatdistance = np.sqrt((sunsat**2).sum(axis=0))

limbangle        = np.arctan2(6378.137, sunearthdistance)

sunlit = []
for idx, value in enumerate(angle):
    sunlit.append(((angle[idx] > limbangle[idx]) or (sunsatdistance[idx] < sunearthdistance[idx])))


if True:
    plt.figure(figsize=(10,11))

    plt.subplot(3, 1, 1)
    plt.plot(minutes, degs*limbangle, '--k')
    plt.plot(minutes, degs*angle, '-r', linewidth=2)
    plt.title('Earth-Sun-ISS and Earth-Sun-limb anges (degs)', fontsize=14)

    plt.subplot(3, 1, 2)
    plt.plot(minutes, sunearthdistance, '--k')
    plt.plot(minutes, sunsatdistance, '-r', linewidth=2)
    plt.title('Earth-Sun Distance and Sat-Sun distance (?)', fontsize=14)

    plt.subplot(3, 1, 3)
    plt.plot(minutes, sunlit, '-k')
    plt.title('ISS is sunlit (boolean)', fontsize=14)
    plt.ylim(-0.1, 1.1)
    plt.xlabel('minutes', fontsize=14)

    plt.show()

Answer 2

MEGA note: My answer is totally wrong and should be unaccepted by any means necessary. Luckily @AllenKummer's answer explains why and gives a revised script. This was posted at about 7 AM local time and must have been pre-coffee; the period shown here is only 40 minutes, seeing that it's wrong would have been a no-brainer, which is indeed the state I'm often in, pre-coffee.

note: I've answered first in Python rather than MathJax but I'll add the corresponding equations later today.

Simplest way is to calculate the opening angle as viewed from the Sun, or the Earth-Sun-Sat angle. Use a home made dot product of normalized vectors to get the angles, then calculate the Earth-Sun-Sat angle.

If the Earth-Sun-Sat angle is greater than the Earth-Sun-limb angle, then Sat is in daylight. If not, then it's in the Earth's shadow.

The transition from "day" to "night" (eclipse) is only a few seconds.

To make a 3D plot of the orbit in Python see this answer (Python-2) to How can I plot a satellite's orbit in 3D from a TLE using Python and Skyfield? and this Python-3 adjustment

Python-3:

TLE = """ISS (ZARYA)             
1 25544U 98067A   19203.81086311  .00000606  00000-0  18099-4 0  9996
2 25544  51.6423 184.5274 0006740 168.1171 264.4057 15.50995519180787"""

name, L1, L2 = TLE.splitlines()

import numpy as np
import matplotlib.pyplot as plt

from skyfield.api import Loader, EarthSatellite

halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]
degs, rads        = 180/pi, pi/180
Re                =   6378.137

load = Loader('~/Documents/fishing/SkyData')  # avoid duplicating large DE files
data = load('de421.bsp')
ts   = load.timescale()

minutes = np.arange(0, 200, 0.1)
times   = ts.utc(2019, 7, 23, 0, minutes)

data    = load('de421.bsp')
Earth   = data['earth']
Sun     = data['sun']
Sat     = Earth + EarthSatellite(L1, L2)

sunpos, earthpos, satpos = [thing.at(times).position.km for thing in (Sun, Earth, Sat)]

sunearth, sunsat         = earthpos-sunpos, satpos-sunpos

sunearthnorm, sunsatnorm = [vec/np.sqrt((vec**2).sum(axis=0)) for vec in (sunearth, sunsat)]

angle = np.arccos((sunearthnorm * sunsatnorm).sum(axis=0))

sunearthdistance = np.sqrt((sunearth**2).sum(axis=0))

limbangle        = np.arctan2(6378.137, sunearthdistance)

sunlit = angle > limbangle  # returns boolean array

if True:
    plt.figure()

    plt.subplot(2, 1, 1)
    plt.plot(minutes, degs*limbangle, '--k')
    plt.plot(minutes, degs*angle, '-r', linewidth=2)
    plt.title('Earth-Sun-ISS and Earth-Sun-limb anges (degs)', fontsize=14)

    plt.subplot(2, 1, 2)
    plt.plot(minutes, sunlit, '-k')
    plt.title('ISS is sunlit (boolean)', fontsize=14)
    plt.ylim(-0.1, 1.1)
    plt.xlabel('minutes', fontsize=14)

    plt.show()