TABLE OF CONTENTS 1.1.1 Earth Centered Inertial Coordinate System 1.1.2 Equatorial Coordinate System (a, d, r) 1.1.3 Horizon Coordinate System 1.2 Calculation of a Ballistic Missile Trajectory 1.2.1 Iterative Method for Solving Eccentricity and Time-of-Flight 1.2.2 Determination of Missile Position at Arbitrary Time DT 1.3.2 Keplerian Orbit Parameters
1 Orbital Dynamics1.1 Coordinate SystemsObject positions in GENEOS are calculated in three coordinate systems: Earth Centered Inertial (ECI), Equatorial (Lat, Long, Alt), and Local Horizon (Az, El, Range). The Equatorial locations are the ones that are stored in GENEOS and used to carry various engagement calculations through the Modules. 1.1.1 Earth Centered Inertial Coordinate SystemThe ECI coordinate system (see Figure 1.1.1) is typically defined as a Cartesian coordinate system, where the coordinates (position) are defined as the distance from the origin along the three orthogonal (mutually perpendicular) axes. The z axis runs along the Earth's rotational axis pointing North, the x axis points in the direction of the vernal equinox (more on this in a moment), and the y axis completes the right-handed orthogonal system. As seen in Figure 1.1.1, the vernal equinox is an imaginary point in space which lies along the line representing the intersection of the Earth's equatorial plane and the plane of the Earth's orbit around the Sun or the ecliptic. Another way of thinking of the x axis is that it is the line segment pointing from the center of the Earth towards the center of the Sun at the beginning of Spring, when the Sun crosses the Earth's equator moving North. The x axis, therefore, lies in both the equatorial plane and the ecliptic. These three axes defining the Earth-Centered Inertial coordinate system are 'fixed' in space and do not rotate with the Earth.
Figure 1.1.1 – The Earth Centered Inertial (ECI) Coordinate System 1.1.2 Equatorial Coordinate System (a, d, r)The direction of the rotation axis of the Earth remains almost constant and so does the equatorial plane perpendicular to this axis. Therefore, the equatorial plane is a suitable reference plane for a coordinate frame that has to be independent of time and position of the observer. This is the most important coordinate system and it mirrors the latitude and longitude system used on Earth; it uses a geocentric spherical coordinate system. Using the celestial equator and celestial poles, we can describe the position of any celestial object with two coordinate angles, called right ascension (a) aka RA, and declination (d) aka Dec, and the distance (r) to the object from Earth (Fig. 1.1.2). The angles tell us what direction to look to observe a celestial object. RA tells us how far east or west it is, and declination how high it is. This is analogous to latitude and longitude, respectively, used to specify the location of objects on the Earth’s surface. The right ascension (RA) of a celestial object is the angular distance on the celestial sphere measured eastward along the celestial equator from the vernal equinox to the hour circle passing through the celestial object (Fig. 1.1.2). Right ascension is usually given in combination with declination. The declination (Dec) of a celestial object is the angular distance on the celestial sphere measured north or south from the celestial equator along the hour circle passing through a celestial object. Declination is usually given in combination with right ascension. For example, from the 2001 Astronomical Almanac the star Cas, which has equatorial coordinates a = 00h09m15.5s and d = +59°09'29" is very close to the meridian for the vernal equinox, so that its local hour angle is nearly equal to the local sidereal time, that is, tS @ h. A plus sign, or absence of a sign, indicates the object is north of the celestial equator. A negative sign indicates an object south of the celestial equator. Since the hour angle and sidereal time vary with time at a constant rate, it is practical to measure them and right ascension in units of time. In its hour angle on the hour angle dial of the telescope. The right ascension can be found in a catalog, which can then be added to the hour angle to find the sidereal time at the moment of observation. For any other time, the sidereal time can be evaluated by adding the time elapsed since the initial observation. When we need to be more accurate, we have to use a sidereal clock to measure time intervals. The orbital motion of the Earth makes celestial objects seem to move faster than the Sun across the sky, so a sidereal clock runs faster by 3m56.56s a day compared to an ordinary solar clock. Since declination and right ascension are almost (see perturbations below) independent of the position of the observer and the motions of the Earth, they can be used in catalogs and star maps. Right ascension is not measured in degrees but in units of time corresponding to the time required for the celestial sphere to rotate once. For example, suppose there is a star at your zenith with right ascension 2h0m0s (Fig. 1.1.2). Three hours from now, there will be a different object at your local meridian where the right ascension will be 5h0m0s. Just as on Earth we measure east and west from the prime meridian at Greenwich in London, we need a fixed point in space to measure right ascension from. The point chosen is called the vernal equinox (Fig. 1.1.2). The vernal equinox is the ascending node of the ecliptic on the celestial equator, it occurs once each year in late March, marking the official beginning of Spring. In addition, it is the time at which the apparent longitude of the Sun is 0°. Alternately, it is the point of intersection of the ascending node of the ecliptic on the equatorial plane, it is also used to refer to the date on which the Sun passes through this intersection. This point is relatively stable, moving backwards only about 50" per year, and thus the line between the vernal equinox and center of the Earth makes a suitable reference axis in the equatorial plane.
Figure 1.1.2 – Illustration of Astronomical Geometry The formulas for transforming between horizon and equatorial coordinate systems are given in Eqn. 1 and 2, which reference Figure 1.1.2. The equatorial to horizon coordinate transformations are (Eqn. 1)
The horizon to equatorial coordinate transformations are (Eqn. 2)
Example: Horizon to equatorial transformation. Convert horizon coordinates AZ = 283°16'16" and EL = 19°20'04" to equatorial coordinates. The observer is at the Greenwich meridian, 52° N, and GST is 0h24m05s. Solution: Using equations (2) Algorithm Results 1. Convert azimuth to decimal degrees. AZ = 283.271 111 2. Convert altitude to decimal degrees. EL = 19.334 444 3. Find sin d = sin a sin f + cos a cos fcos A. sin d = 0.394 256 4. Find sin-1. d = 23.219 571° 5. Find cos h. cos h = 0.036 062 6. Find cos-1. h' = 87.933 328° 7. Find sin A (if h' positive, h = 360° - h', else h = h'). sin A = -0.973 295 h = 87.933 328° 8. Convert h to hours. h = 5.862 222 h 9. Convert h and d to hour angle form. h = 5h51m44s d = 23°13'10" 10. Convert LST to decimal hours. LST = 0.401 389 h 11. Find a = LST – h (if negative, add 24). a = 18.539 167 h 12. Convert a to hour angle form. a = 18h32m21s 1.1.3 Horizon Coordinate SystemThe algorithm for converting from equatorial to horizon coordinates (Reference Paragraph 1.1.2) is: 1.Convert h to decimal hours. 2.Convert h to degrees. 3.Convert d to decimal degrees. 4.Find sin a. 5.Find sin-1 to find a. 6.Find cos A. 7.Find cos-1 to find A'. 8.Find sin h (if positive, A = 360 – A', else A = A'). 9.Convert a and A to hour angle form.
1.2 Calculation of a Ballistic Missile TrajectoryGiven the following input parameters one can solve for a ballistic trajectory through a series of analytical and numerical steps. First, an orbit fitting these given parameters is calculated, and then the position of the missile at a given time DT can be determined.
In order to define a ballistic missile trajectory one must find the orbit’s eccentricity and subsequently the time-of-flight (ToF) from the launch location to the impact location. However, due to the rotation of the earth in an eastwardly direction, while a missile is in flight, the impact location will be increasing it’s longitude with respect to a static Earth frame. This modification to the impact longitude can be calculated as,
where LngIstat is the impact longitude in the static frame and ToF is given in seconds. One should note that special consideration must be given to trajectories that cross the international dateline (180° long.). Therefore one cannot solve for an orbit’s eccentricity without knowing the ToF (to determine proper impact longitude) and conversely the time-of-flight can’t be determined without knowing the eccentricity of the orbit. In order to solve this scenario a numerical iterative method must be employed.
Figure 1.2-A – Ballistic Missile Trajectory Geometry
1.2.1 Iterative Method for Solving Eccentricity and Time-of-FlightThe algorithm for calculating an orbit’s eccentricity and a corresponding ToF relies on the equating of two earth-centered angles (Y). The first is the Ballistic Separation Angle (YB), this is a measurement of the angular separation of the launch and impact radii from earth center given an eccentricity and apogee. The second is the Earth Center Separation Angle (YECSA) which is an angular measure between the launch and impact locations on the Earth’s surface including a ToF which adjusts the impact longitude as shown previously. The algorithm used begins by choosing an eccentricity (e) value of ½ . From this value the following parameters can be calculated,
where Re is the mean radius of the Earth in meters and me is the Earths gravitational parameter in km^3/sec^2. Now the radii of the launch and impact locations must be calculated using the following relations,
where n is an index for launch or impact (L, I) and Fe is a constant describing the oblateness of the earth. Now using the results of these calculations for a fixed e one can find the angular positions of the launch and impact locations in terms of the true anomaly n. The true anomaly is the angle measured at earth center from the perigee of the orbit to the location.
We now can determine the Ballistic Separation Angle YB as the difference between the sum of the launch and impact true anomalies with respect to an entire orbital rotation of 360 degrees.
This YB will be compared to the angle YECSA to determine whether the chosen eccentricity will produce a matching trajectory. The angle YECSA can be computed directly from the launch and impact locations as well as the ToF for a given eccentricity. In order to calculate the time-of-flight one must follow the series of equations shown below,
where the angles used in the trigonometric functions must be in radians.
Figure 1.2-B – Geometry of the Eccentric Anomaly, E With this ToF that corresponds to a chosen eccentricity, one can now calculate the angle YECSA which will lead to the determination of the accuracy of the eccentricity value chosen. The earth center separation angle is calculated by,
Now that both angles (YECSA ,YB) have been calculated, the angle difference is desired to determine the accuracy of the eccentricity value chosen. The angle difference is,
When this value tends towards zero, past some desired lower limit, the accurate values of eccentricity (e) and Time-of-Flight (ToF) have been found. The following algorithm describes the loop used to solve for (e, ToF) for ballistic trajectories: 1. Set eo=.5, emax=1, emin=0 2. Using given en solve for basic orbital parameters (a, TP, A0, RL, RI,) 3. Now solve for angular parameters (nL, nI, YB, EL, EI) 4. Compute the areas swept by ballistic trajectory (AL, AI, AB) 5. Find Time-of-Flight (ToF) and solve for impact longitude in the static frame (LngIstat) 6. Now compute the earth center separation angle YECSA 7. Determine the angle difference DY 8. IF DY@0 THEN (e=en, ToF=ToFn) ®EXIT; ELSE Continue 9. IF DY<0 THEN emin=e ELSE emax=e 10. en=( emin+ emax)/2 11. Return to Step 2 1.2.2 Determination of Missile Position at Arbitrary Time DTOnce the trajectory has been solved by the above method, the exact position of the missile in equatorial coordinates can be computed, at any time in its trajectory, by using a combination of iterative and algebraic methods. This problem is known as Kepler’s Problem and revolves around the solution to a single equation,
Where M is the mean anomaly and the eccentric anomaly E is a function of time. This equation is obviously transcendental in E and therefore requires an iterative method to solve. A graph of M vs. E is shown in Figure 5.5.1.6.2 which demonstrates the equivalence of M and E for a value of p and/or for an eccentricity of zero which corresponds to a circular orbit. The mean anomaly can also be written for a given time as,
Figure 1.2.2 – Mean Anomaly vs. Eccentric Anomaly
where to is the launch time from perigee (usually set to zero for t=DT). In that case M(to) would be the mean anomaly at launch calculated using the eccentric anomaly at launch given in the previous section. Now for an arbitrary time DT the following algorithm must be employed to find the proper value for E at that time:
This value for E(DT) can now be used to find the appropriate true anomaly n and angle YECSA. The true anomaly for this time can be computed from the following equation,
and now the earth center separation angle at that time is simply a difference in the true anomaly from launch,
This value will be used in the calculation of the equatorial coordinates at an arbitrary time. In addition to this angle the Launch Azimuth (b) angle must also be found. This is the direction in which the missile is initially launched and is measured clockwise from true north. The angle b is found through the following calculation,
With this launch azimuth it is now possible to find the missiles latitude, longitude, and altitude at any arbitrary time DT. The following relations provide each coordinate,
These equatorial coordinates provide the necessary information to fix a ballistic missile position at any time DT from launch. In addition, the velocity of the missile can also be calculated if these coordinates are converted to ECI and the incremental derivatives are computed between two adjacent time steps. 1.3 Orbit Parameters1.3.1 Two Line Element SetsA NORAD two-line element set consists of two 69-character lines of data which can be used together with NORAD's SGP4/SDP4 orbital model to determine the position and velocity of the associated satellite. The only valid characters in a two-line element set are the numbers 0-9, the capital letters A-Z, the period, the space, and the plus and minus signs—no other characters are valid. Of course, not all valid characters can be used in all columns within the element set. Figure 1.3.1-A shows what type of character is valid for each column. Columns with a space or period can have no other character. Columns with an 'N' can have any number 0-9 or, in some cases, a space. Columns with an 'A' can have any character A-Z or a space. The column with a 'C' can only have a character representing the classification of the element set—normally either a 'U' for unclassified data or an 'S' for secret data (of course, only unclassified data are publicly available). Columns with a '+' can have either a plus sign, a minus sign, or a space and columns with a '-' can have either a plus or minus sign (if the rest of the field is not blank).
1 NNNNNC NNNNNAAA NNNNN.NNNNNNNN +.NNNNNNNN +NNNNN-N +NNNNN-N N NNNNN 2 NNNNN NNN.NNNN NNN.NNNN NNNNNNN NNN.NNNN NNN.NNNN NN.NNNNNNNNNNNNNN
Figure 1.3.1-A Two-Line Element Set Format
Further restrictions are placed upon the values in each column as the individual fields of data are defined. Tables 1.3.1-1 and 1.3.1-2 define each of the individual fields for lines 1 and 2, respectively. Many of these bear additional explanation. Table 1.3.1-1. Two-Line Element Set Format Definition, Line 1
Column 1 of each line of the two-line element set indicates the line number (and hence the format) for that line. The next field on each line (fields 1.2 and 2.2) indicates the satellite number—actually, the NORAD Catalog Number—of the object the data is for. The NORAD Catalog Number is a unique identifier assigned by NORAD for each earth-orbiting artificial satellite in their SATCAT (Satellite Catalog). For a valid two-line element set, fields 1.2 and 2.2 must be identical. As mentioned above, field 1.3 indicates the security classification of the data—all publicly available data will have a 'U' in this field to indicate unclassified data. The next three fields—fields 1.4 through 1.6—define the International Designator of the object. This identifier is an additional unique designation assigned by the World Data Center-A for Rockets and Satellites (WDC-A-R&S) in accordance with international treaty (1975 Convention on Registration of Objects Launched into Outer Space). The WDC-A-R&S works together with NORAD and NASA's National Space Science Data Center (NSSDC) in maintaining this registry. Although there have been some changes in format since it was first used back in the late 1950s (see "Space Surveillance" in Satellite Times Volume 4 Number 1), the International Designator indicates the year of the launch (field 1.4 only gives the last two digits), the launch of that year (field 1.5), and the piece of that launch (field 1.6) for each object. These three fields can be left blank, but all must be present if any is. Finally, field 1.6 can be either right or left justified—the latter is preferred. As an aside, there are some significant differences between NORAD's Catalog Number and the International Designator. For example, NORAD assigns a catalog number based upon when the object was first observed, whereas the International Designator is always tied to the original launch. For example, the 81st launch of 1968 carried four payloads into orbit: OV2-5, ERS 21 and 28, and LES 6. Together with the Titan 3C transtage rocket body, these objects were assigned International Designators 1968-081A through E and Catalog Numbers 03428 through 03431. Just this past October, however, NORAD cataloged two additional pieces associated with this launch as Catalog Numbers 25000 and 25001—they have the International Designators 1968-081F and G. The next two fields (fields 1.7 and 1.8) together define the reference time for the element set and are jointly referred to as the epoch. Field 1.7 is the two-digit year (more on this later) and field 1.8 is the day of that year. The epoch defines the time to which all of the time-varying fields in the element set are referenced. While talking about the epoch, this is perhaps a good place to answer the other time-related questions. First, how is the epoch time format interpreted? This question is best answered by using an example. An epoch of 98001.00000000 corresponds to 0000 UT on 1998 January 01—in other words, midnight between 1997 December 31 and 1998 January 01. An epoch of 98000.00000000 would actually correspond to the beginning of 1997 December 31—strange as that might seem. Note that the epoch day starts at UT midnight (not noon) and that all times are measured mean solar rather than sidereal time units (the answer to our third question). How will the Year 2000 be handled? Right now, probably the most common question I'm asked is how years beyond 1999 will be handled in the two-line element sets. Actually, the use of a two-digit year affects both fields 1.4 and 1.7, although the impact of the latter is definitely much more important. To date, I have been unable to secure an official response to this question from US Space Command. Informally, I have been told that there will be no change in the two-line element set format to accommodate the arrival of the Year 2000. However, while the format may not change, the interpretation of the format does. Apparently, US Space Command sees no need to change the two-line element set format yet since no artificial earth satellites existed prior to 1957. By their reasoning, two-digit years from 57-99 correspond to 1957-1999 and those from 00-56 correspond to 2000-2056. Therefore, there is no need to change the format for another 50 years! Unfortunately, this reasoning is severely flawed. While changing the format would involve changing large numbers of software packages to accommodate the format changes, not changing the format does not remove the need to modify this software. Instead of incorporating a format change along with the software modifications, modifications must be made now to change the epoch interpretation and later when the format is finally changed. Too bad that when Aerospace Defense Command proposed going from the old five-line to the current two-line format in November 1972 they didn't recommend a four-digit year (they did, however, at least change from a one-digit to a two-digit year). Field 1.9 represents the first derivative of the mean motion divided by two, in units of revolutions per day2, and field 1.10 represents the second derivative of the mean motion divided by six, in units of revolutions per day3. Together, these two fields give a second-order picture of how the mean motion is changing with time. However, these two fields are not used by the SGP4/SDP4 orbital models (only by the simpler SGP model) and, therefore, serve no real purpose. Field 1.11 represents something called B* (BSTAR), which is an SGP4-type drag coefficient. In aerodynamic theory, every object has a ballistic coefficient, B, that is the product of its coefficient of drag, CD, and its cross-sectional area, A, divided by its mass, m.
B = CD A/m
The ballistic coefficient represents how susceptible an object is to drag—the higher the number, the more susceptible. B* is an adjusted value of B using the reference value of atmospheric density, ρo.
B* = B ro/2
B* has units of (earth radii)-1. Fields 1.10 and 1.11 have a somewhat different format that the other fields. In particular, they use a modified exponential notation with an implied leading decimal point. This convention is inherited from FORTRAN where all such numbers range from 0 to less than 1. The first six columns of each field represent the mantissa and the last two represent the exponent. For example, the value -12345-6 corresponds to -0.12345 × 10-6. Each of these two fields can be blank, corresponding to a value of zero. Field 1.12 represents the ephemeris type (i.e., orbital model) used to generate the data. Spacetrack Report Number 3 suggests the following assignments: 1=SGP, 2=SGP4, 3=SDP4, 4=SGP8, 5=SDP8. However, this value is used for internal analysis only—all distributed element sets have a value of zero and are generated using the SGP4/SDP4 orbital model (as appropriate). Field 1.13 represents the element set number. Normally, this number is incremented each time a new element set is generated. In practice, however, this doesn't always happen. When operations switch between the primary and backup Space Control Centers, sometimes the element set numbers get out of sync, with some numbers being reused and others skipped. Unfortunately, this makes it difficult to tell if you have all the element sets for a particular object. The last column on each line (fields 1.14 and 2.10) represents a modulo-10 checksum of the data on that line. To calculate the checksum, simply add the values of all the numbers on each line—ignoring all letters, spaces, periods, and plus signs—and assigning a value of 1 to all minus signs. The checksum is the last digit of that sum. Although this is a very simple error-checking procedure, it should catch 90 percent of all errors. However, many errors can still sneak through. To eliminate these, all data posted on the Celestial WWW site not only pass the checksum test, but must also pass both format and range-checking tests (as described in this article). Line 2 consists primarily of mean elements calculated using the SGP4/SDP4 orbital model. The definitions for fields 2.3 through 2.8 can be seen in Table 1.3.1-2 below. Fields 2.3, 2.4, 2.6, and 2.7 all have units of degrees and can range from 0 up to 360 degrees—field 2.3 (inclination) only goes up to 180 degrees. The eccentricity (field 2.5) is a unitless value with an assumed leading decimal point. For example, a value of 1234567 corresponds to an eccentricity of 0.1234567. The mean motion (field 2.8) is measured in revolutions per day.
Table 1.3.1-2. Two-Line Element Set Format Definition, Line 2
The final field on line 2, prior to the checksum, is the rev number. Since there are several conventions for determining rev numbers, this field also bears some clarification. In NORAD's convention, a revolution begins when the satellite is at the ascending node of its orbit and a revolution is the period between successive ascending nodes. The period from launch to the first ascending node is considered to be Rev 0 and Rev 1 begins when the first ascending node is reached. Since many element sets are generated with epochs that place the satellite near its ascending node, it is important to note whether the satellite has reached the ascending node when calculating subsequent rev numbers. In general, any number smaller than the maximum field size can be padded with either leading spaces or leading zeros. In other words, an epoch can be represented as either 98001.12345678 or 98 1.12345678 or an inclination can be represented as 28.1234 or 028.1234. Convention uses leading zeros for fields 1.5 and 1.8 and leading spaces elsewhere, but either is valid. Obviously, there are a few limitations with the current two-line format. First and foremost is the need for a four-digit year in fields 1.4 and 1.7. Next, there is a need for a more robust form of error checking—perhaps a 16-bit CRC. Such a checksum could be applied to both lines together, not only detecting errors within the data but also mismatched lines 1 and 2. If such changes were made, it might also be wise to increase the field size for the Catalog Number to six or seven digits to support the eventual cataloging of smaller debris. The International Designator format seems to suffice for the foreseeable future, with a four-digit year, up to 999 launches (the most we've had to date in any one year was 129 in 1984), and up to 13,824 pieces (the record holder today is 1994-029 with 672 pieces). Of course, the cataloging of smaller debris—which we may be unable to correlate with the original launch—still presents potential problems. 1.3.2 Keplerian Orbit ParametersOrbital Elements (refer to Figure 1.3.2 below):
The relationships between Apogee, Perigee and Eccentricity may be expressed in various forms:
Figure 1.3.2 – Keplerian Orbital Parameters
|
