Wie kann ich mit pyEphem oder Skyfield (bevorzugt) berechnen, ob ein Satellit derzeit im Sonnenlicht oder in einer Sonnenfinsternis ist?

Als neuer Benutzer kann ich die ursprüngliche Antwort nicht kommentieren oder optimieren, also versuche ich meine eigene. Frohes Update für alle Empfehlungen, die auftauchen.

Die Antwort von @uhoh ist nah dran, aber ein paar Dinge sind zu beachten.

Die TLE sagt uns, dass die Anzahl der Umlaufbahnen pro Tag 15,50995519 beträgt (Zeile 2, Spalten 53–63).

Wir erwarten also, dass unser Satellit alle ~93 Minuten einmal umkreist. Über einen Orbit erwarten wir, dass der Satellit jeweils einmal in und aus der Sonnenfinsternis übergeht.

Die Antwort von @uhoh zeigt zwei solche Übergänge in dieser Zeit. Was fehlt, ist die Überlegung, dass sich der Satellit zweimal im Kegelwinkel des Erdrandes befindet. Sowohl bei Sonnenfinsternis als auch direkt vor der Erde.

Um auf der Antwort von @uhoh aufzubauen, können wir die Entfernung Sonne-Erde und Sonne-Sat berücksichtigen. Der Sat befindet sich in der Sonne, wenn er sich außerhalb des Extremitätenwinkels befindet ODER wenn er näher an der Sonne als an der Erde ist.

Leicht modifizierter Code (Danke an @uhoh, dass du uns gezeigt hast, wie man den Sonne-Erde-Abstand berechnet!)

Ergebnisse:

Um das Ergebnis zu testen, kann man das erste Auftreten des Übergangs True-False und False-True finden über:

sunlit[::].index(False)

und weiter schneiden:

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

Auf diese Weise finden wir Übergänge bei 141 (Wahr-Falsch – Betreten der Sonnenfinsternis), dann 247 (Falsch-Wahr – Beenden der Sonnenfinsternis) und 681 (Wahr-Falsch – Betreten der nächsten Sonnenfinsternis).

247 + 681 / 10 [jeder Zeitschritt ist 1/10 einer Minute] = 92,8 Minuten zwischen Eclipse-Einträgen.

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

MEGA-Hinweis: Meine Antwort ist völlig falsch und sollte mit allen erforderlichen Mitteln nicht akzeptiert werden . Glücklicherweise erklärt die Antwort von @AllenKummer warum und gibt ein überarbeitetes Skript. Dies wurde gegen 7 Uhr Ortszeit veröffentlicht und muss vor dem Kaffee gewesen sein; Der hier gezeigte Zeitraum beträgt nur 40 Minuten, da es falsch ist, wäre ein Kinderspiel gewesen, was in der Tat der Zustand ist, in dem ich mich oft befinde, vor dem Kaffee.

Hinweis: Ich habe zuerst in Python und nicht in MathJax geantwortet , aber ich werde die entsprechenden Gleichungen später heute hinzufügen.

Am einfachsten ist es, den Öffnungswinkel von der Sonne aus gesehen oder den Erde-Sonne-Sat-Winkel zu berechnen. Verwenden Sie ein hausgemachtes Punktprodukt normalisierter Vektoren, um die Winkel zu erhalten, und berechnen Sie dann den Erde-Sonne-Sat-Winkel.

Wenn der Erde-Sonne-Sat-Winkel größer ist als der Erde-Sonne-Glied-Winkel, dann ist Sat im Tageslicht. Wenn nicht, dann liegt es im Schatten der Erde.

Der Übergang von „Tag“ zu „Nacht“ (Finsternis) beträgt nur wenige Sekunden.

Um ein 3D-Diagramm der Umlaufbahn in Python zu erstellen, lesen Sie diese Antwort (Python-2) auf Wie kann ich die Umlaufbahn eines Satelliten in 3D von einem TLE mit Python und Skyfield zeichnen? und diese Python-3-Anpassung

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