Generation of Orbit Number and Equator Crossing Time and Longitude for EOSDIS

Generation of Orbit Number and Equator Crossing Time and Longitude for EOSDIS

Peter D. Noerdlinger

Applied Research Corp

Jan 4, 1996

Draft 2.0

This document derives and specifies algorithms for determining the orbit number and the times and longitudes of downward equator crossing for the TRMM, AM and PM spacecraft. It also gives a FORTRAN implementation derived from an FDF routine for determining osculating Keplerian orbital elements from a single data point and suggests a way to improve the accuracy by using several points. The algorithm for Keplerian elements will not be needed for TRMM, but may be needed for later spacecraft. This draft will emphasize the case of TRMM, giving some guidelines for AM and PM. It should be remembered that there will also be historical data that will be ingested and translated as necessary. The recommended method involves a search and root-finding procedure in the actual data packets. It needs two kinds of initialization: (I) for each equator crossing, it is very desirable to confine the search to a small box in time. (II) At the very outset, a "grand initialization" is needed, after which later equator crossings can be assumed to occur once each orbital period. There are some ambiguities of (+-) 2p that I hope we can resolve simply to shorten the work.

List of Symbols:

References, with acronym designations:

I will refer to the following texts as necessary: Acronym designation follows the entry.

Fundamentals of Astrodynamics, by R. Bate, D. Mueller and J. White (Dover Publications, 1971). (BMW)

Methods of Orbit Determination, by P.R. Escobal (Wiley, NY 1965) (ESC)

Goddard Trajectory Determination System (GTDS) Mathematical Theory, Revision 1 Edited by A.C. Long et al, Goddard Space Flight Center Code 550, Document FDD/552-89/001 or CSC/TR-89/6001, 1989 (GTDS)

Spacecraft Attitude Determination and Control, Ed. J. R. Wertz (D. Reidel, Holland, 1978) (JRW)

Explanatory Supplement to the Astronomical Almanac (Issued by USNO; Published by University Science Books, Mill Valley, CA) 1993 (ES)

The IERS Standards Ed: D. McCarthy (USNO 1993) (IERS)

Reference Coordinate Systems

The incoming spacecraft data are assumed to be in either True of Date or J2000 Earth Centered Inertial (ECI) coordinates. {Recent change - July 1996: The incoming data might even be in True of Reference, as for TRMM. The incoming data must be transformed to J2000 before being sent on to the Toolkit. }{End of 1996 change}

I assume that we can in general obtain mean or osculating Keplerian elements for each orbit. Appendix "GTDS" gives code to generate these if they are not available. Code is included to go the other way (orbital elements to position and velocity) in order to keep the complete package and to provide reversability). The FORTRAN code there is adapted by P. Noerdlinger from code kindly supplied by Dr. Robert deFazio of GSFC, based on FORMULATION FROM DOCUMENT X-582-76-77(APRIL 1976), PAGE 3-49.

C MATH THEORY OF THE G.T.D.S. (Goddard Trajectory Determination System) (1976) In the 1989 version CSC/TR-89/6001 = FDD/552-89/001 the coresponding derivation is on p. 3-61. The changes made were to convert to MKS, to delete unneeded constants such as the geosynchronous orbit radius, and to correct the Gaussian Earth gravity constant to the more recent (IERS) value.

Choice of longitudinal marker for the orbits:

The following issues have surfaced as influencing a choice:

(a) While TSDIS will use the longitude of lowest latitude to label orbits, the longitudes of maximum or minimum latitude are inappropriate for AM and PM, which are nearly polar, so that small orbit variations induce large changes in longitude near the pole.

(b) The longitudes of upward or downward equator crossings are then the "natural" labels for the relative longitude of an orbit. Both AM and PM are descending in daylight, and most observing is in daylight. This means that if we calculate the downward crossing fairly accurately, inspection of the data allow easy checks in daylight, and also we optimize any search for daylight conditions. If there were an orbit-adjust (burn) within the orbit, we would do much better labelling the orbit by the downcrossing than the up, given that at the end we wish to find mostly daylight data.

Therefore, I select the longitude of the descending equator crossing as the standard "fiducial" or "reference" longitude for an orbit. Furthermore, I pair this with the time of crossing, requiring both data to be saved together. This acts not only as a safety check but it permits subsetting the orbit accurately for ascending, descending, or other fragments. I denote the longitude as lon_DOWN and the time as t_DOWNcross. In UTC ASCII I’ll call it asciiUTC_DOWN.

Availability of SDP Toolkit and other software

Certain time translation services would best be performed using the SDP Toolkit. Thus I recommend that at least the following subset of tools and files be in place where the data are handled:

1 the Time and Date tools and lower level functions.

2. the leap seconds file

3. The CSC tools for transforming between True of Date and J2000; precession, nutation

4. (probably) the UT1/polar motion file (we may attempt to circumvent this req’t )

5. scripts and functions to update the data file(s)

In addition, it may be necessary or desirable to have available additional time transformations, for example to and from the FGDC time and date standard, because the metadata standards follow FGDC and considerable historical data may do so as well.

Orbital Elements:

Orbital elements will have to be captured from incoming data (or constructed by some means for historic data) to carry out the program in this document. It is recommended that these elements be kept with the ephemeris files as well, to provide for later error checking or reworking. Osculating elements pertain at the instant of their definition - they represent the fit of a Keplerian ellipse to the current position and velocity. Mean elements are either averages over some time period (such as an orbit) or are found by fitting an orbital model of some kind. For example, the Brouwer-Lyddane method allows for perturbations due to the distortion of the Earth’s gravitational field by its equatorial bulge, and, in fact higher order perturbations up to J5. As it does not allow for air drag - a generally bigger effect, it is longer to set up, test, and compute with, and as it has singularities near 63 degrees and 0 degrees inclination, we do not use it. For our purposes, we will assume initially that it makes little or no difference which kind we have - Brouwer-Lyddane mean elements, Keplerian mean, or Keplerian osculating. But we will follow the Keplerian nomenclature. If we wanted to propagate for more than three or four orbits it might become significant to distinguish between these different sets of elements.

The Keplerian Orbital Elements are defined on pp. 58-59 of BMW. Note that there is a severe problem using these elements for our problem if e or i is zero. The case i=0 is very unlikely but could happen for geosynchronous satellites (actually, for technical reasons having to do with solar/lunar perturbations, station-keeping for these satellites is easiest if i is not quite 0.0, so we might be OK, but best to check that i is not zero and abort with an error message if it is!!) The case e = 0 is also unlikely - it is not planned for any of our missions, but the software had better check and issue a warning in this case. If we actually got an orbit with zero eccentricity it might still be OK so long as the data provider defined w and n consistently - the only problem with this kind of orbit is that the position of perigee is undefined (and it tends to move around very very fast when e ~0.0). If the provider kept a consistent set of elements for a circular orbit so that w , even if virtually meaningless, was stabilized during the data set, and measured in relation to it, all is OK. The problem with the inclination is more serious and different algorithms would have to be used. This problem will arise only if importing data for geosynchronous spacecraft in unusual cases where the orbit drifted out of specifications.

How to Proceed:

The first thing to do is to compare the values in the TRMM and AM specifications with the BMW list and be sure everything agrees. Check all units and convert to SI and radians as needed. It is of paramount importance to check the definition of the sixth orbital element T = elemOrb[5] - call the source if it is unclear - as we’ll have to translate this time to the time streams we know and love best. By checking I mean the definition of the kind of time, and whether ASCII, float, PB5, or whatever format.

Note that TRMM orbital data of all kinds are defined in True of Date "inertial" coordinates. These coordinates have their Z axis aligned with the Earth’s rotation at the current date (that of the data) and do not rotate with the Earth’s diurnal rotation. Their X axis is towards the Vernal Equinox of date. This is felicitous for our problem. AM and PM will have J2000 ephemerides (further verify this for PM, please), which means that the Earth’s axis at Julian Date 2451545.0 will be used as the Z reference axis and the Vernal Equinox at that date as the X reference axis. Accordingly, for the AM and PM platforms we will have to transform the data to True of Date before working the problem of equator-crossing. Otherwise we’d get the crossing of an equator defined by J2000, not by the Earth’s equator at that moment. It is not a big change, but it is significant.

Accuracy issues:

It is virtually certain that we will not have an analytic answer for lon_DOWN and the time t_DOWNcross. Therefore, we’ll need to define a tolerance within which we’ll calculate the time (and from it, derive the longitude). The spacecraft moves about 7 km/s and the Earth turns, at the equator, about 450 m/s. For AM and PM a time error translates to an equator longitude largely through the 450 m/sec, so to get within (say) 25 meters (a very conservative number) we can converge the time to 0.05 seconds. For TRMM we need to be more careful because its component of motion in the XY plane is ~ 6900 (m/s) * cos(35 degrees) or 5600 m/s, away from which we can subtract Earth motion for a relative speed of ~ 5000 m/s . Let’s pick a target accuracy of 50 meters for TRMM, as its orbit is less well determined than AM’s. Then we need to converge to 0.01 s. In summary, these will be our target values: 0.05 sec for AM or PM, 0.01 sec for TRMM. (Leave aside that the TRMM orbit may not be accurate enough for this number to be meaningful - we just don’t want to worsen the situation by being sloppy). With these accuracies in time, all else falls into place. These are very conservative numbers. Please keep the design flexible, in the sense that we may want to loosen these tolerances by a factor (say) 5 if we find we have a problem with them.

If the file $PGSDAT/CSC/utcpole.dat is absent or out of date our final error in longitude could be up to 0.9 sec time equivalent or ~14 arc seconds in angle (440m distance), which is not so hot, but is "OK" if we’re stuck. I am working on some other trick in case this is a problem.

We will have to test the algorithms with a simulator. I have the ASAP simulator from COSMIC, which is adequate. I will give you a copy and work with you on setting up simulations (or if you have test data - that’s fine). The Toolkit simulator is not adequate for testing the software you’ll write, as it does not incorporate orbit perturbations except as mean motion of the node.

Formally, the problem divides into 5 parts, which I now outline. There are some shortcuts, however, that can enable us to jump to an approximate answer in some cases. The parts are:

1. Find time of downward equator crossing

2. Find Right Ascension of spacecraft at that time

3. Find Earth aspect (gast) at the same time

4. Find longitude of equator crossing from items (2) and (3)

5. Find time of upward equator crossing and number the orbit

I treat these one by one. I assume the reader knows standard mathematics (vectors, matrices, trigonometry) and standard algorithms to solve for roots of real equations, such as the method of bisection. The latter can be found, e.g. in Numerical Recipes by Vetterling, Teukolsky et al, or A First Course in Numerical Analysis by Ralston (Mc Graw Hill, 1965).

1. Find time of Equator Crossing

The most accurate value will be found by using the actual ephemeris - look for the transition from Z < 0 to Z > 0 and solve a linear equation between these 2 times for the time when Z = 0.

This requires looking in each packet and testing Z. It is desirable to have an approximate value first, so that the actual search can be started (say) 30 seconds or 1 minute earlier, and continued until Z changes sign. This approach also provides a checking procedure. Thus I give some methods for an approximate solution. If you just want to look in the packets and find the first transition from Z > 0 to Z < 0, I won’t contest that procedure. In that case proceed to step 2. But I think it would be desirable to at least check against the orbital element method as a "sanity check." Note: For AM and PM, the orbit must be transformed from J2000 to TOD before the testing is done. I think especially in that case, to keep costs down, it is better to isolate the interval first as shown below.

1. Possible Approach using the Representation with P and Q Vectors

There are some reasons to be prepared to use this representation (from BMW). It can provide a check on the other methods, and may be most efficient for historic data. Thus it is included here, but for TRMM, AM and PM you can skip to Sect. (1b). We use the representation of the orbit on pp. 72-73 of BMW (see the figure on p. 57). If we had to implement this work from the orbital elements, we’d use the representation of the vectors P and Q from the matrix on p. 82. The first column of this matrix is the vector P and the second is the vector Q. Note that the values on p. 82 do not require n, the true anomaly, nor T, the time of perigee passage. That is because the orbit is relatively fixed in inertial space (it precesses slowly and the perigee moves slowly - see pp. 157-158), and the vectors P and Q serve to locate the orbit in (X,Y,Z) space. Furthermore, the equations on pp., 72-73 handle all the motion round the orbit.

Having found P and Q, we would need to use Eqs. (2.5-1) and (2.5-4) to find the equator crossing. (The first equation gives the crossing, the second the sense - up or down). The brute force method is to increment n in the range (0,2p) and look for the downward passage of Z through 0. The method of bisection could be used to converge on the exact time. This overall procedure is not very efficient. The value n @ M given in Eq.(4) below would suffice for a first cut.

Let us now consider simple ways to obtain a good approximation, so that we can converge fast.

1b. Simpler Approximation Methods for Initializing the Search

A. Using the time of last equator crossing and the period.

Get the period from

P = 2p/n (1a)

with

n = (m/a3)1/2 1.0 + (1.5)J2(A/a)2(1-e2)-3/2(1-1.5 sin2i) (1b)

where A is the Earth’s equatorial radius 6378137.0 m, a the spacecraft’s semi–major axis in meters, and i the orbital inclination. [see JRW, p. 67]. The value of m, the gravitational constant times Earth mass is 3.9860044x1014 m3/sec2, (IERS Standards) and the oblateness parameter J2 is 0.00108263. e is the orbital eccentricity. Thus you can increment the time of last equator crossing by the value of P and be within a second or so, I would guess. Certainly within 4 or 5 seconds. I think the best implementation would be to select a window like 30 seconds plus or minus and converge to the nearest 0.01 or 0.05 second within that. You need two checks in case of failure:

(i) if Z has the same sign on both sides of the chosen interval you must use the sense of variation to move the interval back or forward and also you should issue a warning.

(ii) if you iterate more than 10 bisections, loosen the time tolerance a factor 4 and try again. Issue a warning.

B. Initialization:

The very first time (in a day’s batch - or if you can keep track between batches - the first time ever!) the algorithm is run, you will have to initialize it. After I give a simple initialization method I will give some more equations that will serve as checks on the data for AM and PM (and maybe for TRMM).

From the diagrams on p. 59 of BMW, we see that the argument of latitude

u = w + n (2)

measures the angle within the orbital plane from the ascending node to the spacecraft (note that BMW deal only with u0 - the value of u at the "epoch" time; we useu as a running variable). Obviously , we are interested in the value

u = p (3)

which identifies the descending node. (This will be approximate, due to apsidal motion, but it is close for our spacecraft. This relation is one that can break down as e-> 0). Thus, all we need to do find the time t_DOWNcross in which leads to a value of n in Eq(2) such that (3) is satisfied. Technically, for an orbit of appreciable eccentricity, one can relate the time t and the true anomaly n only via Kepler’s Equation (BMW p. 185). But we know from Eqs. (1a-b) that M = n (t - T) (BMW p. 185). From the fact that the orbital eccentricity is small, however, n @ M. Thus, we can solve (2) and (3) as p = u = w + n @ w + M with M = n (t - T) , and n taken from Eq.(1b), w and T from the orbital elements provided. Specifically, the solution, then is

M @ ( p - w ) (4)

The solution is

t_DOWNcross @ T + ( p - w )/n (1st estimate) (5)

Again, I would guess that this is within 30 seconds and is thus good for start of a convergence search. (Check with simulations using "ASAP"). Of course, you may have to replace p by -p, 3p, 5p, 7p etc. in Eq.(5) to get equator crossing times in the desired ranges (within the orbit time ranges). I think we should save these times as secTAI93, but the metadata authorities are requiring that we keep time and date in FGDC format. Tej will work on translations - but I do not think that anyone will forbid us to imbed setTAI93 as well. (Be careful with units throughout - in Eq.(4), for example, you need "n" in radians per second to get the times in seconds - but be careful with the time reference base. I recommend using all times in secTAI93 internally.

Having this starting place and then continuing orbit by orbit as indicated, I now assume we have the downward crossing times. I assume that the time in Eq(4). is now replaced by the converged value (see "Accuracy Issues"). (See "method of bisection".)

2. Find Right Ascension of spacecraft at that time

The right ascension of the spacecraft at the time of downward crossing is

fSC = atan2(Y,X) (5)

where (X, Y and Z) are the TOD coordinates at that time, from the ephemeris. (For AM and PM, one must transform the spacecraft position from J2000 into TOD before doing this work. ) Note that fSC @ Wdesc_node ,but the latter is a very approximate value from the orbital elements, while fSC is precise, from Eq.(5). In fact, Wdesc_node @ W + p (or - p).

3. Find Earth aspect (gast) at the same time

With the Toolkit present, you can call PGS_TD_TAItoGAST() to get gast. Then see diagram - but note I show Wdesc_node in place of fSC in it. I intend to supply a backup approximation in case of Toolkit problems.

4. Find longitude of equator crossing from items (2) and (3)

lon_DOWN = fSC - gast (6)

This will be kept in the metadata for each orbit, along with the time. I recommend doing it in degrees. (I will check with others on this)

5. Find Time of Upward Equator Crossing and Number the orbit

This is easy to figure out once we have decided that the orbit starts at the ascending or descending node - we ought to ask FDF how a batch will begin - or we could decide some way we’ll number. We need FDF and maybe help from the database people on this. Late note: Rick Hatfield believes that FDF will begin a day’s data with an ascending node crossing. Then, you can break the orbits at these crossings. You can find the approximate times from

M @ ( 2p - w ) (7)

The solution is

t_UPcross @ T + ( 2p - w )/n (1st estimate) (8)

Again, to get times right at start of orbit you might have to replace the 2p with -2p , 0, 2p , 4p etc. Maybe you can check the range in the incoming data, or check with FDF on this ambiguity. Once again, if we want to break at the most accurate times we have to start with Eq.(8) and iterate with bisection looking for the transition Z<0 to Z>0 to get the exact time. Check with Rick Hatfield as to how accurate we want this.

5. Final Results

We’ll embed a list of the following form in the top level metadata.

Exactly how to format these items should be to be checked with the metadata people say, Janet Hylton.

Appendix I : - Spacecraft Planned Orbital Elements

TRMM AM PM

1) inclination 35 degrees 98.19999 deg 98.145 deg

2) semi-major axis 6.72939e6 meters 7.08693e6 m 7.0775892e6 m

3) eccentricity 0.00053998 0.00128162 0.0012

4) RA of ascend. node 0.0 degrees 255.35597113 deg 298.54 deg

5) arg. of perigee 90.0 degrees 69.08696217 deg 90 deg

6) MEAN anomaly 270.0 degrees 290.91292528 deg 270 deg

7)epoch (UTC) 10/1/97 23:00:00 6/30/98 10:51:28.32 12/1/00 10:51:28.32

Appendix II - FORTRAN code

CFrom rdefazio@pop500.gsfc.nasa.gov Wed Nov 8 12:03:46 1995

C Modified by Peter Noerdlinger

C

DOUBLE PRECISION FUNCTION ANOMLY(IND,ANGLE,ECC,IPRINT)

IMPLICIT DOUBLE PRECISION (A-H,O-Z)

C

C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

C FUNCTION COMPUTES TRUE ANOMALY FROM MEAN ANOMALY OR VICE-VERSA.

C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

C

PARAMETER(XMU=3.9860044D14,DEGPRD= 57.29577951D0,

1 ERAD=6378.137D3, PI = 3.141592653589793D0,

2 TWOPI = 2.0D0 * PI)

C

DOUBLE PRECISION MEAN

C

C IND=0 OR POSITIVE GET TRUE ANOMALY FROM MEAN.

C IND=NEGATIVE GET MEAN ANOMALY FROM TRUE.

C

C EQUATIONS ARE FROM THE ORBITAL FLIGHT HANDBOOK(MARTIN CO., 1963.

C BILL

C BRYANT HAS A COPY. 3/76), VOL I, PAGE III-21, EQN 2-7 AND

C PAGE III-23, EQNS 2-57 AND 2-59.

C

IERR=0

IF(ECC.LT.0.D0.OR.ECC.GE.1.D0) IERR=1

IF(IERR.NE.0) WRITE(IPRINT,1500) IND,ANGLE,ECC

1500 FORMAT(' ANOMLY. BIG ERROR. ECCENTRICITY NOT FOR AN ELLIPSE.'/,

1 ' IND,ANGLE,ECC=',I3,2G20.10/)

C

IF(IBUG.EQ.1) WRITE(6,3000) IND,ANGLE,ECC

3000 FORMAT(' ANOMLY. DEBUG PRINTOUT OF INPUT.'/,

1 ' IND(NEG SAYS GET MEAN FROM TRUE, POS SAYS GET TRUE FROM MEAN'

2 ,')=',I3/,' ANGLE,ECC=',2G20.10/)

ANG1=ANGLE

C

IF(IND.LT.0) GO TO 1000

C

C

C GET TRUE ANONALY FROM MEAN ANOMALY BY WAY OF ECCENTRIC ANOMALY.

MEAN=ANG1

C FIND ECANOM THE SATISFIES MEAN=ECANOM-ECC*DSIN(ECANOM)

ECANOM=MEAN

DO 8 I=1,50

ERROR=MEAN-ECANOM+ECC*DSIN(ECANOM)

IF(DABS(ERROR).LE.1.D-12) GO TO 50

TEMP=ECC*DCOS(ECANOM)

HELPER=.01D0 * TEMP**3

IF(TEMP.GT.0.99D0 .OR. I.GT.10) HELPER=0.D0

DERIV=1.D0-TEMP+HELPER

DEL=ERROR/DERIV

IF(DEL.GT.2.D0) DEL=2.D0

ECANOM=ECANOM+DEL

8 CONTINUE

WRITE(IPRINT,1100) MEAN,ECANOM

1100 FORMAT(' ANOMLY. ERROR. ECANOM LOOP NOT CONVERGED.

1MEAN,ECANOM=',2G20.10/)

50 CONTINUE

Y=DSQRT(1.D0-ECC*ECC) * DSIN(ECANOM)

X=DCOS(ECANOM) - ECC

TRUANM=DMOD(DATAN2(Y,X)+TWOPI,TWOPI)

ANG2=TRUANM

GO TO 2000

C

C

1000 CONTINUE

C GET MEAN ANOMALY FROM TRUE ANOMALY BY WAY OF ECCENTRIC ANOMALY.

TRUANM=ANG1

COSTRU=DCOS(TRUANM)

ECANOM=DACOS( (ECC+COSTRU)/(1.D0+ECC*COSTRU) )

C ECANOM IS BETWEEN 0.D0 AND PI.

C TRUE ANOMALY IS BETWEEN 0.D0 AND TWOPI

C PUT ECANOM ON THE PROPER SIDE OF PI.

IF(TRUANM.GT.PI) ECANOM=TWOPI-ECANOM

MEAN = ECANOM - ECC*DSIN(ECANOM)

ANG2=MEAN

C

C

2000 CONTINUE

IF(ANGLE.LT.0.D0) ANG2=ANG2-TWOPI

ANOMLY=ANG2

IF(IBUG.EQ.1) WRITE(6,3100) ANG2

3100 FORMAT(' ANOMLY. DEBUG PRINTOUT OF OUTPUT.'/,

1 ' ANG2=',G20.10/)

RETURN

END

SUBROUTINE CROSS (X,Y,Z)

IMPLICIT DOUBLE PRECISION(A-H,O-Z)

C

C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

C SUBROUTINE COMPUTES CROSS PRODUCT OF TWO INPUT VECTORS.

C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

C

DIMENSION X(3),Y(3),Z(3)

DATA I /0/

I = 1

ENTRY UCROSS (X,Y,Z)

Z(1) = X(2)*Y(3) - X(3)*Y(2)

Z(2) = X(3)*Y(1) - X(1)*Y(3)

Z(3) = X(1)*Y(2) - X(2)*Y(1)

IF (I .EQ. 1) GO TO 1

CMAG = DSQRT (Z(1)**2 + Z(2)**2 + Z(3)**2)

IF(CMAG.GT.1.0D-10) GO TO 5

GO TO 2

5 Z(1) = Z(1)/CMAG

Z(2) = Z(2)/CMAG

Z(3) = Z(3)/CMAG

GO TO 1

2 Z(1) = 0.0D0

Z(2) = 0.0D0

Z(3) = 0.0D0

1 I = 0

RETURN

END

DOUBLE PRECISION FUNCTION DOT (X,Y)

IMPLICIT DOUBLE PRECISION(A-H,O-Z)

C

C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

C FUNCTION COMPUTES DOT PRODUCT OF TWO INPUT VECTORS.

C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

C

DIMENSION X(3),Y(3)

DOT = X(1)*Y(1) + X(2)*Y(2) + X(3)*Y(3)

RETURN

END

Content-Type: text/plain; charset="us-ascii"

DOUBLE PRECISION FUNCTION VNORM(X,Y)

IMPLICIT DOUBLE PRECISION(A-H,O-Z)

C

C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

C FUNCTION NORMALIZES THE INPUT VECTOR.

C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

C

DIMENSION X(3),Y(3)

R = DSQRT(X(1)**2 +X(2)**2 +X(3)**2)

Y(1) = X(1)/R

Y(2) = X(2)/R

Y(3) = X(3)/R

VNORM = R

RETURN

END

SUBROUTINE TOCART(KEPLER,POS,VEL,IPRINT)

IMPLICIT DOUBLE PRECISION (A-H,O-Z)

PARAMETER(XMU=3.9860044D14,DEGPRD= 57.29577951D0,

1 ERAD=6378.137D0, PI = 3.141592653589793D0,

2 TWOPI = 2.0D0 * PI)

C

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

C SUBROUTINE COMPUTES CARTESIAN VECTORS FROM OSCULATING KEPLERIAN

C ELEMENTS.

C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

C

DOUBLE PRECISION INCL,NODE,MEAN,KEPLER(6),P(3,2),

1RTEMP(2),VTEMP(2)

DOUBLE PRECISION POS(3),VEL(3)

C

C

SMA=KEPLER(1)

ECC=KEPLER(2)

INCL=KEPLER(3)

NODE=KEPLER(4)

ARGP=KEPLER(5)

MEAN=KEPLER(6)

MEAN=DMOD(MEAN+TWOPI,TWOPI)

C BE SURE CONVENTIONS ARE SATISFIED.

IF(INCL.EQ.0.D0) NODE=0.D0

IF(ECC.EQ.0.D0) ARGP=0.D0

C

IBUG=0

IERR=0

IF(ECC.LT.0.D0 .OR. ECC.GE.1.D0) IERR=1

IF(SMA.LE.0.D0) IERR=IERR+2

IF(IERR.EQ.0) GO TO 1250

WRITE(IPRINT,1200) SMA,ECC

1200 FORMAT('0TOCART. BIG ERROR. ELLIPSE NOT INPUT. ARBITRARY ',

1 'R/V OUTPUT. IGNORE RESULTS. SMA,ECC=',G12.5,1X,G9.2)

1250 CONTINUE

IF(IBUG.EQ.1) WRITE(6,1000) XMU,KEPLER

1000 FORMAT('0TOCART. DEBUG PRINTOUT OF INPUT.'/,' XMU=',G23.16/,

1 ' KEPLER=',3G20.10/,9X,3G20.10/)

IF(IERR.NE.0) GO TO 1400

C

C FORMULATION FROM DOCUMENT X-582-76-77(APRIL 1976), PAGE 3-49.

C MATH THEORY OF THE G.T.D.S.

C

c in the 1989 version CSC/TR-89/6001 = FDD/552-89/001 it seems to

c be p. 3-61. P. Noerdlinger

c

TRANOM=ANOMLY(+1,MEAN,ECC,IPRINT)

ECANOM=DACOS( (DCOS(TRANOM)+ECC)/(1.D0+ECC*DCOS(TRANOM)) )

IF(MEAN.GT.PI) ECANOM=TWOPI-ECANOM

C REMINDER TO CJP...CHECK RANGE OF ANGLE VALUES.

SINEC=DSIN(ECANOM)

COSEC=DCOS(ECANOM)

DUM1=DSQRT(1.D0-ECC*ECC)

RTEMP(1) = SMA * (COSEC-ECC)

RTEMP(2) = SMA * DUM1 * SINEC

DUM2=DSQRT(XMU/SMA)/(1.D0-ECC*COSEC)

VTEMP(1)=-DUM2*SINEC

VTEMP(2)=DUM2*DUM1*COSEC

C

COSAP=DCOS(ARGP)

SINAP=DSIN(ARGP)

COSNOD=DCOS(NODE)

SINNOD=DSIN(NODE)

COSINC=DCOS(INCL)

SININC=DSIN(INCL)

P(1,1)= COSAP*COSNOD - SINAP*SINNOD*COSINC

P(2,1)= COSAP*SINNOD + SINAP*COSNOD*COSINC

P(3,1)= SINAP*SININC

P(1,2)=-SINAP*COSNOD - COSAP*SINNOD*COSINC

P(2,2)=-SINAP*SINNOD + COSAP*COSNOD*COSINC

P(3,2)= COSAP*SININC

C

POS(1) = P(1,1)*RTEMP(1) + P(1,2)*RTEMP(2)

POS(2) = P(2,1)*RTEMP(1) + P(2,2)*RTEMP(2)

POS(3) = P(3,1)*RTEMP(1) + P(3,2)*RTEMP(2)

VEL(1) = P(1,1)*VTEMP(1) + P(1,2)*VTEMP(2)

VEL(2) = P(2,1)*VTEMP(1) + P(2,2)*VTEMP(2)

VEL(3) = P(3,1)*VTEMP(1) + P(3,2)*VTEMP(2)

IF(IBUG.EQ.1) WRITE(6,1100) P,RTEMP,VTEMP,POS,VEL

1100 FORMAT('0TOCART. DEBUG PRINTOUT OF OUTPUT.'/,' P=',3G20.10/,

1 4X,3G20.10/,' RTEMP=',2G20.10/,' VTEMP=',2G20.10/,

2 ' POS=',3G20.10/,' VEL=',3G20.10/)

GO TO 9000

C

1400 CONTINUE

DO 1401 I=1,3

POS(I)=10000.D0

1401 VEL(I)=1.D0

C

9000 CONTINUE

RETURN

END

SUBROUTINE TOCLAS(POS,VEL,KEPLER,I1,I2,I3,I4,I5,I6,IPRINT,

XITOBIG)

IMPLICIT DOUBLE PRECISION (A-H,O-Z)

PARAMETER(XMU=3.9860044D14,DEGPRD= 57.29577951D0,

1 ERAD=6378.137D0, PI = 3.141592653589793D0,

2 TWOPI = 2.0D0 * PI)

C

C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

C SUBROUTINE COMPUTES OSCULATING KEPLERIAN ELEMENTS FROM CARTESIAN

C VECTORS.

C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

C

LOGICAL CIRCLE,FLAT

DOUBLE PRECISION ANGMOM(3),INCL,MEAN,NODE

DOUBLE PRECISION POS(3),VEL(3),KEPLER(6)

DOUBLE PRECISION HOLD(6),UNODE(3),TEMP(3)

EQUIVALENCE (HOLD(1),A),(HOLD(2),ECC),(HOLD(3),INCL),

1 (HOLD(4),ARGP),(HOLD(5),NODE),(HOLD(6),MEAN)

ANGLE(S,C)=DMOD(DATAN2(S,C)+TWOPI,TWOPI)

C

IBUG=0

ITOBIG=0

IF(IBUG.EQ.1) WRITE(6,1000) XMU,I1,I2,I3,I4,I5,I6,POS,VEL

1000 FORMAT('2_TOCLAS. DEBUG PRINTOUT OF INPUT.'/,

1 ' XMU=',G17.10,' I1 THRU I6=',6I3/,

2 ' POS=',3G20.10/,' VEL=',3G20.10/)

C

C FORMULATION IS FROM THE DOCUMENT X-582-76-77(APRIL 1976), P 3-55,

C MATH THEORY OF THE G.T.D.S.

C

2345 CONTINUE

RSQ=POS(1)**2 + POS(2)**2 + POS(3)**2

R=DSQRT(RSQ)

VSQ=VEL(1)**2 + VEL(2)**2 + VEL(3)**2

RDOTV=DOT(POS,VEL)

CALL CROSS(POS,VEL,ANGMOM)

H=DSQRT(ANGMOM(1)**2+ANGMOM(2)**2+ANGMOM(3)**2)

P=(RSQ*VSQ-RDOTV*RDOTV)/XMU

C

C SEMI-MAJOR AXIS.

A=1.D0/(2.D0/R-VSQ/XMU)

IF(A.LE.0.D0) GO TO 9000

C

C ECCENTRICITY.

ECC=DSQRT(DABS(1.D0-P/A))

IF(ECC.LT.1.D-7) ECC=0.D0

C 1.D-7 IN THE LINE ABOVE IS NEEDED BECAUSE OF ROUNDING ERROR WHEN

C CONVERTING BACK AND FORTH BETWEEN TOCART AND TOCLAS ROUTINES.

IF(ECC.GE.1.D0) GO TO 9000

C

C INCLINATION.

INCL=DACOS(ANGMOM(3)/H)

IF(INCL.LT.1.D-10) INCL=0.D0

IF(INCL.GT.PI-1.D-10) INCL=PI

C

C FLAGS FOR SPECIAL CASES. USUAL FORMULAS BREAK DOWN FOR THESE CASES.

CIRCLE= ECC.EQ.0.D0

FLAT= INCL.EQ.0.D0 .OR. INCL.EQ.PI

C

C

C NODE.

IF(.NOT.FLAT) NODE=ANGLE(ANGMOM(1),-ANGMOM(2))

IF( FLAT ) NODE=0.D0

C

C MEAN ANOMALY.

IF(CIRCLE) GO TO 200

ESINE=RDOTV/DSQRT(A*XMU)

ECOSE=1.D0-R/A

ECANOM=ANGLE(ESINE,ECOSE)

MEAN=DMOD(ECANOM-ESINE+TWOPI,TWOPI)

IF(DABS(MEAN).LT.1.D-13) MEAN=0.D0

GO TO 210

200 CONTINUE

C ORBIT IS A CIRCLE. ARGUMENT OF PERIGEE WILL BE DEFINED AS

C OCCURRING AT THE NODE CROSSING. HENCE, FIND THE ANGLE BETWEEN

C THE S/C POSITION AND THE NODE CROSSING.

UNODE(1)=DCOS(NODE)

UNODE(2)=DSIN(NODE)

UNODE(3)=0.D0

DUM=VNORM(POS,TEMP)

MEAN=DACOS(DOT(UNODE,TEMP))

C MEAN IS NUMBER OF RADIANS BEFORE OR AFTER PERIGEE. RESOLVE.

CALL CROSS(ANGMOM,UNODE,TEMP)

IF(DOT(POS,TEMP).LT.0.D0) MEAN=TWOPI-MEAN

210 CONTINUE

C

C ARGUMENT OF PERIGEE. NOTE USE OF 'MEAN' CALCULATED ABOVE.

IF(CIRCLE) GO TO 300

IF(FLAT) GO TO 310

TRU=ANOMLY(1,MEAN,ECC,IPRINT)

WPLUSF=ANGLE(POS(3),(POS(2)*ANGMOM(1)-POS(1)*ANGMOM(2))/H)

ARGP=DMOD(WPLUSF-TRU+TWOPI,TWOPI)

GO TO 320

310 CONTINUE

C ORBIT IS A (NON-CIRCULAR)ELLIPSE WITH ZERO INCLINATION. USE THE

C INERTIAL S/C LONGITUDE AND TRUE ANOMALY.

ALONG=ANGLE(POS(2),POS(1))

TRU=ANOMLY(1,MEAN,ECC,IPRINT)

ARGP=DMOD(TWOPI-TRU+ALONG,TWOPI)

GO TO 320

300 CONTINUE

C ORBIT IS A CIRCLE. DEFINE PERIGEE AT NODE CROSSING.

ARGP=0.D0

320 CONTINUE

C

C

GO TO 900

9000 CONTINUE

A=40000.D0

ECC=0.0100D0

INCL=1.D0/DEGPRD

NODE=0.D0

ARGP=30.D0/DEGPRD

MEAN=PI

ITOBIG=1

WRITE(IPRINT,1200)

1200 FORMAT('3_TOCLAS. BIG ERROR. ELLIPSE OR CIRCLE NOT INPUT.',

1 ' ARBITRARILY ASSIGNED ELEMENTS OUTPUT. IGNORE RESULTS.')

900 CONTINUE

C

C

IF(I1.NE.0) KEPLER(1)= A

IF(I2.NE.0) KEPLER(2)= ECC

IF(I3.NE.0) KEPLER(3)= INCL

IF(I4.NE.0) KEPLER(4)= NODE

IF(I5.NE.0) KEPLER(5)= ARGP

IF(I6.NE.0) KEPLER(6)= MEAN

C

C

DO 500 I=3,6

500 HOLD(I)=HOLD(I)*DEGPRD

IF(IBUG.EQ.1) WRITE(6,1100) KEPLER,HOLD

1100 FORMAT('5_TOCLAS. DEBUG PRINTOUT OF OUTPUT.'/,

1 ' KEPLER=',3G20.10/,9X,3G20.10/,

2 ' FULL SET OF KEPLERIAN ELEMENTS(ANGLES IN DEG)='/,

3 2(5X,3G20.10/))

RETURN

END