How to determine how long Earth Satellite has been eclipsed?

[bismurphy]

I have a long history with PyEphem but I'm still getting used to Skyfield, and have been working on transitioning for years :)

I have a simple PyEphem script which will tell me whether a satellite is currently eclipsed, and in either case, how long it has been in its current state. I want to use this for the sake of checking the thermals and battery voltages of a satellite currently on-orbit (having sun timing context is useful for making sure these values are nominal).

In PyEphem, the script looks like (abbreviated):

    rise_time, rise_azimuth, max_alt_time, max_alt, set_time, set_azimuth=groundstation.next_pass(satellite)
    initial_eclipse_status = satellite.eclipsed
    while(satellite.eclipsed == initial_eclipse_status):
        #Step back one second and recompute, until we find a change.
        eclipse_time -= ephem.second
        satellite.compute(eclipse_time)
   time_since_eclipse_change = rise_time - eclipse_time

Now the initial eclipse status and the time difference tells me everything I'm looking for.

In Skyfield, I can't seem to find any such method for stepping back in time until it changes. The closest I can do is:

    sf_risetime_past_hour = ts.utc(*list(sf_risetime.tt_calendar()[:-1]),range(-3600,0,60))
    illum_array = sf_satellite.at(sf_risetime_past_hour).is_sunlit(eph)
    print(illum_array)

So once I know a rise time, I can step back for the past hour and check for any changes, but this is a little more clunky, and if the last change in eclipse was over an hour ago (my PyEphem solution is telling me that a particular pass I'm interested in will have been sunlit for 274 minutes continously), I learn nothing. Is there a way to do the Skyfield solution in a manner which would be as flexible as PyEphem?

[glangford]

I don't know if this is the best approach, but here is one idea (I also need to find illumination changes, not quite in the same way).

Could you step backwards multiple hours (up to some maximum n that makes sense) all at once - so range(0, -3600 * n, -60). ie. the range starts from 0 and time is reversed. Then something like

illum_changes = [i for i in range(1,len(illum_array)) if illum_array[i] != illum_array[i-1] ]

...will create an array of indices where the is_sunlit() status changed. We only care about the first value from illum_changes, say its value was 250, then an illumination change occurred 250 minutes ago.

To get finer resolution, additional code could do some date arithmetic and call is_sunlit() again for the block of 60 seconds within that minute.

With the 274 minute example you mention, going back in time second by second would loop over 16,000 times. In the above, if you look back say 5 hours in one time block it will evaluate just 300 different time points. Does that make sense?

EDIT To detect changes in illum_array you could use numpy.diff() rather than a list comprehension, if you prefer

[brandon-rhodes]

Great question, @tj-murphy!

I have never yet settled on a good API that I'll be happy to support forever for incrementing and decrementing times. Until I'm convinced of one, note that you can simply add or subtract increments from a time’s Terrestrial Time value, which is a value in Julian days. Thus:

one_second = 1.0 / 24 / 60 / 60
t = t.ts.tt_jd(t.tt - one_second)

— will decrement t by one second. And t.tt - 1.0 would decrement it by an entire day.

It would be nice to just be able to overwrite the value without creating a new Time object, but so many secondary values are cached on a time object that it’s vastly simpler to instead treat them as each read-only.

Skyfield does have a search routine that will use initial guesses followed by bisection to hone in on the moment when a value changes. My guess is that you could hook the routine up to is_sunlit() to do the search:

https://rhodesmill.org/skyfield/searches.html

Testing should quickly reveal which of these approaches is fastest given your hardware and specific satellite.

Let me know if these suggestions make sense, or if you will need some sample code!

[bismurphy]

Hmmm - your example of updating t by one second seems to make sense, but gives me TypeError: init() missing 1 required positional argument: 'ts'

Which makes sense, since we never fed it a timescale object to use. Given that I originally created my Time object with a ts.utc(example_datetime_object), I tried using ts.tt(t.tt - one_second) but that gives me an error saying that the ephemeris only covers the years 1899-2053 which makes me think somewhere in there the conversion is going wrong where it's misinterpreting the epoch of the two tt's we're dealing with.

[brandon-rhodes]

I should never share code without testing it first! I've updated the example code above, and here's an actual working example:

from skyfield.api import load

ts = load.timescale(builtin=True)
t = ts.utc(2020, 7, 10, 0, 19, 30)

one_second = 1.0 / 24 / 60 / 60
t = t.ts.tt_jd(t.tt - one_second)
print(t.utc_jpl())

Sorry for trying to wing it, I should have taken the extra few seconds to verify the code first!

[bismurphy]

Haha, no worries! I appreciate your willingness to help :) Transitioning from PyEphem to Skyfield has been quite a process for me, but the fact that you've made two awesome tools, and are super helpful along the way, has made it much more manageable. So, thanks!

[bismurphy]

Actually - I seem to have found a massive difference between answers I get with PyEphem code, and Skyfield code, which ought to be doing almost identical things. I ended up making a test script which indeed finds different answers - but only sometimes.

from skyfield.api import Topos,load,utc,Time,EarthSatellite
import ephem

import re #used for converting date_string into a comma-separated list which Skyfield wants

#Works great!
#date_string = '2020/7/13 10:05:44'
#This one fails.
date_string = '2020/7/14 01:12:35'
#ISS TLE.Source: https://spaceflight.nasa.gov/realdata/sightings/SSapplications/Post/JavaSSOP/orbit/ISS/SVPOST.html
tle_line1 = '1 25544U 98067A   20195.20691215  .00016717  00000-0  10270-3 0  9088'
tle_line2 = '2 25544  51.6401 219.6922 0001439 124.7799 235.3486 15.49521317 36048'
#Doing everything in ephem. All variables will start with e_ to prevent mixups.

e_time = ephem.date(date_string)
e_satellite = ephem.readtle('ISS',tle_line1,tle_line2)
e_satellite.compute(e_time)
e_init_eclipse = e_satellite.eclipsed
#This variable holds each candidate time for eclipse to transition. Start with our e_time and work back.
e_test_time = e_time

while(e_satellite.eclipsed == e_init_eclipse):
    e_test_time -= ephem.second
    e_satellite.compute(e_test_time)
#e_test_time now holds the exact moment the most recent eclipse transition happened.
e_min_since_change = (e_time - e_test_time)/ephem.minute
print(("ENTERED" if e_init_eclipse else "EXITED") + " eclipse " + f"{(e_min_since_change):.1f}" + " min ago")

#Now the same, in Skyfield. All variables start with s_ again to prevent crosstalk between the sections' methods.
s_ts = load.timescale(builtin=True)
#Converts the date_string into a list of integers, then pushes that list into s_ts.utc. Somewhat clunky
#but allows for Skyfield and ephem to take in the same input.
s_time = s_ts.utc(*[int(i) for i in re.split('[/ :]',date_string)])
s_eph=load('de421.bsp')
s_satellite = EarthSatellite(tle_line1, tle_line2,'ISS',s_ts)
s_init_illum = s_satellite.at(s_time).is_sunlit(s_eph)
s_test_time = s_time

while(s_init_illum == s_satellite.at(s_test_time).is_sunlit(s_eph)):
    #Need to go through the timescale object. Also, Skyfield does not have constants like ephem.second.
    s_test_time = s_test_time.ts.tt_jd(s_test_time.tt - 1/86400)
s_min_since_change = (s_time - s_test_time) *24*60
print(("EXITED" if s_init_illum else "ENTERED") + " eclipse " + f"{s_min_since_change:.1f}" + " min ago")

When running this as-is, it gives 12.0 and 12.1 - close enough for my purposes. But if you swap the commenting at the top to use the other date string (these are just two times where my location has passes occurring, as found by a separate script - but nothing here is observer-centric, so it doesn't matter where the dates came from), then you'll find that the PyEphem portion gives 178.4 minutes as the answer (which is already strange, since the ISS orbits every 90-ish minutes, and isn't at a high enough inclination to do a full sunlit orbit, as far as I know...), and then Skyfield gives the answer of 363.8 minutes. Furthermore, on my computer (relatively beefy), the PyEphem answer comes effectively instantly while the Skyfield answer takes 50 seconds to find. I'm not complaining about the speed because I'm sure there are time-consuming operations involved with transforming the time objects, and if I iterated every minute instead of every second it would be much faster. It was just striking to me how different the two libraries are in terms of execution time. But the bigger concern I have is that these two methods gave entirely different answers.

Did I mess something up? It doesn't feel like it, especially since the first test case (and actually, I've had multiple test cases) worked, but suddenly the second one doesn't. If the two methods disagree only some of the time, that makes me suspect it's something in the underlying functions.

[glangford]

For comparison purposes, this code produces the same conflicting result but the performance is much better. With 1s resolution looking back 10 hours in time, running both PyEphem and Skyfield takes 2s.

# Modified test
s_ts = load.timescale(builtin=True)
s_satellite = EarthSatellite(tle_line1, tle_line2,'ISS',s_ts)
s_eph=load('de421.bsp')
s_risetime_past_hours = s_ts.utc(2020, 7, 14, 1, 12 , range(35, -3600 * 10, -1) )
illum_array = s_satellite.at(s_risetime_past_hours).is_sunlit(s_eph)
illum_changes = [i for i in range(1,len(illum_array)) if illum_array[i] != illum_array[i-1] ]
if len(illum_changes) > 0:
    print( illum_changes[0] / 60., " minutes, ", )

displays 363.76666666666665 minutes

[glangford]

@brandon-rhodes I am getting weird results depending on how I set the range() - either forwards or backwards to the problem date. Could something be off in the time internals triggered by the fact that the subtraction goes back past midnight to the previous day?

[glangford]

then you'll find that the PyEphem portion gives 178.4 minutes as the answer (which is already strange, since the ISS orbits every 90-ish minutes, and isn't at a high enough inclination to do a full sunlit orbit, as far as I know...)

Apparently ISS can enter a full sun period, but I don't know if we are close to one of them: https://space.stackexchange.com/questions/4686/how-often-does-the-iss-orbit-align-with-the-day-night-terminator

[bismurphy]

Looks like you're right about that!

https://spaceflight.nasa.gov/feedback/expert/answer/mcc/sts-113/11_23_20_01_179.html

This page also talks about ISS beta angles. I'm not certain how to get beta angles in Skyfield, but I punched the ISS into GMAT, and found a beta angle of about 68 degrees. This does indeed indicate that the ISS is in a short period of full-sun.

Still, the disagreement between the two libraries is unusual. I don't suppose it just comes down to the ISS barely grazing the edge of the earth's shadow for a brief moment?

[brandon-rhodes]

@glangford — No, there should, happily, not be any sensitivity to midnight in the time logic. But feel free to investigate further if you do ever print out a time and see a date and time that surprise you.

Thank you for offering vectorized code, as it is indeed much faster than separately building tens of thousands of separate time objects! When I try your suggested code, I see a drop from 22s to 1.1s on my machine.

One substantial cost in either case is the high precision nutation calculation that lets Skyfield produce Sun positions that are accurate enough for modern astronomers. Alas, all that precision is wasted when trying to merely determine whether a satellite is in shadow: getting the Sun's position correctly to a fraction of a milliarcsecond makes no difference if we're merely measuring to a resolution of one second the time when a satellite comes into or out of shadow.

I'm currently working on an API that will let folks control the nutation model expense. But you can secretly get around it even today with the lines:

from skyfield.nutationlib import iau2000b_radians
# ...
s_risetime_past_hours._nutation_angles = iau2000b_radians(s_risetime_past_hours)

With that improvement, your approach more than 3× faster, clocking in at around 0.3s instead of 1.1s.

@tj-murphy — I'll look into how the two libraries compute sun-lit-ness and see if I can learn the difference.

[brandon-rhodes]

@tj-murphy — Computing per-minute to make the script run faster, let's look at Skyfield's is_sunlit() compared to PyEphem’s eclipsed over the 600 minutes considered by the above script. Even better, let's also graph the angle between the Sun's center and the Earth's center from the ISS's point of view over that period.

As you can see, the only issue here seems to be a very slight difference of opinion between Skyfield and PyEphem about at what exact angle the cut-off is that hides the Sun behind the Earth. As the ISS's orbit turns sun-ward, and the "nights" get shorter and shorter — right around the time that the minima in the Earth-Sun separation reach and then surpass 70° — Skyfield is the first to guess that, no, there are no more full sunsets for the ISS. PyEphem, by contrast, takes 2 more orbits before deciding similarly that the Earth is no longer eclipsing the Sun.

Skyfield if anything should be overly generous in assessing sun-lit-ness because, looking at the code, I'm not clear whether it accounts for the Sun’s large size in the sky — it only counts a body as sunlit if the center of the Sun is visible for the body, even though passengers on the body would say the Sun is "up" once any part of the Sun is visible. I guess I figured that, for satellite illumination, the moment when it's halfway illuminated was as good a moment to flip the Boolean as any?

Feel free to dig into how both libraries compute the exact angle. Even better: do you have records from the actual ISS or another satellite for when it measured itself as in sunlight, that we could compare with both libraries?

Here's my script that produced the above plot:

import ephem
import matplotlib.pyplot as plt
from skyfield.api import Topos,load,utc,Time,EarthSatellite
from skyfield.nutationlib import iau2000b_radians

# #ISS TLE.Source: https://spaceflight.nasa.gov/realdata/sightings/SSapplications/Post/JavaSSOP/orbit/ISS/SVPOST.html
date_string = '2020/7/14 01:12:35'
tle_line1 = '1 25544U 98067A   20195.20691215  .00016717  00000-0  10270-3 0  9088'
tle_line2 = '2 25544  51.6401 219.6922 0001439 124.7799 235.3486 15.49521317 36048'

e_time = ephem.date(date_string)
e_satellite = ephem.readtle('ISS',tle_line1,tle_line2)
e_array = []
e_test_time = e_time
for i in range(600):
    e_test_time -= ephem.minute
    e_satellite.compute(e_test_time)
    e_array.append(e_satellite.eclipsed)

s_ts = load.timescale(builtin=True)
s_satellite = EarthSatellite(tle_line1, tle_line2,'ISS',s_ts)
s_eph=load('de421.bsp')
s_risetime_past_hours = s_ts.utc(2020, 7, 14, 1, range(12, 12 - 600, -1))

t = s_risetime_past_hours
s_risetime_past_hours._nutation_angles = iau2000b_radians(s_risetime_past_hours)
illum_array = s_satellite.at(s_risetime_past_hours).is_sunlit(s_eph)

earth = - s_satellite.at(t)
sun = (s_eph['sun'] - s_eph['earth'] - s_satellite).at(t)
sep = earth.separation_from(sun)

minutes = (t.tt - t.tt[0]) * 24.0 * 60.0

fig, (ax, ax2, ax3) = plt.subplots(3, 1)
ax.plot(minutes, sep.degrees, label='Angle between Sun/Earth', linestyle='--')
ax.set(ylabel='Degrees', title='Sun/Earth Separation from the ISS’s viewpoint')
ax.grid()

ax2.plot(minutes, illum_array)
ax2.set(ylabel='is_sunlit')

ax3.plot(minutes, e_array)
ax3.set(xlabel='Minutes', ylabel='eclipsed')

fig.savefig('tmp.png')

[glangford]

Extending the time window out a few days, both libraries agree that there will be a long period where ISS will be sunlit.

[glangford]

...and looking backwards from date_string = '2020/7/19 0:0:00' both libraries line up when ISS returns to 'normal' and is no longer in full sun

[glangford]

@tj-murphy Thinking back to your original use case, I was curious about the longer term view and the seasonal variation of light exposure. For fun I did this plot of ISS over a full year, using is_sunlit(). The blue plot oscillates since it is the raw percentage of sunlight over the orbits in 12 hour blocks, with a smoothed line as well. Not what I expected!

[brandon-rhodes]

Thanks for all of the great diagrams, @glangford!

@tj-murphy — At this point is this issue answered to your satisfaction, or do you have further lingering questions we could address?

[bismurphy]

Go ahead and close the issue.

Thanks everyone for the massive learning experience! I got way more than I bargained for :)

[glangford]

Thanks from me as well, this was a good question and a real learning experience for me too. @tj-murphy you mentioned solar beta and I ended up learning a great deal looking into that too. Here is one final plot showing the agreement between solar beta angle for ISS (in red) versus the theoretical sunlight % (light blue) vs is_sunlit() measured over 24 hour periods (dark blue dots). It doesn't line up exactly because it is measured in time, not in discrete orbits but the agreement is perfect at peaks. It's interesting to see how quickly the transition to full sun happens.

[bismurphy]

@glangford Any chance you'd be inclined to paste source code for that plot? It looks super neat and I might want to play with other satellites/other years and such. Totally cool if you want to let your code remain your own, but thought I'd ask :)

Thanks for all your help!

[glangford]

@tj-murphy Sure, happy to share it - allow me to clean up the code a little, then I will paste it in this thread.

[glangford]

#!/usr/bin/env python3
"""
MIT License

Copyright (c) 2020 Glenn Langford

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
"""

import matplotlib.pyplot as plt
import matplotlib.dates as mdates
from skyfield.api import Topos, load, Time, EarthSatellite
from skyfield.nutationlib import iau2000b_radians
from skyfield.elementslib import osculating_elements_of
from skyfield.constants import ERAD
from numpy import count_nonzero, cos, sin, arcsin, degrees, average, sqrt, square, clip
from scipy.signal import savgol_filter
import datetime
import math

tle_line1 = '1 25544U 98067A   20192.51664795  .00004501  00000-0  89168-4 0  9993'
tle_line2 = '2 25544  51.6438 233.0051 0002955  88.6349  16.2102 15.49210173235636'

# Calculate the Solar beta angle for a satellite at time t (in degrees), and the corresponding
# proportion of sunlit time (expressed as a percentage).
def solar_beta(sat, sun, earth, t):
    sun_astro = earth.at(t).observe(sun)
    sun_ra, sun_dec, sun_distance = sun_astro.radec()
    elements = osculating_elements_of(sat.at(t))
    inclination = elements.inclination.radians
    raan = elements.longitude_of_ascending_node.radians
    pre_beta = cos(sun_dec.radians) * sin(inclination) * sin(raan - sun_ra.radians) + \
                sin(sun_dec.radians) * cos(inclination)
    beta = arcsin(clip(pre_beta, -1, 1))
    sun_prop = sun_proportion(beta, elements)
    return degrees(beta), sun_prop

# https://stackoverflow.com/questions/48235232/valueerror-math-domain-error-from-math-acos-and-nan-from-numpy-arccos
# For a given solar beta angle, return the corresponding theoretical proportion
# of time in sun per orbit
def sun_proportion(solar_beta, osc_elements):
    sm_axis = osc_elements.semi_major_axis.m # in meters, consistent with ERAD
    pre_sun = sqrt(1 - square(ERAD / sm_axis)) / cos(solar_beta)
    pre_sun = clip(pre_sun, -1, 1)
    return (0.5 + 1 / math.pi * arcsin(pre_sun)) * 100

ts = load.timescale(builtin=True)
satellite = EarthSatellite(tle_line1, tle_line2, 'ISS', ts)
eph = load('de421.bsp')
sun = eph['sun']
earth = eph['earth']
sunlit_by_days = []
solar_beta_by_days = []
solar_prop_by_days = []

days = []
for i in range(1, 366):
    dayrange = ts.utc(2020, 1, i, 0, 0, range(0, 60*60*24, 10))
    dayrange._nutation_angles = iau2000b_radians(dayrange)
    illum_array = satellite.at(dayrange).is_sunlit(eph)
    #half = int(len(illum_array) / 2)
    sunlit_pct = count_nonzero(illum_array) / len(illum_array) * 100.0
    #sunlit_pct = count_nonzero(illum_array[:half]) / half * 100.0
    sunlit_by_days.append(sunlit_pct)
    beta, sunprop = solar_beta(satellite, sun, earth, dayrange)
    if i % 10 == 0:
        print(i, beta[0], sunprop[0])
    solar_beta_by_days.append(average(beta))
    solar_prop_by_days.append(average(sunprop))
    days.append(i)

# https://stackoverflow.com/questions/46633544/smoothing-out-curve-in-python
# https://en.wikipedia.org/wiki/Savitzky–Golay_filter
# https://stackoverflow.com/questions/20618804/how-to-smooth-a-curve-in-the-right-way
yhat = savgol_filter(sunlit_by_days, 7, 2) # window size, polynomial order
#
start = datetime.datetime(2020, 1, 1)
end = start + datetime.timedelta(days=365)
year_drange = mdates.drange(start, end, datetime.timedelta(hours=24))
#
fig, ax  = plt.subplots()
fig.set_size_inches(16, 9)
ax.grid()
#
ax.plot(year_drange, sunlit_by_days, '.')
#ax.plot(year_drange, yhat)
ax.plot(year_drange, solar_prop_by_days, color='lightskyblue')

ax.xaxis.set_major_locator(mdates.MonthLocator())
ax.xaxis.set_major_formatter(mdates.DateFormatter('%b'))
plt.xticks(rotation='vertical')
ax.set_xlim(start, end)
#
# Solar Beta plot
#
ax2 = ax.twinx()  # create a second axis that shares the same x
color = 'tab:red'
ax2.set_ylabel('Solar beta angle', color=color)
ax2.plot(year_drange, solar_beta_by_days, color=color)
ax2.tick_params(axis='y', labelcolor=color)
ax2.hlines(y=70, xmin=start, xmax=end, linestyles='dotted', color=color)
ax2.hlines(y=-70, xmin=start, xmax=end, linestyles='dotted', color=color)
ax2.xaxis.set_major_locator(mdates.MonthLocator())
ax2.xaxis.set_major_formatter(mdates.DateFormatter('%b'))
ax.set(ylabel='% Sun', xlabel='', title='2020 ISS Illumination')
fig.savefig('iss-solar.png', orientation='landscape')