Видимая область Skyfield под спутником Земли

Question

Как бы я вычислил площадь под спутником Земли, чтобы я мог построить полосу земли, покрытую при прохождении спутника?

Есть ли в Скайфилде что-нибудь, что могло бы облегчить это?

Редактировать: Просто подумал, что уточнить, что я имею в виду под областью под спутником.Мне нужно изобразить максимальную площадь под спутником, которую можно наблюдать, учитывая, что Земля является сфероидом.Я знаю, как построить траекторию спутника, но теперь мне нужно построить несколько линий, чтобы представить область, видимую этим спутником, когда он летит над землей.

Answer 1

Ваша правка дала понять, что вы хотите.Видимая площадь со спутника может быть легко вычислена (когда земля рассматривается как сфера).Хороший источник для получения некоторого фона на видимой части можно найти здесь .Рассчитать видимую область, когда Земля рассматривается как сплющенный сфероид, будет намного сложнее (и, возможно, даже невозможно).Я думаю, что лучше реформировать эту часть вопроса и опубликовать ее по математике.

Если вы хотите вычислить видимую область, когда Земля рассматривается как сфера, нам нужно внести некоторые корректировки в Skyfield.С помощью спутника, загруженного с помощью API-интерфейса TLE, вы можете легко получить подпункт с положением на земле.Библиотека называет это позицией Geocentric, но на самом деле это позиция Geodetic (где земля рассматривается как сплющенный сфероид).Чтобы исправить это, нам нужно настроить subpoint класса Geocentric, чтобы использовать расчет для позиции Geocentric, а не для позиции Geodetic.Из-за ошибки и отсутствующей информации в функции reverse_terra нам также необходимо заменить эту функцию.И мы должны быть в состоянии получить радиус Земли.Это приводит к следующему:

from skyfield import api
from skyfield.positionlib import ICRF, Geocentric
from skyfield.constants import (AU_M, ERAD, DEG2RAD,
                                IERS_2010_INVERSE_EARTH_FLATTENING, tau)
from skyfield.units import Angle

from numpy import einsum, sqrt, arctan2, pi, cos, sin

def reverse_terra(xyz_au, gast, iterations=3):
    """Convert a geocentric (x,y,z) at time `t` to latitude and longitude.
    Returns a tuple of latitude, longitude, and elevation whose units
    are radians and meters.  Based on Dr. T.S. Kelso's quite helpful
    article "Orbital Coordinate Systems, Part III":
    https://www.celestrak.com/columns/v02n03/
    """
    x, y, z = xyz_au
    R = sqrt(x*x + y*y)

    lon = (arctan2(y, x) - 15 * DEG2RAD * gast - pi) % tau - pi
    lat = arctan2(z, R)

    a = ERAD / AU_M
    f = 1.0 / IERS_2010_INVERSE_EARTH_FLATTENING
    e2 = 2.0*f - f*f
    i = 0
    C = 1.0
    while i < iterations:
        i += 1
        C = 1.0 / sqrt(1.0 - e2 * (sin(lat) ** 2.0))
        lat = arctan2(z + a * C * e2 * sin(lat), R)
    elevation_m = ((R / cos(lat)) - a * C) * AU_M
    earth_R = (a*C)*AU_M
    return lat, lon, elevation_m, earth_R

def subpoint(self, iterations):
    """Return the latitude an longitude directly beneath this position.

    Returns a :class:`~skyfield.toposlib.Topos` whose ``longitude``
    and ``latitude`` are those of the point on the Earth's surface
    directly beneath this position (according to the center of the
    earth), and whose ``elevation`` is the height of this position
    above the Earth's center.
    """
    if self.center != 399:  # TODO: should an __init__() check this?
        raise ValueError("you can only ask for the geographic subpoint"
                            " of a position measured from Earth's center")
    t = self.t
    xyz_au = einsum('ij...,j...->i...', t.M, self.position.au)
    lat, lon, elevation_m, self.earth_R = reverse_terra(xyz_au, t.gast, iterations)

    from skyfield.toposlib import Topos
    return Topos(latitude=Angle(radians=lat),
                    longitude=Angle(radians=lon),
                    elevation_m=elevation_m)

def earth_radius(self):
    return self.earth_R

def satellite_visiable_area(earth_radius, satellite_elevation):
    """Returns the visible area from a satellite in square meters.

    Formula is in the form is 2piR^2h/R+h where:
        R = earth radius
        h = satellite elevation from center of earth
    """
    return ((2 * pi * ( earth_radius ** 2 ) * 
            ( earth_radius + satellite_elevation)) /
            (earth_radius + earth_radius + satellite_elevation))


stations_url = 'http://celestrak.com/NORAD/elements/stations.txt'
satellites = api.load.tle(stations_url)
satellite = satellites['ISS (ZARYA)']
print(satellite)

ts = api.load.timescale()
t = ts.now()

geocentric = satellite.at(t)
geocentric.subpoint = subpoint.__get__(geocentric, Geocentric)
geocentric.earth_radius = earth_radius.__get__(geocentric, Geocentric)

geodetic_sub = geocentric.subpoint(3)

print('Geodetic latitude:', geodetic_sub.latitude)
print('Geodetic longitude:', geodetic_sub.longitude)
print('Geodetic elevation (m)', int(geodetic_sub.elevation.m))
print('Geodetic earth radius (m)', int(geocentric.earth_radius()))

geocentric_sub = geocentric.subpoint(0)
print('Geocentric latitude:', geocentric_sub.latitude)
print('Geocentric longitude:', geocentric_sub.longitude)
print('Geocentric elevation (m)', int(geocentric_sub.elevation.m))
print('Geocentric earth radius (m)', int(geocentric.earth_radius()))
print('Visible area (m^2)', satellite_visiable_area(geocentric.earth_radius(), 
                                                    geocentric_sub.elevation.m))