,
Daniel Ayala [aut] ,
David Ruiz [aut] ,
Lara Selles Vidal [aut] | Title: | Computation of Satellite Position |
|---|---|
| Description: | Provides basic functionalities to calculate the position of satellites given a known state vector. The package includes implementations of the SGP4 and SDP4 simplified perturbation models to propagate orbital state vectors, as well as utilities to read TLE files and convert coordinates between different frames of reference. Several of the functionalities of the package (including the high-precision numerical orbit propagator) require the coefficients and data included in the 'asteRiskData' package, available in a 'drat' repository. To install this data package, run 'install.packages("asteRiskData", repos="https://rafael-ayala.github.io/drat/")'. Felix R. Hoots, Ronald L. Roehrich and T.S. Kelso (1988) <https://celestrak.org/NORAD/documentation/spacetrk.pdf>. David Vallado, Paul Crawford, Richard Hujsak and T.S. Kelso (2012) <doi:10.2514/6.2006-6753>. Felix R. Hoots, Paul W. Schumacher Jr. and Robert A. Glover (2014) <doi:10.2514/1.9161>. |
| Authors: | Rafael Ayala [aut, cre] ,
Daniel Ayala [aut] ,
David Ruiz [aut] ,
Lara Selles Vidal [aut] |
| Maintainer: | Rafael Ayala <[email protected]> |
| License: | GPL-3 |
| Version: | 1.4.3.9999 |
| Built: | 2024-11-11 04:46:46 UTC |
| Source: | https://github.com/rafael-ayala/asterisk |
The horizontal coordinate system, also called azimuth-elevation system, uses the local horizon of an observer as its fundamental plane. In it, a given point is defined by 2 main angles: azimuth and elevation. Azimuth defines the angle of the point around the horizon in the X-Y plane, measured from the true North and usually increasing towards the East. Elevation is the angle between the object and the X-Y plane. Finally, the range defines the distance between the observer and the point.
This function calculates the azimuth, elevation and range given the coordinates of an observed satellite and of an observer. Both sets of coordinates must be provided as Cartesian geocentric coordinates in ITRF.
calculateRazel(geocentricObserver, geocentricSatellite, degreesOutput=TRUE)calculateRazel(geocentricObserver, geocentricSatellite, degreesOutput=TRUE)
geocentricObserver |
Vector with the X, Y and Z components of the position of the observer in ITRF frame. |
geocentricSatellite |
Vector with the X, Y and Z components of the position of the satellite in ITRF frame. |
degreesOutput |
Logical indicating if the output should be in sexagesimal
degrees. If |
A vector with three elements, corresponding to the azimuth and elevation in degrees (or radians if specified) and the range in the same unit as the provided Cartesian coordinates.
https://gssc.esa.int/navipedia/index.php/Transformations_between_ECEF_and_ENU_coordinates
# The following were the coordinates of Italsat-2 in ITRF the 27th of June, 2006 # at 00:58:29.34 UTC, in meters. italsat_ITRF <- c(-37325542.8, 19152438.3, 138384.5) # Let us calculate its azimuth, elevation and range for an observer from Tokyo. # The latitude and longitude of the city are 35.6762 degrees North, 139.6503 # degrees East. Let's assume an observer placed at sea level (0 m) # We first convert these coordinates to ITRF: observer_ITRF <- LATLONtoITRF(c(35.6762, 139.6503, 0), degreesInput=TRUE) # We can now calculate the azimuth and elevation: razel <- calculateRazel(observer_ITRF, italsat_ITRF, degreesOutput=TRUE) razel[1] # Azimuth razel[2] # Elevation# The following were the coordinates of Italsat-2 in ITRF the 27th of June, 2006 # at 00:58:29.34 UTC, in meters. italsat_ITRF <- c(-37325542.8, 19152438.3, 138384.5) # Let us calculate its azimuth, elevation and range for an observer from Tokyo. # The latitude and longitude of the city are 35.6762 degrees North, 139.6503 # degrees East. Let's assume an observer placed at sea level (0 m) # We first convert these coordinates to ITRF: observer_ITRF <- LATLONtoITRF(c(35.6762, 139.6503, 0), degreesInput=TRUE) # We can now calculate the azimuth and elevation: razel <- calculateRazel(observer_ITRF, italsat_ITRF, degreesOutput=TRUE) razel[1] # Azimuth razel[2] # Elevation
The Julian Date (JD) of a given date and time is the number of days since noon of Monday 1st of January 4713 BC, including a fractional part. Modified Julian Date (MJD) are instead the number of days since 00:00 of November 17th, 1858. The difference JD and MJD for a given instant is always 2400000.5, which is the JD of the reference time for MJD.
This function calculates the MJD of a date and time, provided as a date-time character string in UTC time. The output refers by default to the MJD in UTC, but different time systems can be chosen: UTC (Coordinated Universal Time), UT1 (Universal Time), TT (Terrestrial Time) and TDB (Barycentric Dynamical Time).
dateTimeToMJD(dateTime, timeSystem="UTC")dateTimeToMJD(dateTime, timeSystem="UTC")
dateTime |
Date-time string with the date and time in UTC corresponding to the provided geodetic coordinates. |
timeSystem |
Time system into which the MJD should be calculated. Should be one from "UTC" (Coordinated Universal Time; default), "UT1" (Universal Time), "TT" (Terrestrial Time) and "TDB" (Barycentric Dynamical Time). |
The MJD for the specified date and time in the chosen time system.
https://gssc.esa.int/navipedia/index.php/Julian_Date https://gssc.esa.int/navipedia/index.php/Transformations_between_Time_Systems
if(requireNamespace("asteRiskData", quietly = TRUE)) { # Let's calculate the MJD of the 12th of June, 2000 at 10:00:00 UTC time, in UTC MJD_UTC <- dateTimeToMJD("2000-06-12 10:00:00", timeSystem = "UTC") # Let's now calculate the MJD for the same instant in TDB: MJD_TDB <- dateTimeToMJD("2000-06-12 10:00:00", timeSystem = "TDB") # We can now calculate the difference in seconds, which matches the difference # between UTC and TDB for that day: (MJD_UTC - MJD_TDB) * 86400 }if(requireNamespace("asteRiskData", quietly = TRUE)) { # Let's calculate the MJD of the 12th of June, 2000 at 10:00:00 UTC time, in UTC MJD_UTC <- dateTimeToMJD("2000-06-12 10:00:00", timeSystem = "UTC") # Let's now calculate the MJD for the same instant in TDB: MJD_TDB <- dateTimeToMJD("2000-06-12 10:00:00", timeSystem = "TDB") # We can now calculate the difference in seconds, which matches the difference # between UTC and TDB for that day: (MJD_UTC - MJD_TDB) * 86400 }
This function converts an angle in degrees to radians.
deg2rad(degrees)deg2rad(degrees)
degrees |
Value of the angle in degrees. |
The corresponding value of the angle in radians.
deg2rad(180)deg2rad(180)
Keplerian orbital elements are a set of six parameters used to described the orbits of celestial objects, including satellites. While satellites do not follow a perfectly Keplerian orbit, their state at any point can be defined by the orbital parameters that they would have if they were located at the same position with the same velocity following a perfectly Keplerian orbit (i.e., if perturbations were absent). These are called osculating orbital elements.
Keplerian orbital elements can be unequivocally determined from a satellite if its position and velocity are known. This function calculates orbital elements from the position and velocity of a satellite in an ECI (Earth-centered inertial) frame of reference. The elements (such as the equatorial plane) with respect to which the resulting orbital elements will be defined are the same as those used for the ECI frame of reference. The function calculates the six standard orbital elements, plus some alternative elements useful for the characterization of special orbits, such as circular ones or orbits with no inclination.
ECItoKOE(position_ECI, velocity_ECI)ECItoKOE(position_ECI, velocity_ECI)
position_ECI |
Vector with the X, Y and Z components of the position of an object in an ECI frame, in m. |
velocity_ECI |
Vector with the X, Y and Z components of the velocity of an object in an ECI frame, in m/s. |
A list with the following standard and alternative orbital elements:
semiMajorAxis |
Semi-major axis of orbital ellipse in meters. |
eccentricity |
Numerical eccentricity of the orbit. Eccentricity measures how much the orbit deviates from being circular. |
inclination |
Inclination of the orbital plane in radians. Inclination is the angle between the orbital plane and the equator. |
meanAnomaly |
Mean anomaly of the orbit in radians. Mean anomaly indicates where the satellite is along its orbital path, and is defined as the angle between the direction of the perigee and the hypothetical point where the object would be if it was moving in a circular orbit with the same period as its true orbit after the same amount of time since it last crossed the perigee had ellapsed. |
argumentPerigee |
Argument of perigee in radians. This is the angle between the direction of the ascending node and the direction of the perigee (the point of the orbit at which the object is closest to the Earth). |
longitudeAscendingNode |
Longitude of the ascending node (also called right ascension of the ascending node) in radians. This is the angle between the direction of the ascending node (the point where thesatellite crosses the equatorial plane moving north) and the direction of the First Point of Aries (which indicates the location of the vernal equinox). |
trueAnomaly |
True anomaly of the orbit in radians. Unlike mean anomaly, true anomaly is the angle between the direction of the perigee and the actual position of the satellite. |
argumentLatitude |
Argument of latitude of the orbit in radians. Defined as the angle between the equator and the position of the satellite. It is useful to define the position of satellites in circular orbits, where the argument of perigee and true anomaly are not well defined. |
longitudePerigee |
Longitude of perigee of the orbit in radians. Defined as the angle between the vernal equinox and the perigee. It is useful for cases of orbits with 0 inclination, where the longitude of the ascending node and the argument of perigee are not well defined. |
trueLongitude |
Longitude of perigee of the orbit in radians. Defined as the angle between the vernal equinox and the position of the satellite. It is useful for cases of circular orbits with 0 inclination, where the longitude of the ascending node, the argument of perigee and true anomaly are not well defined. |
https://www.gsc-europa.eu/system-service-status/orbital-and-technical-parameters https://celestrak.org/columns/v02n01/ https://www.faa.gov/about/office_org/headquarters_offices/avs/offices/aam/cami/library/online_libraries/aerospace_medicine/tutorial/media/iii.4.1.4_describing_orbits.pdf
# The following were the position and velocity of satellite MOLNIYA 1-83 # the 25th of June, 2006 at 00:33:43 UTC in the GCRF frame (in m and m/s). position_GCRF <- c(-14471729.582, -4677558.558, 9369.461) velocity_GCRF <- c(-3251.691, -3276.008, 4009.228) # Let's calculate the orbital elements of the satellite at that time orbital_elements <- ECItoKOE(position_GCRF, velocity_GCRF)# The following were the position and velocity of satellite MOLNIYA 1-83 # the 25th of June, 2006 at 00:33:43 UTC in the GCRF frame (in m and m/s). position_GCRF <- c(-14471729.582, -4677558.558, 9369.461) velocity_GCRF <- c(-3251.691, -3276.008, 4009.228) # Let's calculate the orbital elements of the satellite at that time orbital_elements <- ECItoKOE(position_GCRF, velocity_GCRF)
SPK (Spacecraft and Planet Kernel) is a binary file format developed by NAIF to store ephemerides (trajectory) of celestial bodies and spacecraft in the Solar System. A detailed description of the SPK file format can be found in NAIF's documentation (https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/req/spk.html).
Each SPK file contains several segments, with each segment comprising a summary (or descriptor), a name and an array of double precision elements. Each segment is conceptually equivalent to an array in the context of generic DAF files. There are several types of SPK segments defined by NAIF, each identified by an SPK type code currently ranging from 1 to 21 (some intermediate values are not used or not available for general public use). Each segment type provides ephemerides information in a different way. Note that the segments stored in a single SPK file can be of different types. A detailed description of the organization of the arrays for each SPK type can be found at https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/req/spk.html
SPK files can be read with readBinSPK. This function allows the evaluation of one segment of an SPK file to any target times. The evaluated segment should be one of the elements in the list of segments returned by calling readBinSPK on an SPK file.
Evaluation at target times of the SPK segment produces state vectors comprising X, Y and Z components for position and velocity for the target body. The target body, center body and frame of reference are those indicated in the corresponding summary of the SPK segment.
The specific algorithm through which an SPK segment is evaluated depends on the type of SPK segment. For a detailed description, see the documentation of readBinSPK or NAIF's documentation.
The output values for all types of segments have been verified to match those
of CSPICE to a precision of sqrt(.Machine$double.eps).
Note that, in addition to position and velocity components, evaluateSPK
will also produce X, Y and Z components of acceleration for segments of types
1, 2, 3, 8, 9, 12, 13, 14, 18, 19, 20 and 21. This is possible because all of
these segments provide some sort of interpolation polynomial coefficients, which
can be differentiated to obtain acceleration values. However, note that said
coefficients are, in principle, not intended to be used to calculate acceleration
values, and SPICE does not return these. Therefore, proceed with caution if
using acceleration values obtained in this way. In the specific case of SPK segments
of types 1 and 21, the provided coefficents are for an interpolation polynomial
for acceleration, and therefore are probably the most reliable. In particular,
for target times exactly matching the reference epoch of any of the data points
included in the segment, the acceleration value should match the original value
used to fit the interpolation polynomial.
evaluateSPKSegment(segment, targetEpochs)evaluateSPKSegment(segment, targetEpochs)
segment |
Single SPK segment to be evaluated. The segment should have the same structure as that of segments returned by the readBinSPK function. Note said function returns all segments contained in the read SPK file; only one of those should be selected and provided here. |
targetEpochs |
Numeric vector indicating the target epochs at which the segment should be evaluated. The epochs should be provided in TDB seconds since J2000, also known as ephemeris time in SPICE. Note that all target epochs should be larger than the start epoch of the segment, and smaller than the end epoch of the segment, both of which are specified in the corresponding segment summary. Segments of type 1 and 21 are an exception. There are segments of these types where the final record covers epochs up to a time larger than the global end epoch of the corresponding segment; in cases where epochs larger than the global end epoch of the segment, but smaller than the end epoch of the last data record, these epochs will also be evaluated. |
A matrix with 7 (for SPK segments of types 5, 10 and 15) or 10 (for all other types) columns. Each row in the matrix represents a target time of evaluation. Column 1 provides the epochs of evaluation. Columns 2-4 provide position components. Columns 5-7 provide velocity components. Column 8-10, when present, provide acceleration components.
https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/req/spk.html https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/req/naif_ids.html https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/req/frames.html Shampine, L. F. and Gordon, M. K., Computer Solution of Ordinary Differential Equations: The Initial Value Problem, 1975 Robert Werner, SPICE spke01 math, 2022. https://doi.org/10.5270/esa-tyidsbu
# The file vgr2_jup230.bsp provided with the package includes information for the # Jupiter flyby of Voyager 2 testSPK <- readBinSPK(paste0(path.package("asteRisk"), "/vgr2_jup230.bsp")) length(testSPK$segments) # It contains a single segment. testSegment <- testSPK$segments[[1]] # Check the initial and end epochs of the interval covered by the segment testSegment$segmentSummary$initialEpoch testSegment$segmentSummary$finalEpoch # Evaluate at target epochs evaluateSPKSegment(testSegment, c(-649364400, -649364374.68, -647364600, -645364800))# The file vgr2_jup230.bsp provided with the package includes information for the # Jupiter flyby of Voyager 2 testSPK <- readBinSPK(paste0(path.package("asteRisk"), "/vgr2_jup230.bsp")) length(testSPK$segments) # It contains a single segment. testSegment <- testSPK$segments[[1]] # Check the initial and end epochs of the interval covered by the segment testSegment$segmentSummary$initialEpoch testSegment$segmentSummary$finalEpoch # Evaluate at target epochs evaluateSPKSegment(testSegment, c(-649364400, -649364374.68, -647364600, -645364800))
The GCRF (Geocentric Celestial Reference Frame) frame of reference is an Earth-centered inertial coordinate frame, where the origin is placed at the center of mass of Earth and the coordinate frame is fixed with respect to the stars (and therefore not fixed with respect to the Earth surface in its rotation). The X-axis is aligned with the mean equinox of Earth at 12:00 Terrestrial Time on the 1st of January, 2000, and the Z-axis is aligned with the Earth´s rotation axis.
It is almost equivalent to the J2000 frame of reference (also called EME2000), and in some contexts it is also referred to as ICRF frame (although in its strict definition, the origin of coordinates is placed at the barycenter of the Solar System).
In the ITRF (International Terrestrial Reference Frame), the origin is also placed at the center of mass of Earth, but the frame rotates with respect to the stars to remain fixed with respect to the Earth surface as it rotates. The Z-axis extends along the true North as defined by the IERS reference pole, and the X-axis extends towards the intersection between the equator and the Greenwich meridian at any time.
The coordinates and velocities input and calculated with the high-precision orbital propagator (hpop) are in the GCRF frame of reference.
This function requires the asteRiskData package, which can be installed
by running install.packages('asteRiskData', repos='https://rafael-ayala.github.io/drat/')
GCRFtoITRF(position_GCRF, velocity_GCRF, dateTime)GCRFtoITRF(position_GCRF, velocity_GCRF, dateTime)
position_GCRF |
Vector with the X, Y and Z components of the position of an object in GCRF frame, in m. |
velocity_GCRF |
Vector with the X, Y and Z components of the velocity of an object in GCRF frame, in m/s. |
dateTime |
Date-time string with the date and time in UTC corresponding to the provided position and velocity vectors. This specifies the time for which the conversion from GCRF to ITRF coordinates will be performed. It is important to provide an accurate value, since the point over the surface of Earth to which a set of GCRF coordinates refers varies with time due to the motion of Earth. |
A list with two elements representing the position and velocity of the satellite in the ITRF (International Terrestrial Reference Frame) frame of reference. Position values are in m, and velocity values are in m/s. Each of the two elements contains three values, corresponding to the X, Y and Z components of position and velocity in this order.
https://celestrak.org/columns/v02n01/
if(requireNamespace("asteRiskData", quietly = TRUE)) { # The following were the position and velocity of satellite MOLNIYA 1-83 # the 25th of June, 2006 at 00:33:43 UTC in the GCRF frame (in m and m/s). position_GCRF <- c(-14471729.582, -4677558.558, 9369.461) velocity_GCRF <- c(-3251.691, -3276.008, 4009.228) # Let´s convert them into the ITRF frame coordinates_ITRF <- GCRFtoITRF(position_GCRF, velocity_GCRF, "2006-06-27 00:58:29.34") }if(requireNamespace("asteRiskData", quietly = TRUE)) { # The following were the position and velocity of satellite MOLNIYA 1-83 # the 25th of June, 2006 at 00:33:43 UTC in the GCRF frame (in m and m/s). position_GCRF <- c(-14471729.582, -4677558.558, 9369.461) velocity_GCRF <- c(-3251.691, -3276.008, 4009.228) # Let´s convert them into the ITRF frame coordinates_ITRF <- GCRFtoITRF(position_GCRF, velocity_GCRF, "2006-06-27 00:58:29.34") }
The GCRF (Geocentric Celestial Reference Frame) frame of reference is an Earth-centered inertial coordinate frame, where the origin is placed at the center of mass of Earth and the coordinate frame is fixed with respect to the stars (and therefore not fixed with respect to the Earth surface in its rotation). The X-axis is aligned with the mean equinox of Earth at 12:00 Terrestrial Time on the 1st of January, 2000, and the Z-axis is aligned with the Earth´s rotation axis. This function converts position in GCRF to geodetic latitude, longitude and altitude, which can be considered to be a non-inertial, Earth-centered frame of reference.
This function requires the asteRiskData package, which can be installed
by running install.packages('asteRiskData', repos='https://rafael-ayala.github.io/drat/')
GCRFtoLATLON(position_GCRF, dateTime, degreesOutput=TRUE)GCRFtoLATLON(position_GCRF, dateTime, degreesOutput=TRUE)
position_GCRF |
Vector with the X, Y and Z components of the position of an object in TEME frame, in m. |
dateTime |
Date-time string with the date and time in UTC corresponding to the provided position vector. This specifies the time for which the conversion from GCRF to geodetic coordinates will be performed. It is important to provide an accurate value, since the point over the surface of Earth to which a set of GCRF coordinates refers varies with time due to the motion of Earth. |
degreesOutput |
Logical indicating if the output should be in sexagesimal
degrees. If |
A vector with three elements, corresponding to the latitude and longitude in degrees (or radians if specified) and the altitude in m.
https://arc.aiaa.org/doi/10.2514/6.2006-6753
if(requireNamespace("asteRiskData", quietly = TRUE)) { # The following orbital parameters correspond to an object with NORAD catalogue # number 24208 (Italsat 2) the 26th of June, 2006 at 00:58:29.34 UTC. n0 <- 1.007781*((2*pi)/(1440)) # Multiplication by 2pi/1440 to convert to radians/min e0 <- 0.002664 # mean eccentricity at epoch i0 <- 3.8536*pi/180 # mean inclination at epoch in radians M0 <- 48.3*pi/180 # mean anomaly at epoch in radians omega0 <- 311.0977*pi/180 # mean argument of perigee at epoch in radians OMEGA0 <- 80.0121*pi/180 # mean longitude of ascending node at epoch in radians Bstar <- 1e-04 # drag coefficient epochDateTime <- "2006-06-26 00:58:29.34" # Let´s calculate the position and velocity of the satellite 1 day later state_1day_TEME <- sgdp4(n0=n0, e0=e0, i0=i0, M0=M0, omega0=omega0, OMEGA0=OMEGA0, Bstar=Bstar, initialDateTime=epochDateTime, targetTime=1440) # We can now convert the results in TEME frame to GCRF frame, previously # multiplying by 1000 to convert the km output of sgdp4 to m state_1day_GCRF <- TEMEtoGCRF(state_1day_TEME$position*1000, state_1day_TEME$velocity*1000, "2006-06-27 00:58:29.34") # Finally, we convert the results in GCRF frame to geodetic latitude, longitude # and altitude state_1day_geodetic <- GCRFtoLATLON(state_1day_GCRF$position, "2006-06-27 00:58:29.34") }if(requireNamespace("asteRiskData", quietly = TRUE)) { # The following orbital parameters correspond to an object with NORAD catalogue # number 24208 (Italsat 2) the 26th of June, 2006 at 00:58:29.34 UTC. n0 <- 1.007781*((2*pi)/(1440)) # Multiplication by 2pi/1440 to convert to radians/min e0 <- 0.002664 # mean eccentricity at epoch i0 <- 3.8536*pi/180 # mean inclination at epoch in radians M0 <- 48.3*pi/180 # mean anomaly at epoch in radians omega0 <- 311.0977*pi/180 # mean argument of perigee at epoch in radians OMEGA0 <- 80.0121*pi/180 # mean longitude of ascending node at epoch in radians Bstar <- 1e-04 # drag coefficient epochDateTime <- "2006-06-26 00:58:29.34" # Let´s calculate the position and velocity of the satellite 1 day later state_1day_TEME <- sgdp4(n0=n0, e0=e0, i0=i0, M0=M0, omega0=omega0, OMEGA0=OMEGA0, Bstar=Bstar, initialDateTime=epochDateTime, targetTime=1440) # We can now convert the results in TEME frame to GCRF frame, previously # multiplying by 1000 to convert the km output of sgdp4 to m state_1day_GCRF <- TEMEtoGCRF(state_1day_TEME$position*1000, state_1day_TEME$velocity*1000, "2006-06-27 00:58:29.34") # Finally, we convert the results in GCRF frame to geodetic latitude, longitude # and altitude state_1day_geodetic <- GCRFtoLATLON(state_1day_GCRF$position, "2006-06-27 00:58:29.34") }
The asteRiskData package provides the data and coefficients required
for calculation of forces for hpop and other functions such certain
conversions between reference frames. Some of the data tables included in the
package are updated periodically with new data. These include Earth orientation
parameters, space weather data and solar and geomagnetic storms. In order to
perform the calculations dependent on such data for the most recent times, the
latest available data must be retrieved.
This function automatically updates the data tables, enabling such calculations for the most recent times.
getLatestSpaceData(targets="all")getLatestSpaceData(targets="all")
targets |
Character vector specifying the data that should be updated. It should be a vector containing one or more of the following strings: "all" (to update all data), "EOP" (Earth orientation parameters), "SW" (space weather), "SS" (solar storms) or "GS" (geomagnetic storms). By default, all data are updated. |
This function is invoked for its side effect, which is updating the data tables
used internally for calculations requiring asteRiskData package, such as
those performed by hpop.
http://www.celestrak.org/SpaceData/EOP-All.txt https://celestrak.org/SpaceData/SW-All.txt https://sol.spacenvironment.net/jb2008/indices.html
if(interactive()) { if(requireNamespace("asteRiskData", quietly = TRUE)) { # The table of Earth orientation parameters distributed with asteRiskData # comprises data up to the 21st of March, 2021 asteRiskData::earthPositions[nrow(asteRiskData::earthPositions),] # The table can be easily updated to include the most recent available data getLatestSpaceData(targets="all") asteRiskData::earthPositions[nrow(asteRiskData::earthPositions),] } }if(interactive()) { if(requireNamespace("asteRiskData", quietly = TRUE)) { # The table of Earth orientation parameters distributed with asteRiskData # comprises data up to the 21st of March, 2021 asteRiskData::earthPositions[nrow(asteRiskData::earthPositions),] # The table can be easily updated to include the most recent available data getLatestSpaceData(targets="all") asteRiskData::earthPositions[nrow(asteRiskData::earthPositions),] } }
Given the position and velocity of a satellite at a given time (in the ICRF
system of coordinates centered on the Solar System Barycenter, any of the main planets,
Earth's Moon or Pluto), propagates its position by calculating its acceleration
(based on a force model) and solving the resulting second-order ODE through
numerical integration. This allows propagation of orbits with considerably
higher accuracy than other propagators such as SGP4 and SDP4, but at the expense
of a much higher computational cost. The forces and effects currently considered
are gravitational attraction by the Earth (using the GGM05C gravity model, with
spherical harmonics up to degree and order of 360); effects of Earth ocean and
solid tides; gravitational attraction by the Moon (using the GRGM1200B gravity
model with spherical harmonics up to degree and order of 1200), effects of solid
Moon tides (currently using an ellastic Moon model), Sun and planets (considered
as point masses); solar radiation pressure; atmospheric drag, and relativistic
effects. The force field is based on the forces described in Satellite Orbits:
Models, Methods and Applications (Oliver Montenbruck and Eberhard Gill) and
Fundamentals of Astrodynamics and Applications (David Vallado).
The NRLMSISE-00 model is used to calculate atmospheric density for the
calculation of atmospheric drag. The FES2014 model is used to calculate Earth
geopotential model corrections due to ocean tides.
As mentioned before, the central body for the frame of reference can be any of the
Solar System Barycenter (SSB), any of the main planets, Earth's Moon or Pluto.
By default, it is assumed to be Earth, corresponding to GCRF (Geocentric ICRF).
The initial position will be checked against the position of said celestial bodies,
to identify if it falls under the Laplacian gravitational sphere of influence of any of them.
If this is the case, and it differs from the specified central body, the coordinate
system will be changed to be centered on the celestial body whose sphere of influence
includes the object of interest. This avoids instability in propagation.
The high-precision numerical orbital propagator requires the asteRiskData
package, which provides the data and coefficients required for calculation of
the modeled forces. asteRiskData can be installed by running
install.packages('asteRiskData', repos='https://rafael-ayala.github.io/drat/')
hpop(position, velocity, dateTime, times, satelliteMass, dragArea, radiationArea, dragCoefficient, radiationCoefficient, earthSphericalHarmonicsDegree=130, solidEarthTides=TRUE, oceanTides=TRUE, moonSphericalHarmonicsDegree=150, solidMoonTides=TRUE, centralBody="Earth", autoCentralBodyChange=TRUE, ...)hpop(position, velocity, dateTime, times, satelliteMass, dragArea, radiationArea, dragCoefficient, radiationCoefficient, earthSphericalHarmonicsDegree=130, solidEarthTides=TRUE, oceanTides=TRUE, moonSphericalHarmonicsDegree=150, solidMoonTides=TRUE, centralBody="Earth", autoCentralBodyChange=TRUE, ...)
position |
Initial position of the satellite in the GCRF system of coordinates. Should be provided as a numeric vector with 3 components that indicate the X, Y and Z components of the position in meters. |
velocity |
Initial velocity of the satellite in the GCRF system of coordinates. Should be provided as a numeric vector with 3 components that indicate the X, Y and Z components of the position in meters/second. |
dateTime |
Date time string in the YYYY-MM-DD HH:MM:SS format indicating the time corresponding to the initial position and velocity, in UTC time. |
times |
Vector with the times at which the position and velocity of the satellite should be calculated, in seconds since the initial time. |
satelliteMass |
Mass of the satellite in kilograms. |
dragArea |
Effective area of the satellite for atmospheric drag in squared meters. If the way that a satellite will orient with respect to its velocity is not known, a mean cross-sectional area should be calculated assuming that the orientation of the satellite with respect to its velocity will vary uniformly. A decent estimate can be obtained with a flat-plate model, where the satellite is considered to be parallelepiped-shaped. The mean effective area can then be calculated as CSA = (S1 + S2 + S3 (+S4))/2, where S1, S2 and S3 are the areas of the three perpendicular surfaces of the model and S4 is an optional term to account for the area of solar panels (potential masking between the solar panels and the main surfaces is not considered; this might be partially accounted for by introducing a factor to reduce the calculated effective area). |
radiationArea |
Effective area of the satellite subject to the effect of radiation pressure in squared meters. |
dragCoefficient |
Drag coefficient (Cd) used for the calculation of atmospheric drag. For low Earth-orbiting satellites, a value of 2.2 is frequently employed if a better approximation is not available. |
radiationCoefficient |
Coefficient for the force resulting from radiation pressure. This parameter is usually referred to as reflectivity coefficient (Cr) and the value varies for different satellites and orbits. If unknown, a value of 1.2 is usually a decent approximation. |
earthSphericalHarmonicsDegree |
Maximum degree and order that should be considered when calculating the Earth geopotential model. The model will be complete up to the specified degree/order, i.e., all zonal, sectorial and tesseral spherical harmonics will be calculated. The maximum possible value is 360, since that is the highest degree and order of the Stokes' coefficients provided in the GGM05C model. Note that spherical harmonics for Earth gravity field will only be used if Earth is the central body for propagation; otherwise, only a point-mass attraction will be calculated. |
solidEarthTides |
Logical indicating if corrections of the Cnm and Snm Stokes' coefficients for the geopotential model due to solid Earth tides should be performed, following IERS 2010 procedures and considering anelasticity of the Earth. |
oceanTides |
Logical indicating if corrections of the Cnm and Snm Stokes' coefficients for the geopotential model due to ocean tides should be performed, using the FES2014 oceanic tides model. |
moonSphericalHarmonicsDegree |
Maximum degree and order that should be considered when calculating the Moon gravity model. The model will be complete up to the specified degree/order, i.e., all zonal, sectorial and tesseral spherical harmonics will be calculated. The maximum possible value is 1200, since that is the highest degree and order of the Stokes' coefficients provided in the GRGM1200B model. Note that spherical harmonics for Moon gravity field will only be used if Moon is the central body for propagation; otherwise, only a point-mass attraction will be calculated. |
solidMoonTides |
Logical indicating if corrections of the Cnm and Snm Stokes' coefficients for the lunar gravity model due to solid Moon tides should be performed, following the procedure described by William and Boggs, 2015 using an elastic Moon model. Corrections are applied to the C20, C21, C22, S21 and S22 coefficients. |
centralBody |
Character string indicating the celestial body on which the supplied initial position (in ICRF) are centered. Should be one of "SSB" (meaning Solar System Barycenter), "Mercury", "Venus", "Earth", "Moon", "Mars", "Jupiter", "Saturn", "Uranus", "Neptune" or "Pluto". The initial position will be checked against the position of said celestial bodies, to identify if it falls under the Laplacian gravitational sphere of influence of any of them. If this is the case, and it differs from the specified central body, the coordinate system will be changed to be centered on the celestial body whose sphere of influence includes the object of interest. |
autoCentralBodyChange |
Logical indicating if the celestial object used as
the center of coordinates should be automatically updated during propagation
based on the radii of the spheres of influence of the main planets, the Moon
and Pluto. By default, |
... |
Additional parameters to be passed to ode to control how numerical integration is performed. By default, the RADAU5 solver is used. |
A data frame with the results of the numerical integration at the requested times.
Each row contains the results for one of the requested times. The data frame
contains 11 columns: time (indicating the time for the corresponding row
in seconds since the initial time), positionX, positionY, positionZ (indicating
the X, Y and Z components of the position for that time in meters), velocityX,
velocityY and velocityZ (indicating the X, Y and Z components of the velocity
for that time in meters/second), accelerationX, accelerationY, accelerationZ
(indicating the X, Y and Z components of the acceleration for that time in
meters/second^2) and centralBody, indicating the central body of the frame of
reference for the results for the corresponding time.
Positions and velocities are returned in the ICRF frame of reference, centered
in the celestial body specified in column centralBody. If autoCentralBodyChange=TRUE,
the celestial body whose sphere of influence includes the object of interest
will be automatically used as the central body. Additionally, if transitions in
or out of the spheres of influence of the main celestial bodies are detected
during propagation of the trajectory, the central body will be automatically
modified accordingly. If autoCentralBodyChange=FALSE, such automatic
changes of the central body will not be performed, and instead the user-specified
central body will be used at all times. Note, however, that it is not recommended
to perform propagation in a frame center at an object different than the celestial
body whose sphere of influence includes the target of propagation, since this
can lead to a substantial loss of accuracy. For details, see M. Vautier, 2008.
Note that, if none of the spheres of influence of the
planets, Moon or Pluto included the object of interest, the center of the ICRF
frame will be placed at the Solar System Barycenter.
Satellite Orbits: Models, Methods and Applications. Oliver Montenbruck and Eberhard Gill. Fundamentals of Astrodynamics and Applications. David Vallado. https://www.mathworks.com/matlabcentral/fileexchange/55167-high-precision-orbit-propagator https://ccmc.gsfc.nasa.gov/modelweb/models/nrlmsise00.php https://digitalcommons.usu.edu/cgi/viewcontent.cgi?article=1144&context=smallsat https://iopscience.iop.org/article/10.1088/1742-6596/911/1/012009/pdf https://www.sciencedirect.com/science/article/pii/S1110016821000016 https://etd.auburn.edu/bitstream/handle/10415/1133/Vautier_Mana_34.pdf?sequence=1 https://agupubs.onlinelibrary.wiley.com/doi/full/10.1002/2014JE004755
if(requireNamespace("asteRiskData", quietly = TRUE)) { # The following are the position and velocity in the GCRF frame of satellite # MOLNIYA 1-83 the 25th of June, 2006 at 00:33:43 UTC. initialPosition <-c(-14568679.5026116, -4366250.78287623, 9417.9289105405) initialVelocity <- c(-3321.17428902497, -3205.49400830455, 4009.26862308806) initialTime <- "2006-06-25 00:33:43" # Molniya satellites have a mass of approximately 1600 kg and a cross-section of # 15 m2. Additionally, let´s use 2.2 and 1.2 as approximately values of the # drag and reflectivity coefficients, respectively. molniyaMass <- 1600 molniyaCrossSection <- 15 molniyaCr <- 1.2 molniyaCd <- 2.2 # Let´s calculate the position and velocity of the satellite for each minute of # the following 10 minutes. targetTimes <- seq(0, 600, by=60) hpop_results <- hpop(initialPosition, initialVelocity, initialTime, targetTimes, molniyaMass, molniyaCrossSection, molniyaCrossSection, molniyaCr, molniyaCd) }if(requireNamespace("asteRiskData", quietly = TRUE)) { # The following are the position and velocity in the GCRF frame of satellite # MOLNIYA 1-83 the 25th of June, 2006 at 00:33:43 UTC. initialPosition <-c(-14568679.5026116, -4366250.78287623, 9417.9289105405) initialVelocity <- c(-3321.17428902497, -3205.49400830455, 4009.26862308806) initialTime <- "2006-06-25 00:33:43" # Molniya satellites have a mass of approximately 1600 kg and a cross-section of # 15 m2. Additionally, let´s use 2.2 and 1.2 as approximately values of the # drag and reflectivity coefficients, respectively. molniyaMass <- 1600 molniyaCrossSection <- 15 molniyaCr <- 1.2 molniyaCd <- 2.2 # Let´s calculate the position and velocity of the satellite for each minute of # the following 10 minutes. targetTimes <- seq(0, 600, by=60) hpop_results <- hpop(initialPosition, initialVelocity, initialTime, targetTimes, molniyaMass, molniyaCrossSection, molniyaCrossSection, molniyaCr, molniyaCd) }
The ITRF (International Terrestrial Reference Frame) is an ECEF (Earth Centered, Earth Fixed) frame of reference, i.e., a non-inertial frame of reference where the origin is placed at the center of mass of Earth, and the frame rotates with respect to the stars to remain fixed with respect to the Earth surface as it rotates. The Z-axis extends along the true North as defined by the IERS reference pole, and the X-axis extends towards the intersection between the equator and the Greenwich meridian at any time.
The GCRF (Geocentric Celestial Reference Frame) frame of reference is an Earth-centered inertial coordinate frame, where the origin is also placed at the center of mass of Earth and the coordinate frame is fixed with respect to the stars (and therefore not fixed with respect to the Earth surface in its rotation). The X-axis is aligned with the mean equinox of Earth at 12:00 Terrestrial Time on the 1st of January, 2000, and the Z-axis is aligned with the Earth´s rotation axis.
This function requires the asteRiskData package, which can be installed
by running install.packages('asteRiskData', repos='https://rafael-ayala.github.io/drat/')
ITRFtoGCRF(position_ITRF, velocity_ITRF, dateTime)ITRFtoGCRF(position_ITRF, velocity_ITRF, dateTime)
position_ITRF |
Vector with the X, Y and Z components of the position of an object in ITRF frame, in m. |
velocity_ITRF |
Vector with the X, Y and Z components of the velocity of an object in ITRF frame, in m/s. |
dateTime |
Date-time string with the date and time in UTC corresponding to the provided position and velocity vectors. This specifies the time for which the conversion from ITRF to GCRF coordinates will be performed. It is important to provide an accurate value, since the point over the surface of Earth to which a set of GCRF coordinates refers varies with time due to the motion of Earth. |
A list with two elements representing the position and velocity of the satellite in the GCRF (Earth-centered non-intertial) frame of reference. Position values are in m, and velocity values are in m/s. Each of the two elements contains three values, corresponding to the X, Y and Z components of position and velocity in this order.
https://celestrak.org/columns/v02n01/
if(requireNamespace("asteRiskData", quietly = TRUE)) { # The following were the position and velocity of satellite MOLNIYA 1-83 # the 25th of June, 2006 at 00:33:43 UTC in the ECEF frame (in m and m/s). position_ITRF <- c(1.734019e+06, -1.510972e+07, 39.08228) velocity_ITRF <- c(1468.832, -3962.464, 4007.039) # Let´s convert them into the GCRF frame coordinates_GCRF <- ITRFtoGCRF(position_ITRF, velocity_ITRF, "2006-06-25 00:33:43") }if(requireNamespace("asteRiskData", quietly = TRUE)) { # The following were the position and velocity of satellite MOLNIYA 1-83 # the 25th of June, 2006 at 00:33:43 UTC in the ECEF frame (in m and m/s). position_ITRF <- c(1.734019e+06, -1.510972e+07, 39.08228) velocity_ITRF <- c(1468.832, -3962.464, 4007.039) # Let´s convert them into the GCRF frame coordinates_GCRF <- ITRFtoGCRF(position_ITRF, velocity_ITRF, "2006-06-25 00:33:43") }
The ITRF (International Terrestrial Reference Frame) is an ECEF (Earth Centered, Earth Fixed) frame of reference, i.e., a non-inertial frame of reference where the origin is placed at the center of mass of Earth, and the frame rotates with respect to the stars to remain fixed with respect to the Earth surface as it rotates. The Z-axis extends along the true North as defined by the IERS reference pole, and the X-axis extends towards the intersection between the equator and the Greenwich meridian at any time. This function converts Cartesian coordinates in the ECEF frame to geodetic latitude, longitude and altitude.
ITRFtoLATLON(position_ITRF, degreesOutput=TRUE)ITRFtoLATLON(position_ITRF, degreesOutput=TRUE)
position_ITRF |
Vector with the X, Y and Z components of the position of an object in ITRF frame, in m. |
degreesOutput |
Logical indicating if the output should be in sexagesimal
degrees. If |
A vector with three elements, corresponding to the latitude and longitude in degrees (or radians if specified) and the altitude in m.
https://arc.aiaa.org/doi/10.2514/6.2006-6753
coordinates_ITRF <- c(5062040.1, -530657.4, 3863936.5) # Let's calculate the geodetic latitude, longitude and altitude geodetic <- ITRFtoLATLON <- (coordinates_ITRF) # A longer application example follows, useful to represent the groundtrack of a # satellite after propagaton with SGP4/SDP4 if(requireNamespace("asteRiskData", quietly = TRUE)) { # The following orbital parameters correspond to an object with NORAD catalogue # number 24208 (Italsat 2) the 26th of June, 2006 at 00:58:29.34 UTC. n0 <- 1.007781*((2*pi)/(1440)) # Multiplication by 2pi/1440 to convert to radians/min e0 <- 0.002664 # mean eccentricity at epoch i0 <- 3.8536*pi/180 # mean inclination at epoch in radians M0 <- 48.3*pi/180 # mean anomaly at epoch in radians omega0 <- 311.0977*pi/180 # mean argument of perigee at epoch in radians OMEGA0 <- 80.0121*pi/180 # mean longitude of ascending node at epoch in radians Bstar <- 1e-04 # drag coefficient epochDateTime <- "2006-06-26 00:58:29.34" # Let´s calculate the position and velocity of the satellite 1 day later state_1day_TEME <- sgdp4(n0=n0, e0=e0, i0=i0, M0=M0, omega0=omega0, OMEGA0=OMEGA0, Bstar=Bstar, initialDateTime=epochDateTime, targetTime=1440) # We can now convert the results in TEME frame to ITRF frame, previously # multiplying by 1000 to convert the km output of sgdp4 to m state_1day_ITRF <- TEMEtoITRF(state_1day_TEME$position, state_1day_TEME$velocity, "2006-06-27 00:58:29.34") # Finally, we can convert the ECEF coordinates to geodetic latitude, longitude # and altitude state_1day_geodetic <- ITRFtoLATLON(state_1day_ITRF$position) }coordinates_ITRF <- c(5062040.1, -530657.4, 3863936.5) # Let's calculate the geodetic latitude, longitude and altitude geodetic <- ITRFtoLATLON <- (coordinates_ITRF) # A longer application example follows, useful to represent the groundtrack of a # satellite after propagaton with SGP4/SDP4 if(requireNamespace("asteRiskData", quietly = TRUE)) { # The following orbital parameters correspond to an object with NORAD catalogue # number 24208 (Italsat 2) the 26th of June, 2006 at 00:58:29.34 UTC. n0 <- 1.007781*((2*pi)/(1440)) # Multiplication by 2pi/1440 to convert to radians/min e0 <- 0.002664 # mean eccentricity at epoch i0 <- 3.8536*pi/180 # mean inclination at epoch in radians M0 <- 48.3*pi/180 # mean anomaly at epoch in radians omega0 <- 311.0977*pi/180 # mean argument of perigee at epoch in radians OMEGA0 <- 80.0121*pi/180 # mean longitude of ascending node at epoch in radians Bstar <- 1e-04 # drag coefficient epochDateTime <- "2006-06-26 00:58:29.34" # Let´s calculate the position and velocity of the satellite 1 day later state_1day_TEME <- sgdp4(n0=n0, e0=e0, i0=i0, M0=M0, omega0=omega0, OMEGA0=OMEGA0, Bstar=Bstar, initialDateTime=epochDateTime, targetTime=1440) # We can now convert the results in TEME frame to ITRF frame, previously # multiplying by 1000 to convert the km output of sgdp4 to m state_1day_ITRF <- TEMEtoITRF(state_1day_TEME$position, state_1day_TEME$velocity, "2006-06-27 00:58:29.34") # Finally, we can convert the ECEF coordinates to geodetic latitude, longitude # and altitude state_1day_geodetic <- ITRFtoLATLON(state_1day_ITRF$position) }
NASA's Jet Propulsion Laboratory (JPL) provides mathematical models of the Solar
System known as Development Ephemerides (DE). The models are given as sets of
Chebyshev coefficients, which cam be used to calculate the position (and its
derivatives) of the Sun, the eight major planets, Pluto and the Moon.
This function employes JPL DE440 to calculate the position (and optionally
velocities also) of the mentioned celestial objects, in ICRF frame.
JPL DE440 covers the period from 1550 to 2650 AC. In addition to the position of
celestial objects, lunar libration angles are also calculated. Internally,
calculations are performed by employing Clenshaw's algorithm together with the
Chebyshev coefficients provided by JPL DE440.
The target time should be specified as a Modified Julian Date (MJD). MJD in different
time systems can be used. Currently, UTC, UT1, TT and TDB are supported.
Additionally, a central body with respect to which positions and velocities are
calculated should be specified. By default, the Solar System Barycenter (SSB) is
used, but additionally Mercury, Venus, Earth, Moon, Mars, Jupiter, Saturn, Uranus,
Neptune or Pluto can be selected.
Note that this function requires the additional package asteRiskData, which
provides the Chebyshev coefficients, and can be installed by running
install.packages("asteRiskData", repos="https://rafael-ayala.github.io/drat/")
JPLephemerides(MJD, timeSystem="UTC", centralBody="SSB", derivatives="acceleration")JPLephemerides(MJD, timeSystem="UTC", centralBody="SSB", derivatives="acceleration")
MJD |
Modified Julian Date of the time for which celestial object ephemerides should be calculated. MJD are fractional number of days since midnight of the 17th of November, 1858. The MJD of a date-time string can be obtained with function dateTimeToMJD. |
timeSystem |
Time system into which the MJD is provided. Should be one from "UTC" (Coordinated Universal Time; default), "UT1" (Universal Time), "TT" (Terrestrial Time) and "TDB" (Barycentric Dynamical Time). |
centralBody |
String indicating the celestial object that will be taken as the center of coordinates to which positions and velocities are referred. Must be one of "SSB" (Solar System Barycenter), "Mercury", "Venus", "Earth", "Moon", "Mars", "Jupiter", "Saturn", "Uranus", "Neptune" or "Pluto". |
derivatives |
String indicating what derivatives of positions should be calculated. Must be one of "none", "velocity" or "acceleration". If "none", only position is calculated. If "velocity", velocities are calculated, as well as first order derivatives of Moon libration angles. If "acceleration", both velocities and accelerations (as well as second order derivatives of Moon libration angles) are calculated. |
A list of vectors providing the positions (in meters), velocities (in m/s; only if requested), accelerations (in m/s^2; only if requested), Moon libration angles (in radians), first derivatives of Moon libration angles (in radians/s; only if velocities were requested) and second derivatives of Moon libration angles (in radians/s^2; only if accelerations were requested) of celestial objects with respect to the specified central body. For position, velocity and acceleration vectors, X, Y and Z components are given in this order. For Moon libration angles and their derivatives, they are given in the following order: phi, theta and psi.
https://gssc.esa.int/navipedia/index.php/Julian_Date https://gssc.esa.int/navipedia/index.php/Transformations_between_Time_Systems
if(requireNamespace("asteRiskData", quietly = TRUE)) { # Let's calculate the MJD of the 12th of June, 2000 at 10:00:00 UTC time, in UTC MJD_UTC <- dateTimeToMJD("2000-06-12 10:00:00", timeSystem = "UTC") # Let's now calculate the JPL ephemerides using Earth as the central body: ephemerides <- JPLephemerides(MJD_UTC, timeSystem = "UTC", centralBody="Earth") # We can now calculate, for example, the exact distance between the barycenters # of Earth and Moon sqrt(sum(ephemerides$positionMoon^2)) }if(requireNamespace("asteRiskData", quietly = TRUE)) { # Let's calculate the MJD of the 12th of June, 2000 at 10:00:00 UTC time, in UTC MJD_UTC <- dateTimeToMJD("2000-06-12 10:00:00", timeSystem = "UTC") # Let's now calculate the JPL ephemerides using Earth as the central body: ephemerides <- JPLephemerides(MJD_UTC, timeSystem = "UTC", centralBody="Earth") # We can now calculate, for example, the exact distance between the barycenters # of Earth and Moon sqrt(sum(ephemerides$positionMoon^2)) }
Keplerian orbital elements are a set of six parameters used to described the orbits of celestial objects, including satellites. While satellites do not follow a perfectly Keplerian orbit, their state at any point can be defined by the orbital parameters that they would have if they were located at the same position with the same velocity following a perfectly Keplerian orbit (i.e., if perturbations were absent). These are called osculating orbital elements.
A complete set of six Keplerian elements defines unequivocally the position and velocity of the satellite in a given frame of reference, and therefore can be used to calculate its cartesian coordinates. This function calculates the coordinates of a satellite in an ECI (Earth-centered inertial) frame of reference from a set of Keplerian orbital elements. The exact ECI frame of the resulting coordinates is the same used to define the supplied orbital elements.
KOEtoECI(a, e, i, M, omega, OMEGA, keplerAccuracy=10e-8, maxKeplerIterations=100)KOEtoECI(a, e, i, M, omega, OMEGA, keplerAccuracy=10e-8, maxKeplerIterations=100)
a |
Semi-major axis of orbital ellipse in meters. |
e |
Numerical eccentricity of the orbit. Eccentricity measures how much the orbit deviates from being circular. |
i |
Inclination of the orbital plane in radians. Inclination is the angle between the orbital plane and the equator. |
M |
Mean anomaly of the orbit in radians. Mean anomaly indicates where the satellite is along its orbital path, and is defined as the angle between the direction of the perigee and the hypothetical point where the object would be if it was moving in a circular orbit with the same period as its true orbit after the same amount of time since it last crossed the perigee had ellapsed. |
omega |
Argument of perigee in radians. This is the angle between the direction of the ascending node and the direction of the perigee (the point of the orbit at which the object is closest to the Earth). |
OMEGA |
Right ascension of the ascending node in radians. This is the angle between the direction of the ascending node (the point where the satellite crosses the equatorial plane moving north) and the direction of the First Point of Aries (which indicates the location of the vernal equinox). |
keplerAccuracy |
Accuracy to consider Kepler's equation solved when calculating eccentric anomaly from mean anomaly. If two consecutive solutions differ by a value lower than this accuracy, integration is considered to have converged. |
maxKeplerIterations |
Maximum number of iterations after which fixed-point integration of Kepler's equation will stop, even if convergence according to the accuracy criterion has not been reached. |
A list with two elements representing the position and velocity of the satellite in the same ECI (Earth Centered, Earth Fixed) frame of reference into which the provided orbital elements were defined. Position values are in m, and velocity values are in m/s. Each of the two elements contains three values, corresponding to the X, Y and Z components of position and velocity in this order.
https://www.gsc-europa.eu/system-service-status/orbital-and-technical-parameters https://celestrak.org/columns/v02n01/ https://downloads.rene-schwarz.com/download/M001-Keplerian_Orbit_Elements_to_Cartesian_State_Vectors.pdf
# Let's calculate the ECI coordinates from the orbital elements provided by a # TLE. It should be noted that this is often not recommended, since the orbital # elements supplied in a TLE are not osculating orbital elements, but instead # mean orbital elements set to fit a range of actual observations. The # recommended procedures are to use TLE only in conjunction with the SGP4/SDP4 # models, and viceversa. # The following orbital parameters correspond to an object with NORAD catalogue # number 24208 (Italsat 2) the 26th of June, 2006 at 00:58:29.34 UTC. n0 <- 1.007781*((2*pi)/(86400)) # Multiplication by 2pi/86400 to convert to radians/s e0 <- 0.002664 # mean eccentricity at epoch i0 <- 3.8536*pi/180 # mean inclination at epoch in radians M0 <- 48.3*pi/180 # mean anomaly at epoch in radians omega0 <- 311.0977*pi/180 # mean argument of perigee at epoch in radians OMEGA0 <- 80.0121*pi/180 # mean longitude of ascending node at epoch in radians # The semi-major axis can be calculated from the mean motion in radians/s # as follows: (mu is the standard gravitational parameter of Earth) mu <- 3.986004418e14 # in units of m3 s-2 a0 <- (mu^(1/3))/(n0^(2/3)) # The ECI coordinates can then be calculated. In this case, they will be in TEME # frame, since the original orbital elements are derived from a TLE coordinates_ECI <- KOEtoECI(a0, e0, i0, M0, omega0, OMEGA0)# Let's calculate the ECI coordinates from the orbital elements provided by a # TLE. It should be noted that this is often not recommended, since the orbital # elements supplied in a TLE are not osculating orbital elements, but instead # mean orbital elements set to fit a range of actual observations. The # recommended procedures are to use TLE only in conjunction with the SGP4/SDP4 # models, and viceversa. # The following orbital parameters correspond to an object with NORAD catalogue # number 24208 (Italsat 2) the 26th of June, 2006 at 00:58:29.34 UTC. n0 <- 1.007781*((2*pi)/(86400)) # Multiplication by 2pi/86400 to convert to radians/s e0 <- 0.002664 # mean eccentricity at epoch i0 <- 3.8536*pi/180 # mean inclination at epoch in radians M0 <- 48.3*pi/180 # mean anomaly at epoch in radians omega0 <- 311.0977*pi/180 # mean argument of perigee at epoch in radians OMEGA0 <- 80.0121*pi/180 # mean longitude of ascending node at epoch in radians # The semi-major axis can be calculated from the mean motion in radians/s # as follows: (mu is the standard gravitational parameter of Earth) mu <- 3.986004418e14 # in units of m3 s-2 a0 <- (mu^(1/3))/(n0^(2/3)) # The ECI coordinates can then be calculated. In this case, they will be in TEME # frame, since the original orbital elements are derived from a TLE coordinates_ECI <- KOEtoECI(a0, e0, i0, M0, omega0, OMEGA0)
Given 2 position vectors and a time difference between them, calculates the velocity that the object should have at the initial and final states to follow a transfer orbit in the target time (Lambert's problem). The transfer orbit can be elliptical, parabolic or hyperbolic, all of which are dealt with. In the case of elliptical orbits, multi-revolution orbits are possible. Additionally, low-path and high-path orbits are also possible. All possible solutions are returned for the case of elliptical transfers. The user must specify if a retrograde transfer is desired. By default, prograde transfers are calculated. Currently, the formulation by Dario Izzo is applied to solve Lambert's problem (https://link.springer.com/article/10.1007/s10569-014-9587-y).
lambert(initialPosition, finalPosition, initialTime, finalTime, retrogradeTransfer=FALSE, centralBody="Earth", maxIterations=2000, atol=0.00001, rtol=0.000001)lambert(initialPosition, finalPosition, initialTime, finalTime, retrogradeTransfer=FALSE, centralBody="Earth", maxIterations=2000, atol=0.00001, rtol=0.000001)
initialPosition |
Vector with the 3 components of the initial position in Cartesian coordinates, in meters. |
finalPosition |
Vector with the 3 components of the initial position in Cartesian coordinates, in meters. |
initialTime |
Either date-time string in UTC indicating the time
corresponding to the initial state vector of the satellite, or numeric value
indicating the starting time in seconds since an arbitrary reference instant.
If provided as a date-time string in UTC, |
finalTime |
Either date-time string in UTC indicating the time
corresponding to the final state vector of the satellite, or numeric value
indicating the final time in seconds.
If |
retrogradeTransfer |
Logical indicating if retrograde transfer orbits
should be calculated, i.e., transfer orbits with an inclination higher than
180º. By default, |
centralBody |
String indicating the central body around which the satellite
is orbiting. Can be one of |
maxIterations |
Maximum number of iterations to perform at the different root-finding steps when solving Lambert's problem. |
atol |
Absolute tolerance value used for the convergence criterion at the different root-finding steps when solving Lambert's problem. |
rtol |
Relative tolerance value used for the convergence criterion at the different root-finding steps when solving Lambert's problem. |
A list with a number of elements equal to the number of possible transfer orbits found. For hyperbolic and parabolic transfer orbits, there will always be a single possible transfer orbit. For elliptic transfer orbits, multi-revolution orbits may exist if the transfer time is large enough. If they are possible, they will be calculated together with the basic, single-revolution orbit. Each element of the top-level list is in itself a list with the following elements:
numberRevs |
Number of complete revolutions that will be performed in this transfer orbit. Always equals 0 for parabolic and hyperbolic transfer orbits, as well as for the non-multirevolution orbit of elliptic transfers. |
path |
String indicating if the transfer orbit is of the high-path type (i.e., has its second focus located beyond the vector connecting the initial and final positions) or of the low-path type (second focus located between this vector and the first focus). |
orbitType |
String indicating the type of transfer orbit (elliptic, parabolic or hyperbolic). |
v1 |
Velocity vector in m/s that the object should have at the start of the transfer orbit. |
v2 |
Velocity vector in m/s that the object should have at the end of the transfer orbit. |
https://link.springer.com/article/10.1007/s10569-014-9587-y
# Consider the following initial and final positions: initialPosition <- c(15945.34, 0, 0) * 1000 finalPosition <- c(12214.83899, 10249.46731, 0) * 1000 # Given a time difference of 76 minutes between the 2 states, calculate the # velocity that the spacecraft should have at the beginning and end of the # transfer orbit lambertSolution <- lambert(initialPosition, finalPosition, 0, 76*60) length(lambertSolution) lambertSolution[[1]]$orbitType lambertSolution[[1]]$v1 lambertSolution[[1]]$v2 # A single transfer orbit is possible (elliptic, single-revolution orbit)# Consider the following initial and final positions: initialPosition <- c(15945.34, 0, 0) * 1000 finalPosition <- c(12214.83899, 10249.46731, 0) * 1000 # Given a time difference of 76 minutes between the 2 states, calculate the # velocity that the spacecraft should have at the beginning and end of the # transfer orbit lambertSolution <- lambert(initialPosition, finalPosition, 0, 76*60) length(lambertSolution) lambertSolution[[1]]$orbitType lambertSolution[[1]]$v1 lambertSolution[[1]]$v2 # A single transfer orbit is possible (elliptic, single-revolution orbit)
The GCRF (Geocentric Celestial Reference Frame) frame of reference is an Earth-centered inertial coordinate frame, where the origin is also placed at the center of mass of Earth and the coordinate frame is fixed with respect to the stars (and therefore not fixed with respect to the Earth surface in its rotation). The X-axis is aligned with the mean equinox of Earth at 12:00 Terrestrial Time on the 1st of January, 2000, and the Z-axis is aligned with the Earth´s rotation axis. This function converts geodetic latitude, longitude and altitude to Cartesian coordinates in the GCRF frame. The WGS84 Earth ellipsoid model is used.
LATLONtoGCRF(position_LATLON, dateTime, degreesInput=TRUE)LATLONtoGCRF(position_LATLON, dateTime, degreesInput=TRUE)
position_LATLON |
Vector with the latitude, longitude and altitude of the object. Latitude and longitude can be provided in sexagesimal degrees or in radians (by default, sexagesimal degrees are asumed). Altitude must be provided in meters. |
dateTime |
Date-time string with the date and time in UTC corresponding to the provided geodetic coordinates. |
degreesInput |
Logical indicating if the input latitude and longitude are
in sexagesimal degrees. If |
A vector with three elements, corresponding to the X, Y and Z components of position in meters in the ECEF frame, in this order.
https://apps.dtic.mil/sti/pdfs/ADA280358.pdf
if(requireNamespace("asteRiskData", quietly = TRUE)) { latitude <- 37.3891 longitude <- -5.9845 altitude <- 20000 # Let´s calculate the corresponding coordinates in GCRF frame, assuming that # the provided coordinates were valid the 20th of October, 2020 at 15:00:00 UTC coordinatesGCRF <- LATLONtoGCRF(c(latitude, longitude, altitude), dateTime="2020-10-20 15:00:00") }if(requireNamespace("asteRiskData", quietly = TRUE)) { latitude <- 37.3891 longitude <- -5.9845 altitude <- 20000 # Let´s calculate the corresponding coordinates in GCRF frame, assuming that # the provided coordinates were valid the 20th of October, 2020 at 15:00:00 UTC coordinatesGCRF <- LATLONtoGCRF(c(latitude, longitude, altitude), dateTime="2020-10-20 15:00:00") }
The ITRF (International Terrestrial Reference Frame) is an ECEF (Earth Centered, Earth Fixed) frame of reference, i.e., a non-inertial frame of reference where the origin is placed at the center of mass of Earth, and the frame rotates with respect to the stars to remain fixed with respect to the Earth surface as it rotates. The Z-axis extends along the true North as defined by the IERS reference pole, and the X-axis extends towards the intersection between the equator and the Greenwich meridian at any time. This function converts geodetic latitude, longitude and altitude to Cartesian coordinates in the ITRF frame. The WGS84 Earth ellipsoid model is used.
LATLONtoITRF(position_LATLON, degreesInput=TRUE)LATLONtoITRF(position_LATLON, degreesInput=TRUE)
position_LATLON |
Vector with the latitude, longitude and altitude of the object. Latitude and longitude can be provided in sexagesimal degrees or in radians (by default, sexagesimal degrees are asumed). Altitude must be provided in meters. |
degreesInput |
Logical indicating if the input latitude and longitude are
in sexagesimal degrees. If |
A vector with three elements, corresponding to the X, Y and Z components of position in meters in the ITRF frame, in this order.
https://apps.dtic.mil/sti/pdfs/ADA280358.pdf
latitude <- 37.3891 longitude <- -5.9845 altitude <- 20000 # Let´s calculate the corresponding coordinates in ECEF frame coordinates_ITRF <- LATLONtoITRF(c(latitude, longitude, altitude))latitude <- 37.3891 longitude <- -5.9845 altitude <- 20000 # Let´s calculate the corresponding coordinates in ECEF frame coordinates_ITRF <- LATLONtoITRF(c(latitude, longitude, altitude))
TLE (Two-/Three- Line Element) is the standard format for representing orbital state vectors. This function parses a character vector where each element represents a line of the TLE. The supplied character vector can have either 2 (for Two Line Elements) or 3 (for Three Line Elements) elements. The two lines of a Two Line Element contain all the information. The additional line in a Three Line Element is optional, and contains just the satellite name. For a detailed description of the TLE format, see https://celestrak.org/columns/v04n03/#FAQ01.
parseTLElines(lines)parseTLElines(lines)
lines |
Character vector where each element is a string corresponding to a line of the TLE. The character vector must have either 2 or 3 elements. |
A list with the following elements that define the orbital state vector of the satellite:
NORADcatalogNumber |
NORAD Catalog Number, also known as Satellite Catalog Number, assigned by United States Space Command to each artificial object orbiting Earth |
classificationLevel |
Classification level of the information for the orbiting object. Can be unclassified, classified, secret or unknown |
internationalDesignator |
International Designator, also known as COSPAR ID, of the object. It consists of the launch year, separated by a hyphen from a three-digit number indicating the launch number for that year and a set of one to three letters indicating the piece for a launch with multiple pieces. |
launchYear |
The launch year of the object |
launchNumber |
The launch number of the object during its launch year |
launchPiece |
The piece for the launch of the object, if it was a launch with multiple pieces |
dateTime |
Date time string to which the orbital state vector corresponds |
elementNumber |
Element number for the object. In principle, every time a new TLE is generated for an object, the element number is incremented, and therefore element numbers could be used to assess if all the TLEs for a certain object are available. However, in practice it is observed that this is not always the case, with some numbers skipped and some numbers repeated. |
inclination |
Mean orbital inclination of the satellite in degrees. This is the angle between the orbital plane of the satellite and the equatorial plane |
ascension |
Mean longitude of the ascending node of the satellite at epoch, also known as right ascension of the ascending node, in degrees. This is the angle between the direction of the ascending node (the point where the satellite crosses the equatorial plane moving north) and the direction of the First Point of Aries (which indicates the location of the vernal equinox) |
eccentricity |
Mean eccentricity of the orbit of the object. Eccentricity is a measurement of how much the orbit deviates from a circular shape, with 0 indicating a perfectly circular orbit and 1 indicating an extreme case of parabolic trajectory |
perigeeArgument |
Mean argument of the perigee of the object in degrees. This is the angle between the direction of the ascending node and the direction of the perigee (the point of the orbit at which the object is closest to the Earth) |
meanAnomaly |
Mean anomaly of the orbit of the object in degrees. This indicates where the satellite is along its orbital path. It is provided as the angle between the direction of the perigee and the hypothetical point where the object would be if it was moving in a circular orbit with the same period as its true orbit after the same amount of time since it last crossed the perigee had ellapsed. Therefore, 0 denotes that the object is at the perigee |
meanMotion |
Mean motion of the satellite at epoch in revolutions/day |
meanMotionDerivative |
First time derivative of the mean motion of the satellite in revolutions/day^2^ |
meanMotionSecondDerivative |
Second time derivative of the mean motion of the satellite in revolutions/day^3^. |
Bstar |
Drag coefficient of the satellite in units of (earth radii)^-1^. Bstar is an adjusted value of the ballistic coefficient of the satellite, and it indicates how susceptible it is to atmospheric drag. |
ephemerisType |
Source for the ephemeris (orbital state vector). Most commonly, it is distributed data obtained by combaining multiple observations with the SGP4/SDP4 models |
epochRevolutionNumber |
Number of full orbital revolutions completed by the object |
objectName |
Name of the object, retrieved from the first line of the TLE if a Three Line Element was provided |
https://celestrak.org/columns/v04n03/#FAQ01
# The following lines correspond to a TLE for Italsat 2 the 26th of June, 2006 # at 00:58:29.34 UTC. italsat2_lines <- c("ITALSAT 2", "1 24208U 96044A 06177.04061740 -.00000094 00000-0 10000-3 0 1600", "2 24208 3.8536 80.0121 0026640 311.0977 48.3000 1.00778054 36119") italsat2_TLE <- parseTLElines(italsat2_lines) italsat2_TLE# The following lines correspond to a TLE for Italsat 2 the 26th of June, 2006 # at 00:58:29.34 UTC. italsat2_lines <- c("ITALSAT 2", "1 24208U 96044A 06177.04061740 -.00000094 00000-0 10000-3 0 1600", "2 24208 3.8536 80.0121 0026640 311.0977 48.3000 1.00778054 36119") italsat2_TLE <- parseTLElines(italsat2_lines) italsat2_TLE
This function converts an angle in radians to degrees.
rad2deg(radians)rad2deg(radians)
radians |
Value of the angle in radians. |
The corresponding value of the angle in degrees.
rad2deg(pi)rad2deg(pi)
DAF (Double Precision Array File) is a binary file architecture designed to store arrays of double precision arrays used by SPICE, NAIF toolkit software library. The architecture forms the basis of multiple file formats used to store different data related to astrodynamics, such as SPK, PCK and CK files (https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/req/daf.html).
DAF files provide a generic architecture onto which more specific file formats are implemented. They are organized in records of fixed length (1024 bytes) containing different information. The first record is called the file record, and contains global metadata for the file. This is followed by an optional block comprising any number of comment records. After this, the file consits of sets of summary records, name records and element records. These are structured as blocks of 1 summary record (which contains multiple array summaries, providing metadata about each array), followed by 1 name record (comprising names for the corresponding arrays whose summaries were in the previous summary record) and finally by as many element records as required to store the arrays described in the corresponding summary records. For a detailed description of the DAF architecture, see NAIF documentation (https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/req/daf.html).
This function allows to read any binary file with a DAF architecture in a generic way, without applying any specific formatting to its contents. Note that this will result in just a list of the file global metadata, comments and different arrays (each one with a summary, a name and a set of elements). The number of elements and meaning of each included in each array, as well as the descriptor integers and doubles contained in the summaries, vary for each specific file type. Therefore, reading a DAF file generically is not likely to bring much meaningful information, unless a precise understanding of the specific file format is taken into account later.
readBinDAF(filename)readBinDAF(filename)
filename |
Path to the binary DAF file. |
A list with three elements. The first element, named metadata, corresponds
to the "file record", which is always the first record in a DAF file, and is a
list with the following metadata elements:
fileType |
String of the format DAF/XXXX indicating the specific file format for the DAF file |
numDoublesSummary |
Number of double precision numbers in each array summary |
numIntsSummary |
Number of integers in each array summary |
summaryRecordSizeDoubles |
Size of each array summary in doubles |
numCharsName |
Number of characters in each array name |
description |
String with an internal name or description of the DAF file |
firstSummaryRecNum |
Integer indicating which is the 1st summary record This can be used to infer the number of comment files. For example, a value of 3 indicates that the 3rd record is the 1st summary record. Since the 1st record is always the file record, this means there is 1 comment record |
lastSummaryRecNum |
Integer indicating which is the last summary record |
firstFreeAddress |
Integer indicating the first free address (in bytes) of the file, i.e., the address of the last byte of the last element record plus 1 |
endianString |
String indicating the endianness of the file. This is automatically taken into account when reading the file. It can be either LTL-IEEE (little endianness) or BIG-IEEE (big endianness), and is determined by the architecture of the system where the file was written |
ftpString |
String used to verify integrity of the DAF file. It should be
exactly equal to |
The second element is named comments, and is a character vector where
each element is a line of comments.
The third element is named arrays, and is a nested list where each top-level
element represents one of the arrays stored in the DAF file and its associated
metadata. Each of the top-level elements is itself a list with the following 3
elements:
arrayName |
String with the name of the array |
arraySummary |
A list with the multiple doubles and integers that are stored
in each array summary of summary records and which provide metadata describing
each array. The elements are named as Double1, Double2, ..., DoubleN; Integer1,
Integer2, ..., Integer(M-2) and finally initialArrayAddress and finalArrayAddress.
N and M are respectively the number of doubles and integers in each array summary,
and are given in elements |
arrayElements |
A numeric vector with all the elements of the array. Note that this includes potentially constants, actual data and additional array metadata. Furthermore, the number, order and meaning of the elements differs greatly between different specific types and subtypes of DAF files, and therefore it is hard to extract any meaningful information without knowledge of the internal organization of each array |
https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/req/daf.html
# The file vgr2_jup230.bsp provided with the package includes information for the # Jupiter flyby of Voyager 2 testDAF <- readBinDAF(paste0(path.package("asteRisk"), "/vgr2_jup230.bsp")) testDAF$metadata # The file seems to be of type SPK testDAF$comments length(testDAF$arrays) # It contains a single array# The file vgr2_jup230.bsp provided with the package includes information for the # Jupiter flyby of Voyager 2 testDAF <- readBinDAF(paste0(path.package("asteRisk"), "/vgr2_jup230.bsp")) testDAF$metadata # The file seems to be of type SPK testDAF$comments length(testDAF$arrays) # It contains a single array
SPK (Spacecraft and Planet Kernel) is a binary file format developed by NAIF to store ephemerides (trajectory) of celestial bodies and spacecraft in the Solar System. The file format is based on the DAF architecture (see readBinDAF). A detailed description of the SPK file format can be found in NAIF's documentation (https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/req/spk.html).
Each SPK file contains several segments, with each segment comprising a summary (or descriptor), a name and an array of double precision elements. Each segment is conceptually equivalent to an array in the context of generic DAF files. There are several types of SPK segments defined by NAIF, each identified by an SPK type code currently ranging from 1 to 21 (some intermediate values are not used or not available for general public use). Each segment type provides ephemerides information in a different way. Note that the segments stored in a single SPK file can be of different types. A detailed description of the organization of the arrays for each SPK type can be found at https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/req/spk.html
This function allows to read SPK binary files of all types except 4, 6, 7 and 16. The data will be presented properly formatted and the meaning of each element is assigned during the reading process of the file. It should be noted that this function just reads SPK kernels; it does not provide any evaluation of ephemerides at arbitrary target times. The process of performing such evaluation differs between SPK types. For example, types 2, 3, 14 and 20 require Chebyshev interpolation; typees 8 and 9 require Lagrange interpolation; type 10 requires the aplication of SGP4/SDP4, etc. Nevertheless, it is still possible to obtain direct ephemerides information just by reading the SPK kernels since many of the types often include reference state vectors between which interpolation is applied, but that can be directly used at the epochs corresponding to said reference state vectors.
readBinSPK(filename)readBinSPK(filename)
filename |
Path to the binary SPK file. |
A list with two elements. The first element, named comments, is a character
vector where each element is a line of comments.
The second element is named segments, and is a nested list where each top-level
element represents one of the segments stored in the SPK file and its associated
metadata. Each of the top-level elements is itself a list with the following 3
elements:
segmentName : String with the name of the segment
segmentSummary : A list with the multiple doubles and integers that are stored in each array summary of summary records and which provide metadata describing each array. In the case of SPK files, these are always the following 9 elements:
SPKType : A description of the type of SPK segment
initialEpoch : The initial epoch for the interval for which ephemeris data are provided in the segment, in ephemeris seconds (equivalent to TDB seconds) past Julian year 2000
finalEpoch : The initial epoch for the interval for which ephemeris data are provided in the segment, in ephemeris seconds (equivalent to TDB seconds) past Julian year 2000
targetNAIFCode : The NAIF integer code for the object for which the segment provides ephemerides. For details, see https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/req/naif_ids.html
targetNAIFCode : The NAIF integer code for the central body of the reference frame in which the segment provides ephemerides. For details, see https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/req/naif_ids.html
frameNAIFCode : The NAIF frame code for the reference frame in which the segment provides ephemerides. For details, see https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/req/frames.html
SPKType : The SPK type code for the segment
initialArrayAddress : The initial address of the array elements corresponding to this segment within the SPK file, in double precision numbers (in order to obtain byte address, multiply by 8 and subtract 7)
finalArrayAddress : The final address of the array elements corresponding to this segment within the SPK file, in double precision numbers (in order to obtain byte address, multiply by 8 and subtract 7)
segmentData : A list with the actual ephemeris data contained in the segment, as well as some type-specific additional metadata
The contents of the last element, segmentData, are different for each SPK
type. Here a summary is provided for each one, but for more detailed descriptions
see NAIF's documentation.
For type 1, which provide Modified Difference Arrays (MDAs), a list where each element is one of the records of the segment. Each of these elements contains the following elements:
referenceEpoch |
The reference epoch that should be used when using this MDA to compute a state vector |
referenceEpoch |
The final epoch for which this MDA should be used to compute a state vector |
stepsizeFunctionVector |
A vector of differences between the reference point and the epochs of all the other data points used to fit the interpolation polynomials from which the MDA is derived. This is of length 15 |
referencePosition |
A numeric vector of length 3 containing the X, Y and Z components of the position at the reference epoch, in km |
referencePosition |
A numeric vector of length 3 containing the X, Y and Z components of the velocity at the reference epoch, in km/s |
MDA |
A matrix with 15 rows and 3 columns providing the constant coefficients to interpolate position, velocity and acceleration at a target epoch. The 1st, 2nd and 3rd columns give the coefficients for the X, Y and Z components respectively. The given coefficients basically are the coefficients of the interpolation polynomial when expressed in Modified Divided Differences (MDD) form. For details, see Shampine and Gordon, 1975. Note that even though the matrix will always have 15 rows (corresponding to 15 coefficients for each components, or an interpolation polynomial degree of 14), some of these can have a value of 0 up to the 15th row, effectively leading to an interpolation polynomial of degree lower than 14. |
maxIntegrationOrderP1 |
The maximum order for the interpolation polynomial
that should be applied amongst all 3 components plus 1. For example, if the interpolation
polynomial orders for the X, Y and Z components are 6, 8 and 7 respectively, this
element will have a value of 9 ( |
integrationOrderArray |
A vector of 3 integers indicating the order of the interpolation polynomials that should be applied for the X, Y and Z components of acceleration respectively. |
A brief description is provided here for the meaning of the coefficients, based on the 'SPICE spke01 math' monograph by Robert Werner. MDAs are a modified version of the more standard coefficients obtained through the divided differences method, which represent an interpolation polynomial in Newton form. In order to calculate a position and velocity, first the basis functions must be computed for both position and velocity. Then each of the components of the basis functions (each of an increasing degree) must be multiplied by the corresponding MDA coefficient, and all the resulting terms are summed.
For type 2, which provide Chebyshev coefficients for position only and at equally spaced time steps, a list with the following elements:
polynomialDegree : An integer indicating the order of the interpolation polynomial that should be applied for all the components.
chebyshevCoefficients : A matrix where each row corresponds to an interpolation interval, and with the following columns:
initialEpoch : Initial epoch of the interpolation intervals, in ephemeris (TDB) seconds since J2000
midPoint : Epoch for the midpoint of the interpolation intervals, in ephemeris (TDB) seconds since J2000
intervalRadius : Radius of the interpolation intervals, in seconds)
positionXCoeffi : A set of N columns, with i ranging from 1 to N, providing
the Chebyshev coefficients for the X component of the position. N is the number
of coefficients for each component, which is equal to polynomialDegree + 1
positionYCoeffi : As positionXCoeffi, but for the Y component of position.
positionZCoeffi : As positionXCoeffi, but for the Z component of position.
For type 3, which provide Chebyshev coefficients for position and velocity and at
equally spaced time steps, the same as for type 2, but chebyshevCoefficients
contains the following additional elements:
velocityXCoeffi |
A set of N columns, with i ranging from 1 to N, providing
the Chebyshev coefficients for the X component of the velocity. N is the number
of coefficients for each component, which is equal to |
velocityYCoeffi |
As |
velocityZCoeffi |
As |
For type 5, which provide discrete state vectors to be propagated following the laws of two-body motion, a list with the following elements:
centralBodyGM The GM parameter (gravitational constant) for the central body, in cube kilometers per square seconds
stateVectors A matrix where each row corresponds to a state vector, and with the following columns:
epoch : Epoch of the state vectors intervals, in ephemeris (TDB) seconds since J2000
positionX : X component of the position, in km
positionY : Y component of the position, in km
positionZ : Z component of the position, in km
velocityX : X component of the velocity, in km/s
velocityY : Y component of the velocity, in km/s
velocityZ : Z component of the velocity, in km/s
For type 8, which provide discrete state vectors at equally spaced time steps to which Lagrange interpolation should be applied to obtain state vectors at arbitrary target times, a list with the following elements:
polynomialDegree : The degree of the interpolation polynomial that should be applied
stateVectors : A matrix where each row corresponds to a state vector, and with the following columns:
epoch : Epoch of the state vectors intervals, in ephemeris (TDB) seconds since J2000
positionX : X component of the position, in km
positionY : Y component of the position, in km
positionZ : Z component of the position, in km
velocityX : X component of the velocity, in km/s
velocityY : Y component of the velocity, in km/s
velocityZ : Z component of the velocity, in km/s
For type 9, which provide discrete state vectors at unequally spaced time steps to which Lagrange interpolation should be applied to obtain state vectors at arbitrary target times, the same as for type 8.
For type 10, which provide TLE that should be propagated with SGP4/SDP4, a list with the following elements:
constants : A list with the following constants used by SGP4/SDP4:
J2 : J2 parameter (dimensionless)
J3 : J3 parameter (dimensionless)
J4 : J4 parameter (dimensionless)
sqrtGM : Square root of the GM parameter, where GM is in cube Earth radii per square minutes
highAltBound : High altitude boundary for atmospheric model (in km)
lowAltBound : Low altitude boundary for atmospheric model (in km)
earthRadius : Equatorial radius of Earth (in km)
distUnitsPerRadius : Distance units per Earth radius. This is usually 1. If different than 1, interprete results with caution
TLEs : A matrix where each row corresponds to a TLE, and with the following columns (nutation angles and their rates are not present in some old SPK files, in which case they will have NULL values):
epoch : Epoch of the TLE, in ephemeris (TDB) seconds since J2000
meanMotionDerivative : First derivative of the mean motion in radians/min^2
meanMotionSecondDerivative : Second derivative of the mean motion in radians/min^3
Bstar : Drag coefficient of the satellite in units of (earth radii)^-1. Bstar is an adjusted value of the ballistic coefficient of the satellite, and it indicates how susceptible it is to atmospheric drag.
inclination : Mean orbital inclination of the satellite in radians
ascension : Mean longitude of the ascending node of the satellite at epoch, also known as right ascension of the ascending node, in radians
eccentricity : Mean eccentricity of the orbit of the object
perigeeArgument : Mean argument of the perigee of the object in radians
meanAnomaly : Mean anomaly of the orbit of the object in radians
meanMotion : Mean motion of the satellite at epoch in radians/min
deltaPsi : Obliquity (psi angle) of the nutation at epoch, in radians
deltaEpsilon : Longitude (epsilon angle) of the nutation at epoch, in radians
deltaPsiDerivative : Derivative of the obliquity (psi angle) of the nutation at epoch, in radians/second
deltaEpsilonDerivative : Derivative of the longitude (epsilon angle) of the nutation at epoch, in radians/second
For type 12, which provide discrete state vectors at equally spaced time steps to which Hermite interpolation should be applied to obtain state vectors at arbitrary target times, the same as for type 8, but the list contains an additional element:
windowSize |
The window size that should be applied during interpolation |
For type 13, which provide discrete state vectors at unequally spaced time steps to which Hermite interpolation should be applied to obtain state vectors at arbitrary target times, the same as for type 12.
For type 14, which provide Chebyshev coefficients for position and velocity and at unequally spaced time steps, the same as for type 3.
For type 15, which provide elements for calculation of ephemerides through the application of a precessing conic propagation model, a list with the following elements:
epochPeriapsis |
Epoch of the periapsis passage in ephemeris (TDB) seconds since J2000 |
unitVectorTrajectoryPoleX |
X component of the unit trajectory pole vector, in km |
unitVectorTrajectoryPoleY |
Y component of the unit trajectory pole vector, in km |
unitVectorTrajectoryPoleZ |
Z component of the unit trajectory pole vector, in km |
unitVectorPeriapsisX |
X component of the unit periapsis vector, in km |
unitVectorPeriapsisY |
Y component of the unit periapsis vector, in km |
unitVectorPeriapsisZ |
Z component of the unit periapsis vector, in km |
semiLatusRectum |
Semi-latus rectum, in km |
eccentricity |
Eccentricity of the orbit |
J2ProcessingFlag |
Flag indicating what J2 corrections should be applied when propagating. If 1, regress line of nodes only. If 2, precess line of apsides only. If 3, don't use any corrections. For any other values, regress line of nodes and precess line of apsides |
unitVectorCentralBodyPoleX |
X component of the unit central body pole vector, in km |
unitVectorCentralBodyPoleY |
Y component of the unit central body pole vector, in km |
unitVectorCentralBodyPoleZ |
Z component of the unit central body pole vector, in km |
centralBodyGM |
The GM parameter (gravitational constant) for the central body, in cube kilometers per square seconds |
centralBodyJ2 |
The J2 parameter for the central body (dimensionless) |
centralBodyRadius |
Radius of the central body, in km |
For type 17, which provide equinoctial elements modelling an object following an elliptic orbit with precessing line of nodes and argument of periapse relative to the equatorial frame of a central body, a list with the following elements:
epochPeriapsis |
Epoch of the periapsis passage in ephemeris (TDB) seconds since J2000 |
semiMajorAxis |
Semi-major axis of the orbit, in km |
equinoctialH |
Value of the equinoctial parameter H at epoch |
equinoctialK |
Value of the equinoctial parameter K at epoch |
meanLongitude |
Mean longitude of the orbit at epoch, in radians |
equinoctialP |
Value of the equinoctial parameter P |
equinoctialQ |
Value of the equinoctial parameter Q |
longitudePeriapsisDerivative |
Derivative of the longitude of periapse at epoch (but it is assumed to be constant at other times), in radians/s |
meanLongitudeDerivative |
Derivative of the mean longitude at epoch (but it is assumed to be constant at other times), in radians/s |
longitudeAscendingNodeDerivative |
Derivative of the longitude of the ascending node at epoch, in radians/s |
equatorialPoleRightAscension |
Right ascension of the pole of the orbital reference system relative to the reference frame of the corresponding SPK segment, in radians |
equatorialPoleDeclination |
Declination of the pole of the orbital reference system relative to the reference frame of the corresponding SPK segment, in radians |
For type 18, which provide ephemerides in the format used by ESA on the Mars Express,
Rosetta, SMART-1 and Venus Express missions (although applicable to any other object),
there are 2 different subtypes: subtype 0 and subtype 1. Subtype 0 should be used
to perform sliding-window Hermite interpolation of position and velocity independently.
Subtype 1 should be used to perform sliding-window Lagrange interpolation of
position and velocity independently. In both cases, segmentData is a list
with the following elements:
subTypeCode : Subtype code for the type 18 SPK segment
polynomialDegree : An integer indicating the order of the interpolation polynomial that should be applied for all the components.
interpolationType : Type of the interpolation that should be applied. Hermite for subtype 0, and Lagrange for subtype 1
windowSize : The window size that should be applied during interpolation
meanLongitude : Mean longitude of the orbit at epoch, in radians
stateVectors : A matrix where each row corresponds to a state vector. The columns differ depending on the subtype. The following will always be present:
epoch : Epoch of the state vectors intervals, in ephemeris (TDB) seconds since J2000
positionX : X component of the position, in km
positionY : Y component of the position, in km
positionZ : Z component of the position, in km
The following are present for subtype 0:
firstVelocityX : X component of the first velocity value, in km/s. The first velocity value should be used together with the reference position to interpolate position values
firstVelocityY : Y component of the first velocity value, in km/s
firstVelocityZ : Z component of the first velocity value, in km/s
secondVelocityX : X component of the second velocity value, in km/s. The second velocity value should be used together with the reference acceleration to interpolate velocity values
secondVelocityY : Y component of the second velocity value, in km/s
secondVelocityZ : Z component of the second velocity value, in km/s
accelerationX : X component of the acceleration, in km/s^2
accelerationY : Y component of the acceleration, in km/s^2
accelerationZ : Z component of the acceleration, in km/s^2
The following are present for subtype 1:
velocityX : X component of the velocity, in km/s
velocityY : Y component of the velocity, in km/s
velocityZ : Z component of the velocity, in km/s
For type 19, which provides the same data as type 18 but condensing multiple type 18 segments into a single 19 segments (only possible if all the segments have the same target object, central body and reference frame; additionally, the coverage of the segments must overlap only at endpoints and leave no gaps), a list with the following 2 elements:
boundaryChoiceFlag : A flag indicating which minisegment should be used for the epochs at which 2 minisegments overlap. If 0, the earlier minisegment that ends at that epoch is used. If 1, the later minisegment that begins at that epoch is used
minisegments : A nested list where each top-level element represents a type 19
subsegment (called minisegments in NAIF's documentation), each of which is a list
with the same elements as a segmentData for a type 18 SPK segment, plus
the following 2 additional elements:
intervalStartEpoch : Beginning of the interpolation interval covered by this minisegment, in ephemeris (TDB) seconds since J2000
intervalEndEpoch : End of the interpolation interval covered by this minisegment,
in ephemeris (TDB) seconds since J2000
Furthermore, a third subtype can be found in type 19 minisegments which is not
found for type 18 segments (at least according to NAIF's documentation). This
has subtype code 2, and should be used to perform sliding-window Hermite
interpolation of position and velocity together (note the difference with subtype 0,
where position and velocity are interpolated independently). In this case,
element interpolationType has a value of "Hermite-joint", and the element
describing the minisegment contains the same elements as a type 18 SPK segment of
subtype 1 (and not of subtype 0).
For type 20, which provide which provide Chebyshev coefficients for velocity only and at equally spaced time steps together with a reference position (so that position can be interpolated by integration of velocity), a list with the following elements:
polynomialDegree : An integer indicating the order of the interpolation polynomial that should be applied for all the components.
dScale : Distance scale used for both position and velocity, in km. For example, if dScale has a value of 149597870.7 (the length of an astronomical unit, AU, in km), it means the distance units are AU.
tScale : Time scale used for velocity, in TDB seconds. For example, a value of 1 means we are using velocity directly in TDB seconds. A value of 86400 means the time units of velocity would be TDB Julian days, etc.
chebyshevCoefficients : A matrix where each row corresponds to an interpolation interval, and with the following columns:
initialEpoch : Initial epoch of the interpolation intervals, in ephemeris (TDB) seconds since J2000
midPoint : Epoch for the midpoint of the interpolation intervals, in ephemeris (TDB) seconds since J2000
intervalRadius : Radius of the interpolation intervals, in seconds)
velocityXCoeffi : A set of N columns, with i ranging from 1 to N, providing
the Chebyshev coefficients for the X component of the velocity. N is the number
of coefficients for each component, which is equal to polynomialDegree + 1
velocityYCoeffiAs velocityXCoeffi, but for the Y component of velocity.
velocityZCoeffiAs velocityXCoeffi, but for the Z component of velocity.
midPointPositionX : X component of the reference position, valid at the midpoint of the interpolation interval. Should be used to interpolate position through integration.
midPointPositionYAs midPointPositionX, but for the Y component of velocity.
midPointPositionZAs midPointPositionX, but for the Z component of velocity.
For type 21, which, like type 1, also provide MDAs, the same nested list as for type 1, with 2 differences. Firstly, unlike for type 1, MDAs of type 21 segments are not limited to a maximum of 15 coefficients per component. Second, each top-level element of the nested list contains an additional element:
numberCoefficients |
The number of coefficients provided for the component with the highest interpolation order |
https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/req/spk.html https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/req/naif_ids.html https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/req/frames.html Shampine, L. F. and Gordon, M. K., Computer Solution of Ordinary Differential Equations: The Initial Value Problem, 1975 Robert Werner, SPICE spke01 math, 2022. https://doi.org/10.5270/esa-tyidsbu
# The file vgr2_jup230.bsp provided with the package includes information for the # Jupiter flyby of Voyager 2 testSPK <- readBinSPK(paste0(path.package("asteRisk"), "/vgr2_jup230.bsp")) length(testSPK$segments) # It contains a single segment testSPK$segments[[1]]$segmentSummary$SPKType testSPK$segments[[1]]$segmentSummary$SPKTypeCode # The segment is of type 1, containing Modified Difference Arrays length(testSPK$segments[[1]]$segmentData) # It contains 566 MDAs# The file vgr2_jup230.bsp provided with the package includes information for the # Jupiter flyby of Voyager 2 testSPK <- readBinSPK(paste0(path.package("asteRisk"), "/vgr2_jup230.bsp")) length(testSPK$segments) # It contains a single segment testSPK$segments[[1]]$segmentSummary$SPKType testSPK$segments[[1]]$segmentSummary$SPKTypeCode # The segment is of type 1, containing Modified Difference Arrays length(testSPK$segments[[1]]$segmentData) # It contains 566 MDAs
OEM (Orbital Ephemeris Message) is one of the three standard file formats defined by the CCSDS for transferring spacecraft orbit information. OEM files contain the position and velocity of a given object at multiple times (epochs). They can also contain optionally acceleration values, covariance matrixes that indicate the uncertainty of the provided state vectors and other additional information. This function reads OEM files, retrieving also the optional fields.
readOEM(filename)readOEM(filename)
filename |
Path to the OEM file. |
A list with two elements. The first element, named header, is a list with
the following elements:
OMVersion |
Version of the OEM format used in the file |
creationDate |
Date of creation of the file |
creator |
Individual or organization that generated the file |
The second element is named dataBlocks, and it contains one element for
each ephemeris data block found in the OEM file. Each of these elements is a
list with the following elements that provide information about the ephemerides
of object (note that some elements are not mandatory and therefore might not be
present in all OEM files; in these cases, their value is set to NULL):
objectName |
Name of the object |
objectID |
Object identifier for the object. Frequently, although not always, the identifier has the format YYYY-NNNPPP, where YYYY is the year of launch, NNN is the three-digit serial number specifying the launch number during year YYYY and PPP is a part specifier comprising 1 to 3 capital letters that indicate the part of the object put into space during the launch. |
referenceFrame |
Frame of reference in which ephemerides are provided. |
refFrameEpoch |
Epoch for the frame of reference, for cases where it is not intrinsic to the frame itself, such as TEME. |
centerName |
Name of the center of coordinates. For example, a celestial body, the barycenter of the entire Solar System or even other spacecraft. |
timeSystem |
Time system used for the ephemerides, covariance matrixes and all other time fields of the data block. |
startTime |
Start time of the time span covered by the ephemerides and covariance matrixes in this data block. |
endTime |
End time of the time span covered by the ephemerides and covariance matrixes in this data block. |
usableStartTime |
Start time of the usable time span covered by the ephemerides and covariance matrixes in this data block. |
usableEndTime |
End time of the usable time span covered by the ephemerides and covariance matrixes in this data block. |
interpolationMethod |
Recommended interpolation method to calculate ephemerides at epochs between those directly given in the data block. |
interpolationOrder |
Recommended interpolation degree to calculate ephemerides at epochs between those directly given in the data block. |
mass |
Mass in kg of the object. |
dragArea |
Effective area of the object subjected to drag, in square meters. |
dragCoefficient |
Drag coefficient of the object. |
solarRadArea |
Effective area of the object subjected to solar radiation pressure, in square meters. |
solarRadCoefficient |
Solar radiation pressure coefficient of the object. |
ephemerides |
Data frame with the 7 or 10 columns providing the ephemerides for the object. The 1st column provides the epochs for each ephemeris; columns 2 to 4 provide the X, Y and Z components of position (in km), and columns 5 to 7 provide the X, Y and Z components of velocity (in km/s). Columns 8 to 10 are optional, and if present provide the X, Y and Z components of acceleration (in km/s2) |
covarianceMatrixes |
List where each element is a 3-element list that provides a covariance matrix for this data block. Each of the 3-element lists corresponding to a covariance matrix contains the following elements:
|
https://public.ccsds.org/Pubs/502x0b2c1e2.pdf https://spotthestation.nasa.gov/trajectory_data.cfm
# The file testOEM.txt provided with the package includes ephemerides data # for the ISS publicly available testOEM_ISS <- readOEM(paste0(path.package("asteRisk"), "/testOEM.txt")) testOEM_ISS$header testOEM_ISS$dataBlocks[[1]]$objectName head(testOEM_ISS$dataBlocks[[1]]$ephemerides)# The file testOEM.txt provided with the package includes ephemerides data # for the ISS publicly available testOEM_ISS <- readOEM(paste0(path.package("asteRisk"), "/testOEM.txt")) testOEM_ISS$header testOEM_ISS$dataBlocks[[1]]$objectName head(testOEM_ISS$dataBlocks[[1]]$ephemerides)
TLE (Two-/Three- Line Element) is a standard structured text file format for representing orbital state vectors. This function reads a TLE file containing one or more TLEs. The TLE file can contain either Two Line Elements or Three Line Elements, but all the TLE in a single file must be of the same type. The two lines of a Two Line Element contain all the ephemeris information. The additional line in a Three Line Element is optional, and contains just the satellite name. For a detailed description of the TLE format, see https://celestrak.com/columns/v04n03/#FAQ01.
readTLE(filename, maxTLEs=NULL)readTLE(filename, maxTLEs=NULL)
filename |
Path to the TLE file. Alternatively, an URL pointing to a TLE file. |
maxTLEs |
Maximum number of TLEs to read, starting from the beginning of the file. By default, all TLEs present in the file are read. |
A list with the following elements that define the orbital state vector of the satellite (or, if the file contained multiple TLE, a nested list, where each element of the top level list represents an orbital state vector, and comprises the following elements):
NORADcatalogNumber |
NORAD Catalog Number, also known as Satellite Catalog Number, assigned by United States Space Command to each artificial object orbiting Earth |
classificationLevel |
Classification level of the information for the orbiting object. Can be unclassified, classified, secret or unknown |
internationalDesignator |
International Designator, also known as COSPAR ID, of the object. It consists of the launch year, separated by a hyphen from a three-digit number indicating the launch number for that year and a set of one to three letters indicating the piece for a launch with multiple pieces. |
launchYear |
The launch year of the object |
launchNumber |
The launch number of the object during its launch year |
launchPiece |
The piece for the launch of the object, if it was a launch with multiple pieces |
dateTime |
Date time string to which the orbital state vector corresponds |
elementNumber |
Element number for the object. In principle, every time a new TLE is generated for an object, the element number is incremented, and therefore element numbers could be used to assess if all the TLEs for a certain object are available. However, in practice it is observed that this is not always the case, with some numbers skipped and some numbers repeated. |
inclination |
Mean orbital inclination of the satellite in degrees. This is the angle between the orbital plane of the satellite and the equatorial plane |
ascension |
Mean longitude of the ascending node of the satellite at epoch, also known as right ascension of the ascending node, in degrees. This is the angle between the direction of the ascending node (the point where the satellite crosses the equatorial plane moving north) and the direction of the First Point of Aries (which indicates the location of the vernal equinox) |
eccentricity |
Mean eccentricity of the orbit of the object. Eccentricity is a measurement of how much the orbit deviates from a circular shape, with 0 indicating a perfectly circular orbit and 1 indicating an extreme case of parabolic trajectory |
perigeeArgument |
Mean argument of the perigee of the object in degrees. This is the angle between the direction of the ascending node and the direction of the perigee (the point of the orbit at which the object is closest to the Earth) |
meanAnomaly |
Mean anomaly of the orbit of the object in degrees. This indicates where the satellite is along its orbital path. It is provided as the angle between the direction of the perigee and the hypothetical point where the object would be if it was moving in a circular orbit with the same period as its true orbit after the same amount of time since it last crossed the perigee had ellapsed. Therefore, 0 denotes that the object is at the perigee |
meanMotion |
Mean motion of the satellite at epoch in revolutions/day |
meanMotionDerivative |
First time derivative of the mean motion of the satellite in revolutions/day^2^ |
meanMotionSecondDerivative |
Second time derivative of the mean motion of the satellite in revolutions/day^3^. |
Bstar |
Drag coefficient of the satellite in units of (earth radii)^-1^. Bstar is an adjusted value of the ballistic coefficient of the satellite, and it indicates how susceptible it is to atmospheric drag. |
ephemerisType |
Source for the ephemeris (orbital state vector). Most commonly, it is distributed data obtained by combaining multiple observations with the SGP4/SDP4 models |
epochRevolutionNumber |
Number of full orbital revolutions completed by the object |
objectName |
Name of the object, retrieved from the first line of the TLE if a Three Line Element was provided |
https://celestrak.org/columns/v04n03/#FAQ01 http://www.celestrak.org/publications/aiaa/2006-6753/AIAA-2006-6753.pdf
# The file testTLE.txt provided with the package includes 29 TLE covering a # variety of satellites, extracted from Revisiting Space Track Report #3 test_TLEs <- readTLE(paste0(path.package("asteRisk"), "/testTLE.txt")) test_TLEs# The file testTLE.txt provided with the package includes 29 TLE covering a # variety of satellites, extracted from Revisiting Space Track Report #3 test_TLEs <- readTLE(paste0(path.package("asteRisk"), "/testTLE.txt")) test_TLEs
This function converts a rotation rate in revolutions per day to radians per minute. This conversion is useful since values in TLEs are given in revolutions per day, but the SGP4 and SDP4 propagators require the mean motion to be provided in radians per minute.
revDay2radMin(revPerDay)revDay2radMin(revPerDay)
revPerDay |
Value of the rotation rate in revolutions per day. |
The corresponding value of the rotation rate in radians per minute.
revDay2radMin(2)revDay2radMin(2)
Given an orbital state vector of a satellite, applies the SDP4 model to
propagate its orbit to the desired time point. This allows the calculation of
the position and velocity of the satellite at different times, both before
and after the time corresponding to the known state vector (referred to as
"epoch"). Kepler's equation is solved through fixed-point integration.
The SDP4 model is a modified version of the SGP4 model that takes into account
the secular and periodic perturbations of the Moon and the Sun on the orbit of
the satellite. It also considers the Earth resonance effects on 24-hour
geostationary and 12-hour Molniya orbits.
Thanks to this, the SDP4 model can correctly propagate the orbit of objects
in deep space (with orbital periods larger than 225 minutes, corresponding to
altitudes higher than 5877.5 km). However, it should be noted that SDP4 employs
only simplified drag equations, and lacks corrections for low-perigee orbits.
Therefore, it is recommended to apply the standard SGP4 model (available through
the sgp4 function) for satellites that are not in deep space. This
implementation is based on a previous SDP4 implementation in Julia
(SatelliteToolbox).
sdp4(n0, e0, i0, M0, omega0, OMEGA0, Bstar, initialDateTime, targetTime, keplerAccuracy=10e-12, maxKeplerIterations=10)sdp4(n0, e0, i0, M0, omega0, OMEGA0, Bstar, initialDateTime, targetTime, keplerAccuracy=10e-12, maxKeplerIterations=10)
n0 |
Mean motion of the satellite at epoch in radians/min. |
e0 |
Mean eccentricity of the orbit of the satellite at epoch. Eccentricity ranges from 0 (perfectly circular orbit) to 1 (parabolic trajectory). |
i0 |
Mean orbital inclination of the satellite at epoch in radians. |
M0 |
Mean anomaly of the satellite at epoch. |
omega0 |
Mean argument of perigee of the satellite at epoch. |
OMEGA0 |
Mean longitude of the ascending node of the satellite at epoch. Also known as right ascension of the ascending node. |
Bstar |
Drag coefficient of the satellite in units of (earth radii)^-1^. Bstar is an adjusted value of the ballistic coefficient of the satellite, and it indicates how susceptible it is to atmospheric drag. |
initialDateTime |
Date-time string in UTC indicating the time corresponding to the known state vector of the satellite. Unlike for SGP4, it must be provided in all cases since it is required to calculate Moon and Sun perturbations. |
targetTime |
Time at which the position and velocity of the satellite should be calculated. It can be provided in two different ways: either as a number corresponding to the time in minutes counting from epoch at which the orbit should be propagated, or as a date-time string in UTC. |
keplerAccuracy |
Accuracy to consider Kepler´s equation solved. If two consecutive solutions differ by a value lower than this accuracy, integration is considered to have converged. |
maxKeplerIterations |
Maximum number of iterations after which fixed-point integration of Kepler's equation will stop, even if convergence according to the accuracy criterion has not been reached. |
A list with three elements. The first two elements represent the position and velocity of the satellite at the target time, in the TEME (True Equator, Mean Equinox) frame of reference. Position values are in km, and velocity values are in km/s. Each of these two elements contains three values, corresponding to the X, Y and Z components of position and velocity in this order. The third element indicates the algorithm used to propagate the orbit (sdp4).
https://juliapackages.com/p/satellitetoolbox https://celestrak.org/NORAD/documentation/spacetrk.pdf http://www.celestrak.org/publications/aiaa/2006-6753/AIAA-2006-6753.pdf
# The following orbital parameters correspond to an object with NORAD catalogue # number 24208 (Italsat 2) the 26th of June, 2006 at 00:58:29.34 UTC. n0 <- 1.007781*((2*pi)/(1440)) # Multiplication by 2pi/1440 to convert to radians/min e0 <- 0.002664 # mean eccentricity at epoch i0 <- 3.8536*pi/180 # mean inclination at epoch in radians M0 <- 48.3*pi/180 # mean anomaly at epoch in radians omega0 <- 311.0977*pi/180 # mean argument of perigee at epoch in radians OMEGA0 <- 80.0121*pi/180 # mean longitude of ascending node at epoch in radians Bstar <- 1e-04 # drag coefficient epochDateTime <- "2006-06-26 00:58:29.34" # Calculation of the orbital period 2*pi/n0 # The period is higher than 225 min, and therefore the SDP4 model should be # applied. Furthermore, from the mean motion in revolutions/day, it can be # seen that it is a geostarionary satellite with a 24-hour period. Let´s # calculate and compare the position and velocity of the satellite at epoch # and 1 day later. state_0 <- sdp4(n0=n0, e0=e0, i0=i0, M0=M0, omega0=omega0, OMEGA0=OMEGA0, Bstar=Bstar, initialDateTime=epochDateTime, targetTime=0) state_1day <- sdp4(n0=n0, e0=e0, i0=i0, M0=M0, omega0=omega0, OMEGA0=OMEGA0, Bstar=Bstar, initialDateTime=epochDateTime, targetTime=1440) state_0 state_1day # The position and velocity are very similar after a full day, in accordance # with the geostationary orbit# The following orbital parameters correspond to an object with NORAD catalogue # number 24208 (Italsat 2) the 26th of June, 2006 at 00:58:29.34 UTC. n0 <- 1.007781*((2*pi)/(1440)) # Multiplication by 2pi/1440 to convert to radians/min e0 <- 0.002664 # mean eccentricity at epoch i0 <- 3.8536*pi/180 # mean inclination at epoch in radians M0 <- 48.3*pi/180 # mean anomaly at epoch in radians omega0 <- 311.0977*pi/180 # mean argument of perigee at epoch in radians OMEGA0 <- 80.0121*pi/180 # mean longitude of ascending node at epoch in radians Bstar <- 1e-04 # drag coefficient epochDateTime <- "2006-06-26 00:58:29.34" # Calculation of the orbital period 2*pi/n0 # The period is higher than 225 min, and therefore the SDP4 model should be # applied. Furthermore, from the mean motion in revolutions/day, it can be # seen that it is a geostarionary satellite with a 24-hour period. Let´s # calculate and compare the position and velocity of the satellite at epoch # and 1 day later. state_0 <- sdp4(n0=n0, e0=e0, i0=i0, M0=M0, omega0=omega0, OMEGA0=OMEGA0, Bstar=Bstar, initialDateTime=epochDateTime, targetTime=0) state_1day <- sdp4(n0=n0, e0=e0, i0=i0, M0=M0, omega0=omega0, OMEGA0=OMEGA0, Bstar=Bstar, initialDateTime=epochDateTime, targetTime=1440) state_0 state_1day # The position and velocity are very similar after a full day, in accordance # with the geostationary orbit
Given an orbital state vector of a satellite, applies the SGP4 or SDP4 model to
propagate its orbit to the desired time point, as appropriate depending on the
orbital period. The model will be automatically chosen depending on the
orbital period. For objects in deep space (with orbital periods larger than
225 minutes, equivalent to altitudes higher than 5877.5 km) the SDP4 model will
be applied. For objects near Earth (orbital periods shorter than 225 minutes, or
altitudes lower than 5877.5 km) the SGP4 model will be used. It is not
recommended to apply SGP4 to objects in deep space or SDP4 to objects near
Earth, but this can be forced by calling directly the sgp4 and
sdp4 functions.
sgdp4(n0, e0, i0, M0, omega0, OMEGA0, Bstar, initialDateTime=NULL, targetTime, keplerAccuracy=10e-12, maxKeplerIterations=10)sgdp4(n0, e0, i0, M0, omega0, OMEGA0, Bstar, initialDateTime=NULL, targetTime, keplerAccuracy=10e-12, maxKeplerIterations=10)
n0 |
Mean motion of the satellite at epoch in radians/min. |
e0 |
Mean eccentricity of the orbit of the satellite at epoch. Eccentricity ranges from 0 (perfectly circular orbit) to 1 (parabolic trajectory). |
i0 |
Mean orbital inclination of the satellite at epoch in radians. |
M0 |
Mean anomaly of the satellite at epoch. |
omega0 |
Mean argument of perigee of the satellite at epoch. |
OMEGA0 |
Mean longitude of the ascending node of the satellite at epoch. Also known as right ascension of the ascending node. |
Bstar |
Drag coefficient of the satellite in units of (earth radii)^-1^. Bstar is an adjusted value of the ballistic coefficient of the satellite, and it indicates how susceptible it is to atmospheric drag. |
initialDateTime |
Date-time string in UTC indicating the time
corresponding to the known state vector of the satellite. It must be provided
for objects in deep space, and also for objects near Earth if |
targetTime |
Time at which the position and velocity of the satellite
should be calculated. It can be provided in two different ways: either as
a number corresponding to the time in minutes counting from epoch at which
the orbit should be propagated, or as a date-time string in UTC, in which case
the date-time string for epoch must be provided through |
keplerAccuracy |
Accuracy to consider Kepler´s equation solved. If two consecutive solutions differ by a value lower than this accuracy, integration is considered to have converged. |
maxKeplerIterations |
Maximum number of iterations after which fixed-point integration of Kepler's equation will stop, even if convergence according to the accuracy criterion has not been reached. |
A list with three elements. The first two elements represent the position and velocity of the satellite at the target time, in the TEME (True Equator, Mean Equinox) frame of reference. Position values are in km, and velocity values are in km/s. Each of these two elements contains three values, corresponding to the X, Y and Z components of position and velocity in this order. The third element indicates the algorithm used to propagate the orbit (sgp4 or sdp4).
https://celestrak.org/NORAD/documentation/spacetrk.pdf http://www.celestrak.org/publications/aiaa/2006-6753/AIAA-2006-6753.pdf
# The following orbital parameters correspond to an object with NORAD catalogue # number 24208 (Italsat 2) the 26th of June, 2006 at 00:58:29.34 UTC. n0 <- 1.007781*((2*pi)/(1440)) # Multiplication by 2pi/1440 to convert to radians/min e0 <- 0.002664 # mean eccentricity at epoch i0 <- 3.8536*pi/180 # mean inclination at epoch in radians M0 <- 48.3*pi/180 # mean anomaly at epoch in radians omega0 <- 311.0977*pi/180 # mean argument of perigee at epoch in radians OMEGA0 <- 80.0121*pi/180 # mean longitude of ascending node at epoch in radians Bstar <- 1e-04 # drag coefficient epochDateTime <- "2006-06-26 00:58:29.34" # Calculation of the orbital period 2*pi/n0 # The period is higher than 225 min, and therefore the SDP4 model should be # applied. Let´s calculatethe position and velocity of the satellite 12 hours # after epoch. italsat_12h <- sgdp4(n0=n0, e0=e0, i0=i0, M0=M0, omega0=omega0, OMEGA0=OMEGA0, Bstar=Bstar, initialDateTime=epochDateTime, targetTime=0) italsat_12h$algorithm # The SDP4 model was correctly chosen.# The following orbital parameters correspond to an object with NORAD catalogue # number 24208 (Italsat 2) the 26th of June, 2006 at 00:58:29.34 UTC. n0 <- 1.007781*((2*pi)/(1440)) # Multiplication by 2pi/1440 to convert to radians/min e0 <- 0.002664 # mean eccentricity at epoch i0 <- 3.8536*pi/180 # mean inclination at epoch in radians M0 <- 48.3*pi/180 # mean anomaly at epoch in radians omega0 <- 311.0977*pi/180 # mean argument of perigee at epoch in radians OMEGA0 <- 80.0121*pi/180 # mean longitude of ascending node at epoch in radians Bstar <- 1e-04 # drag coefficient epochDateTime <- "2006-06-26 00:58:29.34" # Calculation of the orbital period 2*pi/n0 # The period is higher than 225 min, and therefore the SDP4 model should be # applied. Let´s calculatethe position and velocity of the satellite 12 hours # after epoch. italsat_12h <- sgdp4(n0=n0, e0=e0, i0=i0, M0=M0, omega0=omega0, OMEGA0=OMEGA0, Bstar=Bstar, initialDateTime=epochDateTime, targetTime=0) italsat_12h$algorithm # The SDP4 model was correctly chosen.
Given an orbital state vector of a satellite, applies the SGP4 model to
propagate its orbit to the desired time point. This allows the calculation of
the position and velocity of the satellite at different times, both before
and after the time corresponding to the known state vector (referred to as
"epoch"). Kepler's equation is solved through fixed-point integration.
The SGP4 model can only accurately propagate the orbit of objects near Earth
(with an orbital period shorter than 225 minutes, corresponding approximately to
an altitude lower than 5877.5 km). For propagation of objects in deep space, the
SDP4 model should be used, available through the sdp4 function.
This implementation is based on the theory and implementation described in
Space Track Report #3, and includes the corrections summarized in Revisiting
Space Track Report #3.
sgp4(n0, e0, i0, M0, omega0, OMEGA0, Bstar, initialDateTime=NULL, targetTime, keplerAccuracy=10e-12, maxKeplerIterations=10)sgp4(n0, e0, i0, M0, omega0, OMEGA0, Bstar, initialDateTime=NULL, targetTime, keplerAccuracy=10e-12, maxKeplerIterations=10)
n0 |
Mean motion of the satellite at epoch in radians/min. |
e0 |
Mean eccentricity of the orbit of the satellite at epoch. Eccentricity ranges from 0 (perfectly circular orbit) to 1 (parabolic trajectory). |
i0 |
Mean orbital inclination of the satellite at epoch in radians. |
M0 |
Mean anomaly of the satellite at epoch. |
omega0 |
Mean argument of perigee of the satellite at epoch. |
OMEGA0 |
Mean longitude of the ascending node of the satellite at epoch. Also known as right ascension of the ascending node. |
Bstar |
Drag coefficient of the satellite in units of (earth radii)^-1^. Bstar is an adjusted value of the ballistic coefficient of the satellite, and it indicates how susceptible it is to atmospheric drag. |
initialDateTime |
Optional date-time string in UTC indicating the time
corresponding to the known state vector of the satellite. It must be provided
if |
targetTime |
Time at which the position and velocity of the satellite
should be calculated. It can be provided in two different ways: either as
a number corresponding to the time in minutes counting from epoch at which
the orbit should be propagated, or as a date-time string in UTC, in which case
the date-time string for epoch must be provided through |
keplerAccuracy |
Accuracy to consider Kepler´s equation solved. If two consecutive solutions differ by a value lower than this accuracy, integration is considered to have converged. |
maxKeplerIterations |
Maximum number of iterations after which fixed-point integration of Kepler's equation will stop, even if convergence according to the accuracy criterion has not been reached. |
A list with three elements. The first two elements represent the position and velocity of the satellite at the target time, in the TEME (True Equator, Mean Equinox) frame of reference. Position values are in km, and velocity values are in km/s. Each of these two elements contains three values, corresponding to the X, Y and Z components of position and velocity in this order. The third element indicates the algorithm used to propagate the orbit (sgp4).
https://celestrak.org/NORAD/documentation/spacetrk.pdf http://www.celestrak.org/publications/aiaa/2006-6753/AIAA-2006-6753.pdf
# The following orbital parameters correspond to an object with NORAD catalogue # number 88888 the 1st of October, 1980 at 23:41:24 UTC. n0 <- 16.05824518*((2*pi)/(1440)) # Multiplication by 2pi/1440 to convert to radians/min e0 <- 0.0086731 # mean eccentricity at epoch i0 <- 72.8435*pi/180 # mean inclination at epoch in radians M0 <- 110.5714*pi/180 # mean anomaly at epoch in radians omega0 <- 52.6988*pi/180 # mean argument of perigee at epoch in radians OMEGA0 <- 115.9689*pi/180 # mean longitude of ascending node at epoch in radians Bstar <- 0.66816e-4 # drag coefficient # Calculation of the orbital period 2*pi/n0 # The period is lower than 225 min, and therefore the SGP4 model is valid. # Let´s calculate the position and velocity of the satellite 40 minutes after # epoch new_state <- sgp4(n0=n0, e0=e0, i0=i0, M0=M0, omega0=omega0, OMEGA0=OMEGA0, Bstar=Bstar,targetTime = 40) new_state# The following orbital parameters correspond to an object with NORAD catalogue # number 88888 the 1st of October, 1980 at 23:41:24 UTC. n0 <- 16.05824518*((2*pi)/(1440)) # Multiplication by 2pi/1440 to convert to radians/min e0 <- 0.0086731 # mean eccentricity at epoch i0 <- 72.8435*pi/180 # mean inclination at epoch in radians M0 <- 110.5714*pi/180 # mean anomaly at epoch in radians omega0 <- 52.6988*pi/180 # mean argument of perigee at epoch in radians OMEGA0 <- 115.9689*pi/180 # mean longitude of ascending node at epoch in radians Bstar <- 0.66816e-4 # drag coefficient # Calculation of the orbital period 2*pi/n0 # The period is lower than 225 min, and therefore the SGP4 model is valid. # Let´s calculate the position and velocity of the satellite 40 minutes after # epoch new_state <- sgp4(n0=n0, e0=e0, i0=i0, M0=M0, omega0=omega0, OMEGA0=OMEGA0, Bstar=Bstar,targetTime = 40) new_state
The TEME (True Equator, Mean Equinox) and GCRF (Geocentric Celestial Reference Frame) are both ECI frames of reference, i.e., Earth-centered inertial coordinate frames, where the origin is placed at the center of mass of Earth and the coordinate frame is fixed with respect to the stars (and therefore not fixed with respect to the Earth surface in its rotation).
The difference between the two resides in the fact that in the GCRF frame, the X-axis and Z-axis are aligned respectively with the mean equinox and rotation axis of Earth at 12:00 Terrestrial Time on the 1st of January, 2000, while in the TEME frame they are aligned with the mean equinox and rotation axis at the time of the corresponding TLE. Due to the change of the direction of the vernal equinox and the rotation axis over time, coordinates in the two frames differ slightly.
The function implement 2 algorithms. The default one applies first a transformation to ITRF using the algorithm implemented by David A. Vallado in the teme2ecef routine, followed by a transformation to GCRF using nutation-precesion-bias and polar motion matrices (transformation for velocity takes into account variations in Earth's rotation speed through length of day data provided by IERS). The second one follows the same procedure implemented by NAIF in the SPICE routine ZZTEME, and is provided for cases where users aim to reproduce SPICE's output.
This function requires the asteRiskData package, which can be installed
by running install.packages('asteRiskData', repos='https://rafael-ayala.github.io/drat/')
TEMEtoGCRF(position_TEME, velocity_TEME, dateTime, SPICEAlgorithm, ephemerisTime)TEMEtoGCRF(position_TEME, velocity_TEME, dateTime, SPICEAlgorithm, ephemerisTime)
position_TEME |
Vector with the X, Y and Z components of the position of an object in TEME frame, in m. |
velocity_TEME |
Vector with the X, Y and Z components of the velocity of an object in TEME frame, in m/s. |
dateTime |
Date-time string with the date and time in UTC corresponding
to the provided position and velocity vectors. This specifies the time for
which the conversion from TEME to GCRF coordinates will be performed. It is
required due to the change in the exact position of the rotation axis of
Earth due to precesion, nutation and polar motion. Either |
SPICEAlgorithm |
Logical indicating if the algorithm implemented in
SPICE's ZZTEME routine should be used instead of the default algorithm. By
default, |
ephemerisTime |
Time in TDB seconds since J2000 for which the conversion
from TEME to GCRF coordinates will be performed. This is also known as
"ephemeris time" in SPICE. Either |
A list with two elements representing the position and velocity of the satellite in the ECEF (Earth Centered, Earth Fixed) frame of reference. Position values are in m, and velocity values are in m/s. Each of the two elements contains three values, corresponding to the X, Y and Z components of position and velocity in this order.
https://celestrak.org/columns/v04n03/#FAQ01
if(requireNamespace("asteRiskData", quietly = TRUE)) { # The following orbital parameters correspond to an object with NORAD catalogue # number 24208 (Italsat 2) the 26th of June, 2006 at 00:58:29.34 UTC. n0 <- 1.007781*((2*pi)/(1440)) # Multiplication by 2pi/1440 to convert to radians/min e0 <- 0.002664 # mean eccentricity at epoch i0 <- 3.8536*pi/180 # mean inclination at epoch in radians M0 <- 48.3*pi/180 # mean anomaly at epoch in radians omega0 <- 311.0977*pi/180 # mean argument of perigee at epoch in radians OMEGA0 <- 80.0121*pi/180 # mean longitude of ascending node at epoch in radians Bstar <- 1e-04 # drag coefficient epochDateTime <- "2006-06-26 00:58:29.34" # Let´s calculate the position and velocity of the satellite 1 day later state_1day_TEME <- sgdp4(n0=n0, e0=e0, i0=i0, M0=M0, omega0=omega0, OMEGA0=OMEGA0, Bstar=Bstar, initialDateTime=epochDateTime, targetTime=1440) # We can now convert the results in TEME frame to GCRF frame, previously # multiplying by 1000 to convert the km output of sgdp4 to m state_1day_GCRF <- TEMEtoGCRF(state_1day_TEME$position*1000, state_1day_TEME$velocity*1000, "2006-06-27 00:58:29.34") }if(requireNamespace("asteRiskData", quietly = TRUE)) { # The following orbital parameters correspond to an object with NORAD catalogue # number 24208 (Italsat 2) the 26th of June, 2006 at 00:58:29.34 UTC. n0 <- 1.007781*((2*pi)/(1440)) # Multiplication by 2pi/1440 to convert to radians/min e0 <- 0.002664 # mean eccentricity at epoch i0 <- 3.8536*pi/180 # mean inclination at epoch in radians M0 <- 48.3*pi/180 # mean anomaly at epoch in radians omega0 <- 311.0977*pi/180 # mean argument of perigee at epoch in radians OMEGA0 <- 80.0121*pi/180 # mean longitude of ascending node at epoch in radians Bstar <- 1e-04 # drag coefficient epochDateTime <- "2006-06-26 00:58:29.34" # Let´s calculate the position and velocity of the satellite 1 day later state_1day_TEME <- sgdp4(n0=n0, e0=e0, i0=i0, M0=M0, omega0=omega0, OMEGA0=OMEGA0, Bstar=Bstar, initialDateTime=epochDateTime, targetTime=1440) # We can now convert the results in TEME frame to GCRF frame, previously # multiplying by 1000 to convert the km output of sgdp4 to m state_1day_GCRF <- TEMEtoGCRF(state_1day_TEME$position*1000, state_1day_TEME$velocity*1000, "2006-06-27 00:58:29.34") }
The TEME (True Equator, Mean Equinox) frame of reference is an Earth-centered inertial coordinate frame, where the origin is placed at the center of mass of Earth and the coordinate frame is fixed with respect to the stars (and therefore not fixed with respect to the Earth surface in its rotation). The coordinates and velocities calculated with the SGP4 and SDP4 models are in the TEME frame of reference. This function converts positions and velocities in TEME to the ITRF (International Terrestrial Reference Frame), which is an ECEF (Earth Centered, Earth Fixed) frame of reference. In the ITRF, the origin is also placed at the center of mass of Earth, but the frame rotates with respect to the stars to remain fixed with respect to the Earth surface as it rotates. The Z-axis extends along the true North as defined by the IERS reference pole, and the X-axis extends towards the intersection between the equator and the Greenwich meridian at any time.
This function requires the asteRiskData package, which can be installed
by running install.packages('asteRiskData', repos='https://rafael-ayala.github.io/drat/')
TEMEtoITRF(position_TEME, velocity_TEME, dateTime)TEMEtoITRF(position_TEME, velocity_TEME, dateTime)
position_TEME |
Vector with the X, Y and Z components of the position of an object in TEME frame, in m. |
velocity_TEME |
Vector with the X, Y and Z components of the velocity of an object in TEME frame, in m/s. |
dateTime |
Date-time string with the date and time in UTC corresponding to the provided position and velocity vectors. This specifies the time for which the conversion from TEME to ITRF coordinates will be performed. It is important to provide an accurate value, since the point over the surface of Earth to which a set of TEME coordinates refers varies with time due to the motion of Earth. |
A list with two elements representing the position and velocity of the satellite in the ITRF (International Terrestrial Reference Frame) frame of reference. Position values are in m, and velocity values are in m/s. Each of the two elements contains three values, corresponding to the X, Y and Z components of position and velocity in this order.
https://celestrak.org/columns/v04n03/#FAQ01
if(requireNamespace("asteRiskData", quietly = TRUE)) { # The following orbital parameters correspond to an object with NORAD catalogue # number 24208 (Italsat 2) the 26th of June, 2006 at 00:58:29.34 UTC. n0 <- 1.007781*((2*pi)/(1440)) # Multiplication by 2pi/1440 to convert to radians/min e0 <- 0.002664 # mean eccentricity at epoch i0 <- 3.8536*pi/180 # mean inclination at epoch in radians M0 <- 48.3*pi/180 # mean anomaly at epoch in radians omega0 <- 311.0977*pi/180 # mean argument of perigee at epoch in radians OMEGA0 <- 80.0121*pi/180 # mean longitude of ascending node at epoch in radians Bstar <- 1e-04 # drag coefficient epochDateTime <- "2006-06-26 00:58:29.34" # Let´s calculate the position and velocity of the satellite 1 day later state_1day_TEME <- sgdp4(n0=n0, e0=e0, i0=i0, M0=M0, omega0=omega0, OMEGA0=OMEGA0, Bstar=Bstar, initialDateTime=epochDateTime, targetTime=1440) # We can now convert the results in TEME frame to ITRF frame, previously # multiplying by 1000 to convert the km output of sgdp4 to m state_1day_ITRF <- TEMEtoITRF(state_1day_TEME$position*1000, state_1day_TEME$velocity*1000, "2006-06-27 00:58:29.34") }if(requireNamespace("asteRiskData", quietly = TRUE)) { # The following orbital parameters correspond to an object with NORAD catalogue # number 24208 (Italsat 2) the 26th of June, 2006 at 00:58:29.34 UTC. n0 <- 1.007781*((2*pi)/(1440)) # Multiplication by 2pi/1440 to convert to radians/min e0 <- 0.002664 # mean eccentricity at epoch i0 <- 3.8536*pi/180 # mean inclination at epoch in radians M0 <- 48.3*pi/180 # mean anomaly at epoch in radians omega0 <- 311.0977*pi/180 # mean argument of perigee at epoch in radians OMEGA0 <- 80.0121*pi/180 # mean longitude of ascending node at epoch in radians Bstar <- 1e-04 # drag coefficient epochDateTime <- "2006-06-26 00:58:29.34" # Let´s calculate the position and velocity of the satellite 1 day later state_1day_TEME <- sgdp4(n0=n0, e0=e0, i0=i0, M0=M0, omega0=omega0, OMEGA0=OMEGA0, Bstar=Bstar, initialDateTime=epochDateTime, targetTime=1440) # We can now convert the results in TEME frame to ITRF frame, previously # multiplying by 1000 to convert the km output of sgdp4 to m state_1day_ITRF <- TEMEtoITRF(state_1day_TEME$position*1000, state_1day_TEME$velocity*1000, "2006-06-27 00:58:29.34") }
The TEME (True Equator, Mean Equinox) frame of reference is an Earth-centered inertial coordinate frame, where the origin is placed at the center of mass of Earth and the coordinate frame is fixed with respect to the stars (and therefore not fixed with respect to the Earth surface in its rotation). The coordinates and velocities calculated with the SGP4 and SDP4 models are in the TEME frame of reference. This function converts position in TEME to geodetic latitude, longitude and altitude, which can be considered to be a non-inertial, Earth-centered frame of reference.
TEMEtoLATLON(position_TEME, dateTime, degreesOutput=TRUE)TEMEtoLATLON(position_TEME, dateTime, degreesOutput=TRUE)
position_TEME |
Vector with the X, Y and Z components of the position of an object in TEME frame, in m. |
dateTime |
Date-time string with the date and time in UTC corresponding to the provided position vector. This specifies the time for which the conversion from TEME to geodetic coordinates will be performed. It is important to provide an accurate value, since the point over the surface of Earth to which a set of TEME coordinates refers varies with time due to the motion of Earth. |
degreesOutput |
Logical indicating if the output should be in sexagesimal
degrees. If |
A vector with three elements, corresponding to the latitude and longitude in degrees (or radians if specified) and the altitude in m.
https://arc.aiaa.org/doi/10.2514/6.2006-6753
if(requireNamespace("asteRiskData", quietly = TRUE)) { # The following orbital parameters correspond to an object with NORAD catalogue # number 24208 (Italsat 2) the 26th of June, 2006 at 00:58:29.34 UTC. n0 <- 1.007781*((2*pi)/(1440)) # Multiplication by 2pi/1440 to convert to radians/min e0 <- 0.002664 # mean eccentricity at epoch i0 <- 3.8536*pi/180 # mean inclination at epoch in radians M0 <- 48.3*pi/180 # mean anomaly at epoch in radians omega0 <- 311.0977*pi/180 # mean argument of perigee at epoch in radians OMEGA0 <- 80.0121*pi/180 # mean longitude of ascending node at epoch in radians Bstar <- 1e-04 # drag coefficient epochDateTime <- "2006-06-26 00:58:29.34" # Let´s calculate the position and velocity of the satellite 1 day later state_1day_TEME <- sgdp4(n0=n0, e0=e0, i0=i0, M0=M0, omega0=omega0, OMEGA0=OMEGA0, Bstar=Bstar, initialDateTime=epochDateTime, targetTime=1440) # We can now convert the results in TEME frame to geodetic latitude, longitude # and altitude, previously multiplying by 1000 to convert the km output of # sgdp4 to m state_1day_geodetic <- TEMEtoLATLON(state_1day_TEME$position*1000, "2006-06-27 00:58:29.34") }if(requireNamespace("asteRiskData", quietly = TRUE)) { # The following orbital parameters correspond to an object with NORAD catalogue # number 24208 (Italsat 2) the 26th of June, 2006 at 00:58:29.34 UTC. n0 <- 1.007781*((2*pi)/(1440)) # Multiplication by 2pi/1440 to convert to radians/min e0 <- 0.002664 # mean eccentricity at epoch i0 <- 3.8536*pi/180 # mean inclination at epoch in radians M0 <- 48.3*pi/180 # mean anomaly at epoch in radians omega0 <- 311.0977*pi/180 # mean argument of perigee at epoch in radians OMEGA0 <- 80.0121*pi/180 # mean longitude of ascending node at epoch in radians Bstar <- 1e-04 # drag coefficient epochDateTime <- "2006-06-26 00:58:29.34" # Let´s calculate the position and velocity of the satellite 1 day later state_1day_TEME <- sgdp4(n0=n0, e0=e0, i0=i0, M0=M0, omega0=omega0, OMEGA0=OMEGA0, Bstar=Bstar, initialDateTime=epochDateTime, targetTime=1440) # We can now convert the results in TEME frame to geodetic latitude, longitude # and altitude, previously multiplying by 1000 to convert the km output of # sgdp4 to m state_1day_geodetic <- TEMEtoLATLON(state_1day_TEME$position*1000, "2006-06-27 00:58:29.34") }