The First Day of a Satellite
The first day of a satellite in orbit requires a lot of planning (even though the plans go out of the window rather quickly in many cases). The critical inputs to this planning are: - When will we be able to talk to the satellite (groundstation contact times) - Where should I point my groundstation antenna (groundstation antenna azimuth-elevation angles) - When will the satellite see the sun and generate power (satellite illumination times)
This how-to guide shows how these bits of crucial information can be generated with satkit. The steps are: 1. Set up the orbit 2. Set up the groundstation 3. Set up the search interval and the propagator 4. Run the analyses
However, as always, the very first thing to do is to initialise satkit and Orekit.
[ ]:
from pathlib import Path
# If satkit import fails, try to locate the module
import os
try:
import satkit
except ModuleNotFoundError:
os.chdir(os.path.join("..", ".."))
os.getcwd()
from satkit import init_satkit, u
init_satkit(Path("data", "orekit-data", "orekit-data-reference.zip"), Path(".."))
# Orekit / satkit init complete
Then the user groundstation location is initialised. The QuantityInit class provides the geodetic_point method to initialise the groundstation location safely with units.
[ ]:
from satkit.utils.quantity_init import QuantityInit
# Set your groundstation location (geodetic coordinates)
longitude = 45.0 * u.deg
latitude = 40 * u.deg
altitude = 0 * u.m
gnd_station = QuantityInit.geodetic_point(latitude, longitude, altitude)
The next step is to set up the search interval. The search start start and search duration parameters are defined by the target analysis duration and are to be set by the user.
[ ]:
from org.orekit.time import TimeScalesFactory
from satkit.time.time import AbsoluteDateExt
from satkit.time.timeinterval import TimeInterval
# Shorthand for UTC
utc = TimeScalesFactory.getUTC()
# Set your search interval
search_start = AbsoluteDateExt("2020-06-10T00:00:00.000", utc)
duration = 1 * u.day
search_interval = TimeInterval.from_duration(search_start, duration)
Finally it is time to set up the orbit and the propagator.
Usually, the orbit (in the form of TLE) is given in Spacetrack or Celestrak websites, even though for small satellites launched in large batches, they can’t easily identify which satellite is which, at least at the very beginning. Nevertheless, UHF comms is very forgiving and in many cases even an approximate TLE would work. In the first part of this example, we assume that a reasonably good TLE is known. The TLE goes with the SGP4 analytical propagator, and very little set up is necessary.
The other alternative is to start with the cartesian coordinates (for example, separation coordinates provided by the launcher). In this case, the full set of satellite parameters (orbit, physical properties etc.) as well as the numerical propagator with a proper force model should be set up. Clearly, this is much more involved than the TLE and SGP4 combination. A detailed gravity model and a simple drag model (NRLMSISE00 atmosphere and a spherical satellite model) is shown here. For a propagation for a few days, the solar or lunar gravity is probably not needed. Solar radiation pressure is also too small to have a real effect for most smallsats or cubesats within this timeframe and for the purposes of illumination and pass computation.
The using_tle flag below chooses between the TLE based operations and the cartesian initial conditions. The end result is the propagator object that is used for the analyses later.
[ ]:
from org.orekit.propagation.analytical.tle import TLE, TLEPropagator
from org.orekit.propagation.numerical import NumericalPropagator
from org.hipparchus.ode.nonstiff import DormandPrince853Integrator
from org.orekit.propagation import SpacecraftState
from org.orekit.orbits import CartesianOrbit, OrbitType
from org.orekit.forces.gravity.potential import GravityFieldFactory
from org.orekit.forces.gravity import HolmesFeatherstoneAttractionModel
from orekit import JArray_double
from org.orekit.models.earth.atmosphere.data import CssiSpaceWeatherData
from org.orekit.models.earth.atmosphere import NRLMSISE00
from org.orekit.forces.drag import IsotropicDrag
from org.orekit.forces.drag import DragForce
from org.orekit.bodies import CelestialBodyFactory, OneAxisEllipsoid
from org.orekit.utils import IERSConventions, TimeStampedPVCoordinates, PVCoordinatesProvider, Constants
from org.orekit.frames import FramesFactory
from org.hipparchus.geometry.euclidean.threed import Vector3D
# True if a TLE is used for the orbit, False if a cartesian initial condition is used
using_tle = True
if using_tle:
# Insert your TLE here
line1 = "1 40697U 15028A 20164.50828565 .00000010 00000-0 20594-4 0 9999"
line2 = "2 40697 98.5692 238.8182 0001206 86.9662 273.1664 14.30818200259759"
# Init TLE object
tle = TLE(line1, line2)
# Set up the SGP4 propagator
propagator = TLEPropagator.selectExtrapolator(tle)
else:
# Init useful frames
itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, True)
gcrf = FramesFactory.getGCRF()
teme = FramesFactory.getTEME()
earth = OneAxisEllipsoid(
Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
Constants.WGS84_EARTH_FLATTENING,
itrf,
)
sun_coords = PVCoordinatesProvider.cast_(CelestialBodyFactory.getSun())
# User-defined initial state with cartesian coordinates
# (e.g., spacecraft separation vector)
initialDate = AbsoluteDateExt("2020-06-12T12:11:55.88016Z", utc)
position = Vector3D(-3712324.066353917, -6134153.57973047, 64.00521472066569)
velocity = Vector3D(-945.60658471178, 582.3975490543295, 7374.596184023312)
pos_vel_time = TimeStampedPVCoordinates(initialDate, position, velocity)
initial_orbit = CartesianOrbit(pos_vel_time, gcrf, Constants.WGS84_EARTH_MU)
# User-defined spacecraft mass [kg]
satellite_mass = 10.0
# initial state
initial_state = SpacecraftState(initial_orbit, satellite_mass)
# Set up the numerical propagator
min_step = 0.001
max_step = 1000.0
init_step = 60.0
# spatial tolerance (meters)
position_tolerance = 1.0
tolerances = NumericalPropagator.tolerances(position_tolerance,
initial_orbit,
initial_orbit.getType())
integrator = DormandPrince853Integrator(min_step, max_step,
# Double array of doubles needs to be cast in Python
JArray_double.cast_(tolerances[0]),
JArray_double.cast_(tolerances[1]))
integrator.setInitialStepSize(init_step)
propagator = NumericalPropagator(integrator)
propagator.setOrbitType(OrbitType.CARTESIAN)
propagator.setInitialState(initial_state)
# Gravity Model for the propagation (degree:8, order: 8)
gravityProvider = GravityFieldFactory.getNormalizedProvider(8, 8)
propagator.addForceModel(HolmesFeatherstoneAttractionModel(itrf, gravityProvider))
# Drag Model for the propagation
# User-defined satellite properties
cross_section = 0.1 * 0.1 # cross-section area in m^2
cd = 2.0 # drag coefficient
# Initialise atmosphere
cswl = CssiSpaceWeatherData("SpaceWeather-All-v1.2.txt")
atmosphere = NRLMSISE00(cswl, sun_coords, earth)
isotropic_drag = IsotropicDrag(cross_section, cd)
# Add the drag model to the forces
propagator.addForceModel(DragForce(atmosphere, isotropic_drag))
# Finally, run the propagation
end_state = propagator.propagate(search_interval.start, search_interval.end)
print(f"end state: {end_state}")
After all the preparations, the analyses can finally be run. The first analysis is to generate the “satellite illumination events”, showing when the satellite sees the sunlight, which is critical for thermal analyses and power generation.
This is also useful to relate the received telemetry to the physical events. For example, the telemetry should show the activation of some of the solar arrays around illumination start time - the absence of any activity may therefore indicate an attitude determination or control failure, erroneous orbit definition or a solar array failure.
[ ]:
from satkit.events.eventfinders import sat_illum_finder
# ----- find sat illumination times -----
illum_times = sat_illum_finder(search_interval, propagator, use_total_eclipse=True)
print("Satellite Illumination Intervals:")
print("----------------------------------------------------------------------------")
print(illum_times)
The second analysis is to find the passes above a groundstation. This should give us a detailed list of all pass intervals, maximum elevation value (may determine the quality of the communications as well help plan the tasks for the pass) and azimuth-elevation-range as well as range-rate values, to plan the communications (antenna pointing and doppler corrections).
[ ]:
from satkit.events.gnd_access import GroundPassList
# ----- find passes -----
# minimum elevation definition
min_elev = 5 * u.deg
# stepsize for the azimuth
az_el_timestep = 60 * u.s
# find the passes
passes = GroundPassList(search_interval,
gnd_station,
min_elev,
propagator,
planet=None,
refraction_model=None,
az_el_timestep=az_el_timestep)
print(f"Groundstation Contact Intervals (above {min_elev:~P} elevation):")
print("----------------------------------------------------------------------------")
for gnd_pass in passes.pass_list:
print(f"{gnd_pass.pass_interval} -- max elev: {gnd_pass.max_elev_time} ({gnd_pass.max_elev_aer.el.to('deg'):~P})")
Finally, the azimuth-elevation list of the first pass is given as an example. Note that the stepsize was given in the previous step.
[ ]:
# ----- find az-el-range -----
print()
print("Groundstation contact az-el list (first pass only, as a sample):")
print("----------------------------------------------------------------------------")
for az_el in passes.pass_list[0].az_el_list:
print(f"t: {az_el.t} az: {az_el.az.to('deg'):~P} el: {az_el.el.to('deg'):~P}")