4
$\begingroup$

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

$\endgroup$
  • 1
    $\begingroup$ My answer is wrong! It would be great if you can un-accept it and then accept the correct answer. Thanks! $\endgroup$ – uhoh Dec 3 '19 at 5:13
7
+200
$\begingroup$

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: ISS In Sun considering Distance

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()
| improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ Holy granola! My answer is totally wrong, and I've edited it and made a note of it. I understand that you tried to edit my answer, but I think your first post is quite welcomed and valuable. Welcome to Space! $\endgroup$ – uhoh Dec 3 '19 at 5:13
2
$\begingroup$

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

ISS is in Sunshine!

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()
| improve this answer | |
$\endgroup$
  • 2
    $\begingroup$ To make the code slightly briefer, you might be able to use the separation_from() method to compute the angle between two vectors. $\endgroup$ – Brandon Rhodes Jul 27 '19 at 0:49
  • 1
    $\begingroup$ @BrandonRhodes good point! You'd mentioned that here quite a while ago, and then I'd used it here but after that it seems that I'd forgotten about it. Thanks for the pointer! $\endgroup$ – uhoh Jul 27 '19 at 6:50
  • $\begingroup$ @AllenKummer I see you've made an edit but i's been rejected by my and another user as an "attempt to reply". The better way to go about this is to leave a comment explaining what you think may be wrong, rather than overwriting another user's answer. In this case though, unfortunately, being a new user with less than 50 reputation points I think Stack Exchange doesn't allow you to leave comments, so you can ping me in The Pod Bay chatroom. $\endgroup$ – uhoh Dec 2 '19 at 0:11
  • $\begingroup$ @AllenKummer If you'd like to propose using a modified script, you can share it using Pastebin or other code-sharing tools. Thanks! $\endgroup$ – uhoh Dec 2 '19 at 0:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.