Try out a rival approach to is_sunlit()
This commit is contained in:
Brandon Rhodes
parent
dac6d847d1
commit
06ee83aa66
|
|
@ -394,53 +394,50 @@ See :doc:`positions` to learn more about these possibilities.
|
|||
Finding when a satellite will be illuminated
|
||||
--------------------------------------------
|
||||
|
||||
It is sometimes important to understand if the sun is
|
||||
currently illuminating the satellite.
|
||||
For instance, if it's after sundown in your location,
|
||||
and the satellite is still illuminated by the sun,
|
||||
you might be able to see that sunlight reflected with the naked eye.
|
||||
Additionally, this can be helpful when understanding
|
||||
satellite power and thermal cycles as it goes in and
|
||||
out of eclipse.
|
||||
A satellite is generally only visible to a ground observer
|
||||
when there is still sunlight up at its altitude.
|
||||
The satellite will visually disappear
|
||||
when it enters the Earth’s shadow
|
||||
and reappear when it comes out of eclipse.
|
||||
If you are planning to observe a satellite visually,
|
||||
rather than with radar or radio,
|
||||
you will want to know which satellite passes are in sunlight.
|
||||
Knowing a satellite’s sunlit periods
|
||||
is also helpful when modeling satellite power and thermal cycles
|
||||
as it goes in and out of eclipse.
|
||||
|
||||
Skyfield provides a simple geometric estimation for this
|
||||
with the :meth:`~skyfield.sgp4lib.EarthSatellite.is_sunlit()` method,
|
||||
which provides a boolean response for each corresponding time provided.
|
||||
This requires knowledge of the earth and sun,
|
||||
and will require providing ephemeris data.
|
||||
You can check if the satellite will be illuminated at whatever
|
||||
times you desire (when the TLE is still accurate):
|
||||
Skyfield provides a simple geometric estimate for this
|
||||
through the :meth:`~skyfield.sgp4lib.EarthSatellite.is_sunlit()` method.
|
||||
Given an ephemeris with which it can compute the Sun’s position,
|
||||
it will return ``True`` when the satellite is in sunlight
|
||||
and ``False`` otherwise.
|
||||
|
||||
.. testcode::
|
||||
|
||||
from skyfield.api import load
|
||||
satellite = by_name['ISS (ZARYA)']
|
||||
de421 = load('de421.bsp')
|
||||
eph = load('de421.bsp')
|
||||
satellite = by_name['ISS (ZARYA)']
|
||||
|
||||
# Define the times you are interested in
|
||||
minutes = range(0, 90, 10)
|
||||
times = ts.utc(2020, 4, 29, 0, minutes)
|
||||
two_hours = ts.utc(2014, 1, 20, 0, range(0, 120, 20))
|
||||
|
||||
# Calculate the sunlit vector
|
||||
sunlit = satellite.is_sunlit(ephemeris=de421, times=times)
|
||||
for idx, time in enumerate(times):
|
||||
print(time.utc_jpl(), "{} is sunlit: {}".format(satellite.name,sunlit[idx]))
|
||||
sunlit1 = satellite.is_sunlit(ephemeris=eph, times=two_hours)
|
||||
sunlit2 = satellite.at(two_hours).is_sunlit(eph)
|
||||
|
||||
This produces a `sunlit` vector of booleans you can reference alongside your times:
|
||||
for ti, sunlit1_i, sunlit2_i in zip(two_hours, sunlit1, sunlit2):
|
||||
print('{} {} is in {} {}'.format(
|
||||
ti.utc_strftime('%Y-%m-%d %H:%M'),
|
||||
satellite.name,
|
||||
'sunlight' if sunlit1_i else 'shadow',
|
||||
'sunlight' if sunlit2_i else 'shadow',
|
||||
))
|
||||
|
||||
.. testoutput::
|
||||
|
||||
A.D. 2020-Apr-29 00:00:00.0000 UT ISS (ZARYA) is sunlit: True
|
||||
A.D. 2020-Apr-29 00:10:00.0000 UT ISS (ZARYA) is sunlit: True
|
||||
A.D. 2020-Apr-29 00:20:00.0000 UT ISS (ZARYA) is sunlit: False
|
||||
A.D. 2020-Apr-29 00:30:00.0000 UT ISS (ZARYA) is sunlit: False
|
||||
A.D. 2020-Apr-29 00:40:00.0000 UT ISS (ZARYA) is sunlit: False
|
||||
A.D. 2020-Apr-29 00:50:00.0000 UT ISS (ZARYA) is sunlit: False
|
||||
A.D. 2020-Apr-29 01:00:00.0000 UT ISS (ZARYA) is sunlit: True
|
||||
A.D. 2020-Apr-29 01:10:00.0000 UT ISS (ZARYA) is sunlit: True
|
||||
A.D. 2020-Apr-29 01:20:00.0000 UT ISS (ZARYA) is sunlit: True
|
||||
|
||||
You can see when it is sunlit at the given times!
|
||||
2014-01-20 00:00 ISS (ZARYA) is in sunlight sunlight
|
||||
2014-01-20 00:20 ISS (ZARYA) is in sunlight sunlight
|
||||
2014-01-20 00:40 ISS (ZARYA) is in shadow shadow
|
||||
2014-01-20 01:00 ISS (ZARYA) is in shadow shadow
|
||||
2014-01-20 01:20 ISS (ZARYA) is in sunlight sunlight
|
||||
2014-01-20 01:40 ISS (ZARYA) is in sunlight sunlight
|
||||
|
||||
Avoid calling the observe method
|
||||
--------------------------------
|
||||
|
|
|
|||
|
|
@ -3,27 +3,24 @@
|
|||
from numpy import nan, where
|
||||
from .functions import length_of
|
||||
|
||||
def intersect_line_and_sphere(line_endpoint, sphere_center, sphere_radius):
|
||||
def intersect_line_and_sphere(endpoint, center, radius):
|
||||
"""Compute distance to intersections of a line and a sphere.
|
||||
|
||||
Given a line that passes through the origin (0,0,0) and an (x,y,z)
|
||||
``line_endpoint``, and a sphere with the (x,y,z) ``sphere_center``
|
||||
and scalar ``sphere_radius``, return the distance from the origin to
|
||||
their two intersections. If the line is tangent to the sphere, the
|
||||
two intersections will be at the same distance. If the line does
|
||||
not intersect the sphere, two ``nan`` values will be returned.
|
||||
Given a line through the origin (0,0,0) and an (x,y,z) ``endpoint``,
|
||||
and a sphere with the (x,y,z) ``center`` and scalar ``radius``,
|
||||
return the distance from the origin to their two intersections.
|
||||
|
||||
If the line is tangent to the sphere, the two intersections will be
|
||||
at the same distance. If the line does not intersect the sphere,
|
||||
two ``nan`` values will be returned.
|
||||
|
||||
"""
|
||||
# See http://paulbourke.net/geometry/circlesphere/index.html#linesphere
|
||||
# Names "b" and "c" designate the familiar values from the quadratic
|
||||
# formula; happily, a = 1 because we use a unit vector for the line.
|
||||
|
||||
unit_vector = line_endpoint / length_of(line_endpoint)
|
||||
minus_b = 2.0 * (unit_vector * sphere_center).sum(axis=0)
|
||||
r_squared = sphere_radius * sphere_radius
|
||||
c = (sphere_center * sphere_center).sum(axis=0) - r_squared
|
||||
minus_b = 2.0 * (endpoint / length_of(endpoint) * center).sum(axis=0)
|
||||
c = (center * center).sum(axis=0) - radius * radius
|
||||
discriminant = minus_b * minus_b - 4 * c
|
||||
exponent = where(discriminant < 0, nan, 0.5)
|
||||
dsqrt = discriminant ** exponent
|
||||
print(discriminant, (minus_b - dsqrt) / 2.0, (minus_b + dsqrt) / 2.0)
|
||||
dsqrt = discriminant ** where(discriminant < 0, nan, 0.5) # avoid sqrt(<0)
|
||||
return (minus_b - dsqrt) / 2.0, (minus_b + dsqrt) / 2.0
|
||||
|
|
|
|||
|
|
@ -1,10 +1,11 @@
|
|||
# -*- coding: utf-8 -*-
|
||||
"""Classes representing different kinds of astronomical position."""
|
||||
|
||||
from numpy import arccos, array, clip, einsum, exp, full, nan
|
||||
from .constants import ANGVEL, DAY_S, DEG2RAD, RAD2DEG, tau
|
||||
from numpy import arccos, array, clip, einsum, exp, full, nan, nan_to_num
|
||||
from .constants import ANGVEL, ERAD, DAY_S, DEG2RAD, RAD2DEG, tau
|
||||
from .data.spice import inertial_frames
|
||||
from .earthlib import compute_limb_angle, refract, reverse_terra
|
||||
from .geometry import intersect_line_and_sphere
|
||||
from .functions import (
|
||||
_mxv, _mxm, dots, from_spherical, length_of, rot_x, rot_z,
|
||||
to_spherical, _to_array,
|
||||
|
|
@ -677,6 +678,14 @@ class Geocentric(ICRF):
|
|||
longitude=Angle(radians=lon),
|
||||
elevation_m=elevation_m)
|
||||
|
||||
def is_sunlit(self, ephemeris):
|
||||
"""Return whether a position in Earth orbit is in sunlight."""
|
||||
sun_m = (ephemeris['sun'] - ephemeris['earth']).at(self.t).position.m
|
||||
earth_m = - self.position.m
|
||||
near, far = intersect_line_and_sphere(sun_m + earth_m, earth_m, ERAD)
|
||||
return nan_to_num(far) <= 0
|
||||
# BUT: what about normal satellite positions relative to an observer?
|
||||
|
||||
|
||||
def _to_altaz(position_au, observer_data, temperature_C, pressure_mbar):
|
||||
"""Compute (alt, az, distance) relative to the observer's horizon.
|
||||
|
|
|
|||
Loading…
Reference in New Issue