C$Procedure LSRGP C SUBROUTINE LSRGPJ ( ORB0 , DATE , TSTEP , J , LTOP , LSFLAG , & ORB1 ) C C$ Log C C 19-JUN-1990 - Eric Cannell - creation of LSRGP. The Earth gravity C model is taken directly from the RGP C library. All I did was to add luni-solar C perturbations. C 11-JUL-1990 - Eric Cannell - added LSFLAG C 13-Sep-1991 - Bruce Shapiro - Make J a parameter & change name to LSRGPJ C C$ Purpose C C Please read the Restrictions Section is you intend to modify LSRGP. C C LSRGP (Luni-Solar Recurrence Gravity Perturbations) is an analytic tool C for propagating the orbital parameters of an Earth orbiting satellite C under the influence of both Earth gravity and lunar-solar perturbations. C LSRGP specifically deals with secular and long periodic perturbations. C C Although LSRGP can model the Earth gravity up to J29; luni-solar C gravitational attraction is modeled using a simplified version the C luni-solar gravitational potential. The luni-solar model uses the C following index ranges: C C L = 2 C M = 0 , 1 , 2 C H = 0 , 1 , 2 C C Here, it is assumed that: C C If L = 2, then P = 1. C C J = -1 , 0 , 1 C C Since Q = 2P - L, Q = 0. C C Note that both models use similar indices; however, the models are never C mixed together in the source code. LSRGP computes Earth gravity effects C first. Only after that does LSRGP compute luni-solar pertubations. C C LSRGP can be invoked with LSFLAG equal false in order to ignore luni-solar C effects. In such a case, LSRGP is equivalent to RGP. C C LSRGP uses non-singular orbital elements to avoid singularities near the C critical inclination and other problems that are associated with near C circular orbits. C C LSRGP does not model atmospheric drag and solar radiation pressure; C however, these effects can be easily modelled with other routines. C C$ Input_Arguments C C Name Type Dim Units Description C ----------------------------------------------------------------------------- C ORB0 DP 6 km,deg the initial classical orbital elements. C The six elements are (a,e,i,LAN,w,M). C DATE C*(*) 1 --> epoch of ORB0. DATE is used to determine C the geometry of the Sun and Moon. DATE's C format is 'dd-mmm-yyyy hh:mm:ss.ffff'. C DATE can be abbreviated as per TIMETRANS. C TSTEP DP 1 sec length of time step to propagate orbit C ORB0 using gravity model up to LTOP C LTOP I 1 - the maximum L index to consider in C the Earth gravity model. For example, C LTOP = 15 would use up to J15. C LTOP must be in the range of 2..LMAX. C LSFLAG L 1 - if true, luni-solar effects are included. C if false, luni-solar effects are ignored. C J DP (2:29) - earth gravity field zonals J2, J3, ... C C$ Output_Arguments C C Name Type Dim Units Description C ----------------------------------------------------------------------------- C ORB0 DP 6 km,deg the final classical orbital elements. C The six elements are (a,e,i,LAN,w,M). C C$ References C C 1] Ram Bhat's handout titled "Analytical Tool for Propagating Orbital C Parameters". C C 2] Ram Bhat's handout titled "The MTARG Algorithm", which is a superset C of Reference 1 and the Luni-Solar Gravitational Attraction Algorithm. C C 3] R.H. Merson, "The Dynamic Model of PROP, A Computer Program for the C Refinement of Orbital Parameters of Earth Orbiters", Royal Aircraft C Establishment, Technical Report #66255. C C 4] W.M. Kaula, "Theory of Satellite Geodesy". C C$ Restrictions C C 1] LSRGP is only for orbits about the Earth. C C 2] LSRGP has been designed to be fast, not necessarily structured or C modular. There is a high degree of coupling between the equations and C simple variables. Keep this in mind when attempting to modifiy LSRGP. C C The Earth gravity model and the Luni-Solar perturbation model are C contained in separate areas of the source code. This implies that one C can modify either model without having to worry about messing up the C other. C C Extending the Earth model is pretty much a matter of modifying the C Earth model parameter statements. You may run into an overflow problem C with large factorials. In that case, massage the order of computation C to avoid overflow. At some point, of course, the VAX's maximum double C precision value of about 1.7D+38 is going bite you. C C Extending the Luni-Solar model is also a matter of modifying the C Luni-Solar model parameter statements. However, you will have to C add the additional statements to compute the contents of new elements C in the now larger matricies (i.e., in GSC, DGDE, GMOON, GSUN, FSC, C DFDI, FMOON, and FSUN). Also, you may have to add a Do Loop for the C Q index if it ends up having a range of values. Finally, keep in C mind that I have optimized the loops to be as fast as possible by C moving nested loop invariate computations outward as far as possible. C C Please be careful!! C C 3] LSRGP is not valid for geo-stationery satellites. C C 4] The inclination of the satellite cannot be zero. C C 5] LSRGP does not currently include tesseral harmonics. C C 6] 2 <= LTOP <= LMAX for the Earth model. C C$ Parameters C C These parameters are for the Earth Gravity Model. C INTEGER LMAX PARAMETER ( LMAX = 29 ) INTEGER LBTM PARAMETER ( LBTM = 2 ) INTEGER KBTM PARAMETER ( KBTM = 0 ) INTEGER KMAX PARAMETER ( KMAX = LMAX ) INTEGER KBTM1 PARAMETER ( KBTM1 = KBTM + 1 ) INTEGER KNUM PARAMETER ( KNUM = KMAX - KBTM + 1 ) INTEGER KNUM1 PARAMETER ( KNUM1 = KNUM - 1 ) INTEGER LKNUM PARAMETER ( LKNUM = ( LMAX - LBTM + 1 ) * KNUM ) INTEGER LKNUM1 PARAMETER ( LKNUM1 = LMAX * KNUM ) DOUBLE PRECISION GMRTH PARAMETER ( GMRTH = 398600.44807345D0 ) DOUBLE PRECISION RE PARAMETER ( RE = 6378.14D0 ) C These parameters are for the Luni-Solar Gravity Model. INTEGER HLO PARAMETER ( HLO = 0 ) INTEGER HHI PARAMETER ( HHI = 2 ) INTEGER JLO PARAMETER ( JLO = -1 ) INTEGER JHI PARAMETER ( JHI = 1 ) INTEGER LLO PARAMETER ( LLO = 2 ) INTEGER LHI PARAMETER ( LHI = 2 ) INTEGER MLO PARAMETER ( MLO = 0 ) INTEGER MHI PARAMETER ( MHI = 2 ) INTEGER PLO PARAMETER ( PLO = 1 ) INTEGER PHI PARAMETER ( PHI = 1 ) INTEGER QLO PARAMETER ( QLO = 0 ) INTEGER QHI PARAMETER ( QHI = 0 ) DOUBLE PRECISION GMMOON PARAMETER ( GMMOON = 4902.7927809104D0 ) DOUBLE PRECISION GMSUN PARAMETER ( GMSUN = 132712441933.00783456D0 ) DOUBLE PRECISION NMOON PARAMETER ( NMOON = 2.649899946D-6 ) DOUBLE PRECISION NSUN PARAMETER ( NSUN = 1.990968723D-7 ) C These parameters are useful anywhere. DOUBLE PRECISION PI PARAMETER ( PI = 3. 14159 26535 89793 23846 D0 ) DOUBLE PRECISION PI2 PARAMETER ( PI2 = PI / 2D0 ) DOUBLE PRECISION D2R PARAMETER ( D2R = PI / 180.0D0 ) C C$ Declarations_of_Input_and_Output_Arguments C CHARACTER*(*) DATE LOGICAL LSFLAG INTEGER LTOP DOUBLE PRECISION ORB0 ( 6 ) DOUBLE PRECISION ORB1 ( 6 ) DOUBLE PRECISION TSTEP C C$ Declarations_of_Local_Variables C DOUBLE PRECISION A0 DOUBLE PRECISION A1 DOUBLE PRECISION A2 DOUBLE PRECISION AKPCK DOUBLE PRECISION AKPSK DOUBLE PRECISION ALMOON (LLO:LHI) DOUBLE PRECISION ALSUN (LLO:LHI) DOUBLE PRECISION AM DOUBLE PRECISION AM2 DOUBLE PRECISION AM3 DOUBLE PRECISION AS DOUBLE PRECISION AS2 DOUBLE PRECISION AS3 DOUBLE PRECISION BB DOUBLE PRECISION BDC DOUBLE PRECISION CEC DOUBLE PRECISION CI DOUBLE PRECISION CI2 DOUBLE PRECISION CI4 DOUBLE PRECISION CIM DOUBLE PRECISION CIQSI DOUBLE PRECISION CIS DOUBLE PRECISION COSINT (KMAX) DOUBLE PRECISION CW DOUBLE PRECISION CKW DOUBLE PRECISION CWT1 DOUBLE PRECISION CY0M DOUBLE PRECISION CY0S DOUBLE PRECISION DE DOUBLE PRECISION DEBAR DOUBLE PRECISION DEMH DOUBLE PRECISION DEMJ DOUBLE PRECISION DEML DOUBLE PRECISION DEMM DOUBLE PRECISION DEMOON DOUBLE PRECISION DEMP DOUBLE PRECISION DESH DOUBLE PRECISION DESJ DOUBLE PRECISION DESL DOUBLE PRECISION DESM DOUBLE PRECISION DESP DOUBLE PRECISION DESUN DOUBLE PRECISION DFDI (LLO:LHI,MLO:MHI,PLO:PHI) DOUBLE PRECISION DGDE (LLO:LHI,PLO:PHI,QLO:QHI) DOUBLE PRECISION DI DOUBLE PRECISION DIFFOM DOUBLE PRECISION DIFFOS DOUBLE PRECISION DIMH DOUBLE PRECISION DIMJ DOUBLE PRECISION DIML DOUBLE PRECISION DIMM DOUBLE PRECISION DIMOON DOUBLE PRECISION DIMP DOUBLE PRECISION DISH DOUBLE PRECISION DISJ DOUBLE PRECISION DISL DOUBLE PRECISION DISM DOUBLE PRECISION DISP DOUBLE PRECISION DISUN DOUBLE PRECISION DL DOUBLE PRECISION DLBAR DOUBLE PRECISION DLDUM DOUBLE PRECISION DLM (LLO:LHI,MLO:MHI) DOUBLE PRECISION DLMH DOUBLE PRECISION DLMJ DOUBLE PRECISION DLML DOUBLE PRECISION DLMM DOUBLE PRECISION DLMOON DOUBLE PRECISION DLMP DOUBLE PRECISION DLP DOUBLE PRECISION DLSH DOUBLE PRECISION DLSJ DOUBLE PRECISION DLSL DOUBLE PRECISION DLSM DOUBLE PRECISION DLSP DOUBLE PRECISION DLSUN DOUBLE PRECISION DO DOUBLE PRECISION DODUM DOUBLE PRECISION DOMH DOUBLE PRECISION DOMJ DOUBLE PRECISION DOML DOUBLE PRECISION DOMM DOUBLE PRECISION DOMOON DOUBLE PRECISION DOMP DOUBLE PRECISION DOSH DOUBLE PRECISION DOSJ DOUBLE PRECISION DOSL DOUBLE PRECISION DOSM DOUBLE PRECISION DOSP DOUBLE PRECISION DOSUN DOUBLE PRECISION DX DOUBLE PRECISION DXBAR DOUBLE PRECISION DXMH DOUBLE PRECISION DXMJ DOUBLE PRECISION DXML DOUBLE PRECISION DXMM DOUBLE PRECISION DXMOON DOUBLE PRECISION DXMP DOUBLE PRECISION DXP DOUBLE PRECISION DXSH DOUBLE PRECISION DXSJ DOUBLE PRECISION DXSL DOUBLE PRECISION DXSM DOUBLE PRECISION DXSP DOUBLE PRECISION DXSUN DOUBLE PRECISION E0 DOUBLE PRECISION E2 DOUBLE PRECISION EM DOUBLE PRECISION EM2 DOUBLE PRECISION EM3 DOUBLE PRECISION EM4 DOUBLE PRECISION ES DOUBLE PRECISION ES2 DOUBLE PRECISION ES3 DOUBLE PRECISION ES4 DOUBLE PRECISION ESI (KMAX) DOUBLE PRECISION ETA0 DOUBLE PRECISION ETA1 DOUBLE PRECISION F DOUBLE PRECISION F1 DOUBLE PRECISION F1M DOUBLE PRECISION F1S DOUBLE PRECISION F2 DOUBLE PRECISION F2M DOUBLE PRECISION F2S DOUBLE PRECISION FFY DOUBLE PRECISION FMOON (LLO:LHI,MLO:MHI,HLO:HHI) DOUBLE PRECISION FSC (LLO:LHI,MLO:MHI,PLO:PHI) DOUBLE PRECISION FSUN (LLO:LHI,MLO:MHI,HLO:HHI) DOUBLE PRECISION GMMNA2 DOUBLE PRECISION GMMNQS DOUBLE PRECISION GMNICM DOUBLE PRECISION GMOON (LLO:LHI,HLO:HHI,JLO:JHI) DOUBLE PRECISION GMSNA2 DOUBLE PRECISION GMSNQS DOUBLE PRECISION GSC (LLO:LHI,PLO:PHI,QLO:QHI) DOUBLE PRECISION GSCDF DOUBLE PRECISION GSCFSC DOUBLE PRECISION GSNICS DOUBLE PRECISION GSUN (LLO:LHI,HLO:HHI,JLO:JHI) INTEGER H DOUBLE PRECISION I0 DOUBLE PRECISION I1 DOUBLE PRECISION ICM DOUBLE PRECISION ICS INTEGER IDX DOUBLE PRECISION ISM DOUBLE PRECISION ISS INTEGER IV INTEGER JJ DOUBLE PRECISION J2SQR DOUBLE PRECISION JGA DOUBLE PRECISION JGABA INTEGER K INTEGER KTOP DOUBLE PRECISION KWT INTEGER L DOUBLE PRECISION L0 DOUBLE PRECISION L1 DOUBLE PRECISION L2H DOUBLE PRECISION L2HJ DOUBLE PRECISION LDOT1 DOUBLE PRECISION LKBYKK DOUBLE PRECISION LNRPL DOUBLE PRECISION LUNSTT ( 6 ) INTEGER M DOUBLE PRECISION M0 DOUBLE PRECISION MDLM DOUBLE PRECISION MDP DOUBLE PRECISION MM DOUBLE PRECISION MS DOUBLE PRECISION N DOUBLE PRECISION N0 DOUBLE PRECISION N1 DOUBLE PRECISION NA2 DOUBLE PRECISION NA2QSI DOUBLE PRECISION NJ DOUBLE PRECISION NJR DOUBLE PRECISION NRPL (LMAX) DOUBLE PRECISION NRPL4 DOUBLE PRECISION O0 DOUBLE PRECISION O1 DOUBLE PRECISION ODOT1 DOUBLE PRECISION ODOTJ INTEGER P DOUBLE PRECISION PP DOUBLE PRECISION Q DOUBLE PRECISION Q2 DOUBLE PRECISION QE DOUBLE PRECISION QSA DOUBLE PRECISION QSI DOUBLE PRECISION RPL1 DOUBLE PRECISION SAS DOUBLE PRECISION SEC DOUBLE PRECISION SI DOUBLE PRECISION SI2 DOUBLE PRECISION SICI DOUBLE PRECISION SIM DOUBLE PRECISION SIM2 DOUBLE PRECISION SININT (KMAX) DOUBLE PRECISION SIS DOUBLE PRECISION SIS2 DOUBLE PRECISION SUNSTT ( 6 ) DOUBLE PRECISION SW DOUBLE PRECISION SKW DOUBLE PRECISION SWT DOUBLE PRECISION SY0M DOUBLE PRECISION SY0S DOUBLE PRECISION W0 DOUBLE PRECISION WDOT1 DOUBLE PRECISION WDOTJ DOUBLE PRECISION WM DOUBLE PRECISION WS DOUBLE PRECISION WT DOUBLE PRECISION XI0 DOUBLE PRECISION XI1 DOUBLE PRECISION XM DOUBLE PRECISION XM2 DOUBLE PRECISION XM4 DOUBLE PRECISION XM6 DOUBLE PRECISION XM8 DOUBLE PRECISION XM10 DOUBLE PRECISION XM12 DOUBLE PRECISION XM14 DOUBLE PRECISION XS DOUBLE PRECISION XS2 DOUBLE PRECISION XS4 DOUBLE PRECISION XS6 DOUBLE PRECISION XS8 DOUBLE PRECISION XS10 DOUBLE PRECISION XS12 DOUBLE PRECISION XS14 DOUBLE PRECISION Y0M DOUBLE PRECISION Y0S C C$ Save_Statements C DOUBLE PRECISION A (LBTM:LMAX,KBTM:KMAX) SAVE A DOUBLE PRECISION AK (LBTM:LMAX,KBTM:KMAX) SAVE AK DOUBLE PRECISION AKP (LBTM:LMAX,KBTM:KMAX) SAVE AKP DOUBLE PRECISION ALPHA (1:LMAX,KBTM:KMAX) SAVE ALPHA DOUBLE PRECISION B (LBTM:LMAX,KBTM:KMAX) SAVE B DOUBLE PRECISION BETA (LBTM:LMAX,KBTM:KMAX) SAVE BETA DOUBLE PRECISION BK (LBTM:LMAX,KBTM:KMAX) SAVE BK DOUBLE PRECISION CK (LBTM:LMAX,KBTM:KMAX) SAVE CK DOUBLE PRECISION DK (LBTM:LMAX,KBTM:KMAX) SAVE DK DOUBLE PRECISION FACTRL (0:LMAX) SAVE FACTRL LOGICAL FIRST SAVE FIRST DOUBLE PRECISION GAMMA (LBTM:LMAX,KBTM:KMAX) SAVE GAMMA DOUBLE PRECISION J (2:29) C SAVE J DOUBLE PRECISION T (LBTM:LMAX,KBTM:KMAX) SAVE T DOUBLE PRECISION U2KK (0:KMAX) SAVE U2KK DOUBLE PRECISION UK (KBTM:KMAX) SAVE C C$ Data_Statements C DATA A / LKNUM * 0D0 / DATA AK / LKNUM * 0D0 / DATA AKP / LKNUM * 0D0 / DATA ALPHA / LKNUM1 * 0D0 / DATA B / LKNUM * 0D0 / DATA BETA / LKNUM * 0D0 / DATA BK / LKNUM * 0D0 / DATA CK / LKNUM * 0D0 / DATA COSINT / KMAX * 0D0 / DATA DK / LKNUM * 0D0 / DATA FIRST / .TRUE. / DATA GAMMA / LKNUM * 0D0 / C DATA J / 0.10826258D-02,-0.25338975D-05, C & -0.16238211D-05,-0.22963180D-06, C & 0.54309576D-06,-0.35775823D-06, C & -0.20980278D-06,-0.12070701D-06, C & -0.24441360D-06, 0.23257025D-06, C & -0.19259423D-06,-0.22280385D-06, C & 0.11011675D-06,-0.15949534D-07, C & 0.41544562D-07,-0.87909637D-07, C & -0.70704603D-07,-0.45577223D-08, C & -0.17217485D-06, 0.11404873D-09, C & 0.16061855D-07, 0.14036748D-06, C & 0.19771236D-07,-0.10809433D-07, C & -0.59183191D-07,-0.81329176D-07, C & 0.17087477D-06, 0.86549317D-07 C & / DATA SININT / KMAX * 0D0 / DATA T / LKNUM * 0D0 / DATA UK / 1D0 , KNUM1 * 2D0 / C C$ Method C-& C1 Check if LTOP is in range of 2..LMAX. IF ( ( LTOP .LT. 2 ) .OR. ( LTOP .GT. LMAX ) ) THEN WRITE(*,301) LTOP , LMAX 301 FORMAT(//,1X,'LSRGP: LTOP(',I,') not in range of 2..',I2,'.') STOP END IF C1 Compute SAVED variables only once. IF ( FIRST ) THEN C2 Negate first time through flag. FIRST = .FALSE. C2 KTOP is LTOP. KTOP = LTOP C2 Compute factorial array and scaling array U2KK. U2KK is used to C2 compute GAMMA. FACTRL(0) = 1D0 U2KK (0) = -1D0 DO 100 IDX = 1 , LTOP FACTRL(IDX) = DFLOAT(IDX) * FACTRL(IDX-1) 100 CONTINUE DO 101 IDX = 1 , KTOP U2KK (IDX) = -UK(IDX) & / 2D0 ** DFLOAT( 2 * IDX ) 101 CONTINUE C2 Compute identities ALPHA, BETA, and T. Also save the C2 lower subdiagonal of array B (the diagonal is DATA'ed to 1D0). ALPHA(1,0) = 1D0 DO 102 L = LBTM , LTOP B(L,L-1) = 1D0 DO 104 K = KBTM , L-1 ALPHA(L,K) = DFLOAT( ( L + K + 1 ) * ( L - K ) ) & / 2D0 & / DFLOAT( K + 1 ) BETA (L,K) = ALPHA(L-1,K) 104 CONTINUE DO 105 K = KBTM , L IF ( MOD(L-K,2) .EQ. 0 ) THEN IV = ( L - K ) / 2 C3 The following two statements re-arrange computations C3 in order to avoid overflow errors. Note that the DO C3 loop will NOT be executed if K=0 (standard Fortran-77). LKBYKK = 1D0 / FACTRL(K) / FACTRL(K) DO 1051 IDX = L+1 , L+K LKBYKK = LKBYKK * DFLOAT( IDX ) 1051 CONTINUE T(L,K) = DFLOAT( (-1)**IV ) & * FACTRL( L ) & * LKBYKK & / DFLOAT( 2**L ) & / FACTRL(IV) & / FACTRL(L-IV) END IF 105 CONTINUE 102 CONTINUE END IF C1 Compute constants based upon classical keplerian elements. A0 = ORB0(1) E0 = ORB0(2) I0 = ORB0(3) * D2R O0 = ORB0(4) * D2R W0 = ORB0(5) * D2R M0 = ORB0(6) * D2R A2 = A0 * A0 E2 = E0 * E0 PP = A0 * ( 1D0 - E2 ) N0 = DSQRT( GMRTH / A0**3 ) Q = DSQRT( 1D0 - E2 ) Q2 = Q * Q CI = DCOS( I0 ) CI2 = CI * CI CI4 = CI2 * CI2 SI = DSIN( I0 ) SI2 = SI * SI F = SI2 SICI = SI * CI CW = DCOS( W0 ) SW = DSIN( W0 ) RPL1 = RE / PP NRPL(1) = RPL1 C1 Update N0 with with J2 term. N = N0 * ( 1.0D0 & + 1.5D0 & * J(2) & * Q & * RPL1 * RPL1 & * ( 1D0 - 1.5D0 * SI2 ) & ) DO 106 IDX = 2 , LTOP NRPL(IDX ) = NRPL(IDX-1) * RPL1 NRPL(IDX-1) = NRPL(IDX-1) * N 106 CONTINUE NRPL(LTOP) = NRPL(LTOP) * N NRPL4 = N * RPL1 * RPL1 * RPL1 * RPL1 ESI(1) = E0 * SI DO 107 IDX = 2 , KTOP ESI(IDX) = ESI(IDX-1) * ESI(1) 107 CONTINUE J2SQR = J(2) * J(2) NJR = (3D0/32D0) * J2SQR * NRPL4 C1 Compute non-singular orbital elements. XI0 = E0 * CW ETA0 = E0 * SW L0 = W0 + M0 C1 ********************************************************************* C1 ********************************************************************* C1 *** *** C1 *** Earth Gravity Model *** C1 *** *** C1 ********************************************************************* C1 ********************************************************************* C1 Many factors are a function of the index L. Nest their computation C1 with a single L loop in order to avoid duplicating overhead for C1 for separate L do loops. ODOTJ = 0D0 WDOTJ = 0D0 NJ = 0D0 DO 108 L = LBTM , LTOP C2 Compute GAMMA(L,K). LNRPL = FACTRL(L-1) * NRPL(L) DO 109 K = KBTM , L-2 GAMMA(L,K) = U2KK(K) * LNRPL * T(L,K) / FACTRL(L-K-1) 109 CONTINUE C2 Compute inclination and eccentricity functions, A(L,K) AND B(L,K). A(L,L ) = 1D0 A(L,L-1) = CI DO 110 K = L-2 , KBTM , -1 A(L,K) = CI * A(L,K+1) - ALPHA(L,K+1) * SI2 * A(L,K+2) & / ( 2D0 * (K+1) ) B(L,K) = B(L,K+1) + BETA(L,K+1) * E2 * B(L,K+2) & / ( 2D0 * (K+1) ) 110 CONTINUE C2 Compute secular functions dLAN/dt, dw/dt, and N (ODOTJ, WDOTJ, and NJ). JGABA = J(L) * GAMMA(L,0) * ALPHA(L,0) * B(L,0) * A(L,1) JGA = J(L) * GAMMA(L,0) * A(L,0) BB = B(L,1) * BETA(L,0) ODOTJ = ODOTJ - JGABA WDOTJ = WDOTJ + JGABA * CI + JGA * ( L * B(L,0) + BB ) NJ = NJ + Q * JGA * ( ( L - 1 ) * B(L,0) - BB ) 108 CONTINUE C1 Now add to secular variations the J2**2 perturbations. ODOT1 = ODOTJ & + NJR * CI & * ( 4D0 + 12D0 * Q - 9D0 * E2 & - ( 40D0 + 36D0 * Q - 5D0 * E2 ) * CI2 & ) WDOT1 = WDOTJ & + -0.25D0 * NJR & * ( 10D0 - 24D0 * Q + 25D0 * E2 & + ( 36D0 + 192D0 * Q - 126D0 * E2 ) * CI2 & - ( 430D0 + 360D0 * Q - 45D0 * E2 ) * CI4 & ) N1 = NJ & + 0.25D0 * NJR * Q & * ( ( 10D0 + 16D0 * Q - 25D0 * E2 ) & - ( 60D0 + 96D0 * Q - 90D0 * E2 ) * CI2 & + ( 130D0 + 144D0 * Q - 25D0 * E2 ) * CI4 & ) LDOT1 = N1 + WDOT1 C1 Compute the integrals for sine and cosine. DO 111 K = 1 , KTOP KWT = K * WDOT1 * TSTEP F1 = DSIN( KWT ) / KWT F2 = ( 1D0 - DCOS( KWT ) ) / KWT**2 CKW = DCOS( K * ( W0 - PI2 ) ) SKW = DSIN( K * ( W0 - PI2 ) ) COSINT(K) = TSTEP * ( F1 * CKW - KWT * F2 * SKW ) SININT(K) = TSTEP * ( F1 * SKW + KWT * F2 * CKW ) 111 CONTINUE C1 Once again, nest computations with a single L loop in order to avoid C1 duplicating overhead for for separate L do loops. AKPCK = 0D0 AKPSK = 0D0 BDC = 0D0 DODUM = 0D0 DLDUM = 0D0 DO 112 L = LBTM + 1 , LTOP C2 In fact, we can use a single K loop also! As it turns out, GAMMA(L,K) C2 is a factor of each coefficient and series inside this K loop. C2 Since GAMMA(L,K) = 0 when L - K is odd, only execute the body of the C2 K loop when L - K is even. Of course, one must set the skipped C2 elements to zero in a DATA statement. SAVEing them is probably C2 a good idea also. DO 113 K = KBTM1 , L-2 IF ( MOD( L - K , 2 ) .EQ. 0 ) THEN C3 Compute long periodic coefficients. AK (L,K) = J(L) * GAMMA(L,K) * ESI(K) * A(L,K) * B(L,K) AKP(L,K) = K * AK(L,K) / ESI(1) BK (L,K) = L * AK(L,K) CK (L,K) = AK(L,K) * A(L,K+1) * ALPHA(L,K) / A(L,K) DK (L,K) = AK(L,K) * B(L,K+1) * BETA (L,K) / B(L,K) C3 Compute series for long periodics. AKPCK = AKPCK + AKP(L,K) * COSINT(K) AKPSK = AKPSK + AKP(L,K) * SININT(K) BDC = & BDC & + ( BK(L,K) + DK(L,K) + CK(L,K) * CI ) * COSINT(K) DODUM = & DODUM & + ( E0 * CI * AKP(L,K) / SI - CK(L,K) ) * COSINT(K) DLDUM = & DLDUM & + COSINT(K) & * ( ( 1D0 + Q ) * BK(L,K) & - E0 * CI2 * AKP(L,K) / SI & + CK(L,K) * CI & - AK(L,K) & ) END IF 113 CONTINUE 112 CONTINUE C1 Compute long periodic perturbations. SEC = 2D0 + 5D0 * E2 & - ( 32D0 + 112D0 * E2 ) * CI2 & + ( 30D0 + 135D0 * E2 ) * CI4 CEC = E0 * Q2 * ( 1D0 - 16D0 * CI2 + 15D0 * CI4 ) FFY = 1D0 - (15D0/14D0) * F QSA = Q2 * SI * AKPSK SAS = ( SI2 - E2 * CI2 ) * AKPCK / SI DX = CW * QSA & - SW * SAS & - ETA0 * BDC & + 0.5D0 * E0 * NJR * SW * SEC * COSINT(2) & - NJR * CW * CEC * SININT(2) DE = SW * QSA & + CW * SAS & + XI0 * BDC & - 0.5D0 * E0 * NJR * CW * SEC * COSINT(2) & - NJR * SW * CEC * SININT(2) DI = - E0 * CI * AKPSK & - 14D0 * NJR * E2 * SI * CI * FFY * SININT(2) DO = DODUM & - NJR * ( E2 * ( 16D0 - 30D0 * CI2 ) * CI ) * COSINT(2) DL = DLDUM & - 0.5D0 * NJR * ( SEC & + ( 5D0 * E2 - 2D0 ) * CEC / Q / E0 & ) * COSINT(2) C1 Compute secular perturbations. WT = WDOT1 * TSTEP CWT1 = DCOS( WT ) - 1D0 SWT = DSIN( WT ) DXBAR = XI0 * CWT1 - ETA0 * SWT DEBAR = ETA0 * CWT1 + XI0 * SWT DLBAR = ( N +LDOT1 ) * TSTEP C1 ********************************************************************* C1 ********************************************************************* C1 *** *** C1 *** Luni-Solar Gravity Model *** C1 *** *** C1 ********************************************************************* C1 ********************************************************************* C1 If LSFLAG = .FALSE., then ignore luni-solar effects. IF ( .NOT. LSFLAG ) THEN DXMOON = 0D0 DXSUN = 0D0 DEMOON = 0D0 DESUN = 0D0 DIMOON = 0D0 DISUN = 0D0 DOMOON = 0D0 DOSUN = 0D0 DLMOON = 0D0 DLSUN = 0D0 GO TO 901 END IF C1 Determine the classical orbital elements for the Sun and Moon C1 at the given epoch. CALL LUNORB( DATE , LUNSTT ) CALL SUNORB( DATE , SUNSTT ) C1 Compute required functions of elements of Sun and Moon. AM = LUNSTT(1) AM2 = AM * AM AM3 = AM2 * AM AS = SUNSTT(1) AS2 = AS * AS AS3 = AS2 * AS EM = LUNSTT(2) EM2 = EM * EM EM3 = EM2 * EM EM4 = EM3 * EM ES = SUNSTT(2) ES2 = ES * ES ES3 = ES2 * ES ES4 = ES3 * ES CIM = DCOS( LUNSTT(3) * D2R ) CIS = DCOS( SUNSTT(3) * D2R ) SIM = DSIN( LUNSTT(3) * D2R ) SIM2 = SIM * SIM SIS = DSIN( SUNSTT(3) * D2R ) SIS2 = SIS * SIS DIFFOM = O0 - LUNSTT(4) * D2R DIFFOS = O0 - SUNSTT(4) * D2R WM = LUNSTT(5) * D2R WS = SUNSTT(5) * D2R MM = LUNSTT(6) * D2R MS = SUNSTT(6) * D2R C1 Compute the spacecraft's eccentricity function G(lpq). Only compute C1 G for the required range of indices. GSC(2,1,0) = ( 1D0 - E2 ) ** -1.5D0 C1 Compute the spacecraft's eccentricity function dG(lpq)/de. Only C1 compute dG/de for the required range of indices. DGDE(2,1,0) = 3D0 * E0 * ( 1D0 - E2 ) ** -2.5D0 C1 Compute the moon's eccentricity function G(lhj). Only compute C1 G for the required range of indices. GMOON(2,0,-1) = -.5D0 * EM + EM ** 3D0 / 16D0 GMOON(2,0, 0) = 1D0 - 2.5D0 * EM2 + (13D0/16D0) * EM4 GMOON(2,0, 1) = 3.5D0 * EM - (123D0/16D0) * EM3 GMOON(2,1,-1) = 1.5D0 * EM + (27D0/16D0) * EM3 GMOON(2,1, 0) = ( 1D0 - EM2 ) ** -1.5D0 GMOON(2,1, 1) = GMOON(2,1,-1) GMOON(2,2,-1) = GMOON(2,0, 1) GMOON(2,2, 0) = GMOON(2,0, 0) GMOON(2,2, 1) = GMOON(2,0,-1) C1 Compute the Sun's eccentricity function G(lhj). Only compute C1 G for the required range of indices. GSUN(2,0,-1) = -.5D0 * ES + ES ** 3D0 / 16D0 GSUN(2,0, 0) = 1D0 - 2.5D0 * ES2 + (13D0/16D0) * ES4 GSUN(2,0, 1) = 3.5D0 * ES - (123D0/16D0) * ES3 GSUN(2,1,-1) = 1.5D0 * ES + (27D0/16D0) * ES3 GSUN(2,1, 0) = ( 1D0 - ES2 ) ** -1.5D0 GSUN(2,1, 1) = GSUN(2,1,-1) GSUN(2,2,-1) = GSUN(2,0, 1) GSUN(2,2, 0) = GSUN(2,0, 0) GSUN(2,2, 1) = GSUN(2,0,-1) C1 Compute the spacecraft's inclination function F(lmp). Only compute C1 F for the required range of indices. FSC(2,0,1) = .75D0 * SI2 - .5D0 FSC(2,1,1) = -1.50D0 * SICI FSC(2,2,1) = 1.50D0 * SI2 C1 Compute the spacecraft's inclination function dF(lmp)/di. Only C1 compute dF/di for the required range of indices. DFDI(2,0,1) = 1.5D0 * SICI DFDI(2,1,1) = -1.5D0 * ( 2D0 * CI2 - 1.D0 ) DFDI(2,2,1) = 3.0D0 * SICI C1 Compute the moon's inclination function F(lmh). Only compute C1 F for the required range of indices. FMOON(2,0,0) = -.375D0 * SIM2 FMOON(2,0,1) = .750D0 * SIM2 - .5D0 FMOON(2,0,2) = -.375D0 * SIM2 FMOON(2,1,0) = .750D0 * SIM * ( 1D0 + CIM ) FMOON(2,1,1) = -1.500D0 * SIM * CIM FMOON(2,1,2) = -.750D0 * SIM * ( 1D0 - CIM ) FMOON(2,2,0) = .750D0 * ( 1D0 + CIM ) ** 2D0 FMOON(2,2,1) = 1.500D0 * SIM2 FMOON(2,2,2) = .750D0 * ( 1D0 - CIM ) ** 2D0 C1 Compute the Sun's inclination function F(lmh). Only compute C1 F for the required range of indices. FSUN(2,0,0) = -.375D0 * SIS2 FSUN(2,0,1) = .750D0 * SIS2 - .5D0 FSUN(2,0,2) = -.375D0 * SIS2 FSUN(2,1,0) = .750D0 * SIS * ( 1D0 + CIS ) FSUN(2,1,1) = -1.500D0 * SIS * CIS FSUN(2,1,2) = -.750D0 * SIS * ( 1D0 - CIS ) FSUN(2,2,0) = .750D0 * ( 1D0 + CIS ) ** 2D0 FSUN(2,2,1) = 1.500D0 * SIS2 FSUN(2,2,2) = .750D0 * ( 1D0 - CIS ) ** 2D0 C1 Besides the eccentricity and inclination functions, we might as well C1 compute other products outside of the nested index loops, even though C1 this may not save any execution time. C2 Compute products from the spacecraft's orbital elements. QE = Q / E0 QSI = Q * SI CIQSI = CI / QSI NA2 = N * A2 NA2QSI = NA2 * QSI C2 Compute the power ratio for the expected range of the L index of the C2 spacecraft's semi-major axis to that of the sun and moon. ALMOON(2) = A2 / AM3 ALSUN (2) = A2 / AS3 C2 Compute the partial product D*(L-M)!/(L+M)!, where D=1 when M=0 and C2 D=2 when not M=0. DLM(2,0) = 1D0 DLM(2,1) = 1D0 / 3D0 DLM(2,2) = 1D0 / 12D0 C1 Compute the change in the spacecraft's non-singular elements due C1 to both the Sun and Moon. For execution speed, the nested loops C1 over L, M, P, H, and J and intermingled for both bodies. C1 C1 The non-singular elements computed are xi, eta, inclination, longitude C1 of the ascending node, and latitude (w + M). C1 C1 Luni-solar effects do not change the semi-major axis of the spacecraft. C1 C1 Note that loop ranges for indices M and H only go up to L, not C1 MHI or HHI. However, indices P and J range explicitly from their C1 low to high values. C1 C1 If LSRGP is modified, then the you should make sure that the loop C1 range for each index reflects the new index range. C2 Set L index accumulator values to zero for both Sun and Moon. DEML = 0D0 DIML = 0D0 DLML = 0D0 DOML = 0D0 DXML = 0D0 DESL = 0D0 DISL = 0D0 DLSL = 0D0 DOSL = 0D0 DXSL = 0D0 C2 Begin the L index loop. DO 201 L = LLO , LHI C3 Set M index accumulator values to zero for both Sun and Moon. DEMM = 0D0 DIMM = 0D0 DLMM = 0D0 DOMM = 0D0 DXMM = 0D0 DESM = 0D0 DISM = 0D0 DLSM = 0D0 DOSM = 0D0 DXSM = 0D0 C3 Begin the M index loop. DO 202 M = MLO , L C4 Convert M from integer to double precision. MDP = DFLOAT( M ) C4 Set P index accumulator values to zero for both Sun and Moon. DEMP = 0D0 DIMP = 0D0 DLMP = 0D0 DOMP = 0D0 DXMP = 0D0 DESP = 0D0 DISP = 0D0 DLSP = 0D0 DOSP = 0D0 DXSP = 0D0 C4 Begin the P index loop. DO 203 P = PLO , PHI C5 Set H index accumulator values to zero for both Sun and Moon. DEMH = 0D0 DIMH = 0D0 DLMH = 0D0 DOMH = 0D0 DXMH = 0D0 DESH = 0D0 DISH = 0D0 DLSH = 0D0 DOSH = 0D0 DXSH = 0D0 C5 Begin H index loop. DO 204 H = HLO , L C6 Compute L - 2H for later use in computing integrals. L2H = DFLOAT( L - 2 * H ) C6 Set J index accumulator values to zero for both Sun and Moon. DEMJ = 0D0 DIMJ = 0D0 DLMJ = 0D0 DOMJ = 0D0 DXMJ = 0D0 DESJ = 0D0 DISJ = 0D0 DLSJ = 0D0 DOSJ = 0D0 DXSJ = 0D0 C6 Begin J index loop. (Using JJ, since J is for J2, J3, ...) DO 205 JJ = JLO , JHI C7 Compute the integrals of sin(y0) and cos(y0) for C7 both the Sun and the Moon. L2HJ = L2H + DFLOAT( JJ ) Y0M = MDP * DIFFOM - L2H * WM - L2HJ * MM Y0S = MDP * DIFFOS - L2H * WS - L2HJ * MS CY0M = DCOS( Y0M ) SY0M = DSIN( Y0M ) CY0S = DCOS( Y0S ) SY0S = DSIN( Y0S ) XM = TSTEP * ( MDP * ODOT1 - L2HJ * NMOON ) XM2 = XM * XM XM4 = XM2 * XM2 XM6 = XM4 * XM2 XM8 = XM6 * XM2 XM10 = XM8 * XM2 XM12 = XM10 * XM2 XM14 = XM12 * XM2 XS = TSTEP * ( MDP * ODOT1 - L2HJ * NSUN ) XS2 = XS * XS XS4 = XS2 * XS2 XS6 = XS4 * XS2 XS8 = XS6 * XS2 XS10 = XS8 * XS2 XS12 = XS10 * XS2 XS14 = XS12 * XS2 F1M = 1D0 & - XM2 / 6D0 & + XM4 / 120D0 & - XM6 / 5040D0 & + XM8 / 362880D0 & - XM10 / 39916800D0 & + XM12 / 6227020800D0 & - XM14 / 1307674368000D0 F2M = .5D0 & - XM2 / 24D0 & + XM4 / 720D0 & - XM6 / 40320D0 & + XM8 / 3628800D0 & - XM10 / 479001600D0 & + XM12 / 87178291200D0 & - XM14 / 20922789888000D0 F1S = 1D0 & - XS2 / 6D0 & + XS4 / 120D0 & - XS6 / 5040D0 & + XS8 / 362880D0 & - XS10 / 39916800D0 & + XS12 / 6227020800D0 & - XS14 / 1307674368000D0 F2S = .5D0 & - XS2 / 24D0 & + XS4 / 720D0 & - XS6 / 40320D0 & + XS8 / 3628800D0 & - XS10 / 479001600D0 & + XS12 / 87178291200D0 & - XS14 / 20922789888000D0 ICM = TSTEP * ( F1M * CY0M - XM * F2M * SY0M ) ISM = TSTEP * ( XM * F2M * CY0M + F1M * SY0M ) ICS = TSTEP * ( F1S * CY0S - XS * F2S * SY0S ) ISS = TSTEP * ( XS * F2S * CY0S + F1S * SY0S ) C7 Compute common terms for J index. GMNICM = GMOON(L,H,JJ) * ICM GSNICS = GSUN (L,H,JJ) * ICS C7 Update xi element wrt Moon for J index. DXMJ = DXMJ + GMNICM C7 Update xi element wrt Sun for J index. DXSJ = DXSJ + GSNICS C7 Update eta element wrt Moon for J index. DEMJ = DEMJ + GMNICM C7 Update eta element wrt Sun for J index. DESJ = DESJ + GSNICS C7 Update inclination element wrt Moon for J index. DIMJ = DIMJ + GMOON(L,H,JJ) * ISM C7 Update inclination element wrt Sun for J index. DISJ = DISJ + GSUN(L,H,JJ) * ISS C7 Update LAN element wrt Moon for J index. DOMJ = DOMJ + GMNICM C7 Update LAN element wrt Sun for J index. DOSJ = DOSJ + GSNICS C7 Update argument of latitude element wrt Moon for J index. DLMJ = DLMJ + GMNICM C7 Update argument of latitude element wrt Sun for J index. DLSJ = DLSJ + GSNICS C6 End J index loop. 205 CONTINUE C6 Update xi element wrt Moon for H index. DXMH = DXMH + DXMJ * FMOON(L,M,H) C6 Update xi element wrt Sun for H index. DXSH = DXSH + DXSJ * FSUN(L,M,H) C6 Update eta element wrt Moon for H index. DEMH = DEMH + DEMJ * FMOON(L,M,H) C6 Update eta element wrt Sun for H index. DESH = DESH + DESJ * FSUN(L,M,H) C6 Update inclination element wrt Moon for H index. DIMH = DIMH + DIMJ * FMOON(L,M,H) C6 Update inclination element wrt Sun for H index. DISH = DISH + DISJ * FSUN(L,M,H) C6 Update LAN element wrt Moon for H index. DOMH = DOMH + DOMJ * FMOON(L,M,H) C6 Update LAN element wrt Sun for H index. DOSH = DOSH + DOSJ * FSUN(L,M,H) C6 Update argument of latitude element wrt Moon for H index. DLMH = DLMH + DLMJ * FMOON(L,M,H) C6 Update argument of latitude element wrt Sun for H index. DLSH = DLSH + DLSJ * FSUN(L,M,H) C5 End H index loop. 204 CONTINUE C5 Compute common term for xi and eta. DXP = ( CIQSI * DFDI(L,M,P) * GSC(L,P,0) & - QE * FSC(L,M,P) * DGDE(L,P,0) & ) C5 Update xi element wrt Moon for P index. DXMP = DXMP + DXMH * DXP C5 Update xi element wrt SUN for P index. DXSP = DXSP + DXSH * DXP C5 Update eta element wrt Moon for P index. DEMP = DEMP - DEMH * DXP C5 Update eta element wrt SUN for P index. DESP = DESP - DESH * DXP C5 Compute common term for inclination and LAN. GSCFSC = GSC(L,P,0) * FSC (L,M,P) GSCDF = GSC(L,P,0) * DFDI(L,M,P) C5 Update inclination element wrt Moon for P index. DIMP = DIMP + DIMH * GSCFSC C5 Update inclination element wrt SUN for P index. DISP = DISP + DISH * GSCFSC C5 Update LAN element wrt Moon for P index. DOMP = DOMP + DOMH * GSCDF C5 Update LAN element wrt SUN for P index. DOSP = DOSP + DOSH * GSCDF C5 Compute common term for argument of latitude. DLP = GSC(L,P,0) & * ( CIQSI * DFDI(L,M,P) & + DFLOAT( 2 * L ) * FSC(L,M,P) & ) C5 Update argument of latitude element wrt Moon for P index. DLMP = DLMP + DLMH * DLP C5 Update argument of latitude element wrt SUN for P index. DLSP = DLSP + DLSH * DLP C4 End the P index loop. 203 CONTINUE C4 Update xi element wrt Moon for M index. DXMM = DXMM + DXMP * DLM(L,M) C4 Update xi element wrt Sun for M index. DXSM = DXSM + DXSP * DLM(L,M) C4 Update eta element wrt Moon for M index. DEMM = DEMM + DEMP * DLM(L,M) C4 Update eta element wrt Sun for M index. DESM = DESM + DESP * DLM(L,M) C4 Compute common term for inclination. MDLM = MDP * DLM(L,M) C4 Update inclination element wrt Moon for M index. DIMM = DIMM + DIMP * MDLM C4 Update inclination element wrt Sun for M index. DISM = DISM + DISP * MDLM C4 Update LAN element wrt Moon for M index. DOMM = DOMM + DOMP * DLM(L,M) C4 Update LAN element wrt Sun for M index. DOSM = DOSM + DOSP * DLM(L,M) C4 Update argument of latitude element wrt Moon for M index. DLMM = DLMM + DLMP * DLM(L,M) C4 Update argument of latitude element wrt Sun for M index. DLSM = DLSM + DLSP * DLM(L,M) C3 End the M index loop. 202 CONTINUE C3 Update xi element wrt Moon for L index. DXML = DXML + DXMM * ALMOON(L) C3 Update xi element wrt Sun for L index. DXSL = DXSL + DXSM * ALSUN(L) C3 Update eta element wrt Moon for L index. DEML = DEML + DEMM * ALMOON(L) C3 Update eta element wrt Sun for L index. DESL = DESL + DESM * ALSUN(L) C3 Update inclination element wrt Moon for L index. DIML = DIML + DIMM * ALMOON(L) C3 Update inclination element wrt Sun for L index. DISL = DISL + DISM * ALSUN(L) C3 Update LAN element wrt Moon for L index. DOML = DOML + DOMM * ALMOON(L) C3 Update LAN element wrt Sun for L index. DOSL = DOSL + DOSM * ALSUN(L) C3 Update argument of latitude element wrt Moon for L index. DLML = DLML + DLMM * ALMOON(L) C3 Update argument of latitude element wrt Sun for L index. DLSL = DLSL + DLSM * ALSUN(L) C2 End the L index loop. 201 CONTINUE C2 Compute common terms for Sun and Moon. GMMNA2 = GMMOON / NA2 GMSNA2 = GMSUN / NA2 GMMNQS = GMMOON / NA2QSI GMSNQS = GMSUN / NA2QSI C2 Complete change to xi due to Moon. DXMOON = DXML * ETA0 * GMMNA2 C2 Complete change to xi due to Sun. DXSUN = DXSL * ETA0 * GMSNA2 C2 Complete change to eta due to Moon. DEMOON = DEML * XI0 * GMMNA2 C2 Complete change to eta due to Sun. DESUN = DESL * XI0 * GMSNA2 C2 Complete change to inclination due to Moon. DIMOON = DIML * GMMNQS C2 Complete change to inclination due to Sun. DISUN = DISL * GMSNQS C2 Complete change to LAN due to Moon. DOMOON = DOML * GMMNQS C2 Complete change to LAN due to Sun. DOSUN = DOSL * GMSNQS C2 Complete change to argument of latitude due to Moon. DLMOON = - DLML * GMMNA2 C2 Complete change to argument of latitude due to Sun. DLSUN = - DLSL * GMSNA2 C1 ********************************************************************* C1 ********************************************************************* C1 *** *** C1 *** Add Effects of Both Earth and Luni-Solar Gravity *** C1 *** *** C1 ********************************************************************* C1 ********************************************************************* C1 Target for GO TO if LSFLAG is false. 901 CONTINUE C1 Propagate non-singular elements using both Earth and Luni-Solar C1 models (luni-solar effects may be zero if LSFLAG is false). A1 = A0 XI1 = XI0 + DX + DXBAR + DXMOON + DXSUN ETA1 = ETA0 + DE + DEBAR + DEMOON + DESUN I1 = I0 + DI + DIMOON + DISUN O1 = O0 + DO + ODOT1 * TSTEP + DOMOON + DOSUN L1 = L0 + DL + DLBAR + DLMOON + DLSUN C1 Convert non-singular elements to Keplerian classical orbital elements. ORB1(1) = A1 ORB1(2) = DSQRT( XI1**2 + ETA1**2 ) ORB1(3) = I1 / D2R ORB1(4) = DMOD( O1 / D2R , 360D0 ) IF ( ORB1(4) .LT. 0D0 ) ORB1(4) = ORB1(4) + 360D0 ORB1(5) = ATAN2( ETA1 , XI1 ) / D2R IF ( ORB1(5) .LT. 0D0 ) ORB1(5) = ORB1(5) + 360D0 ORB1(6) = DMOD( L1 / D2R - ORB1(5) , 360D0 ) IF ( ORB1(6) .LT. 0D0 ) ORB1(6) = ORB1(6) + 360D0 C1 Before we go, let me post mean motion on the LSRGP Bulletin Board. CALL LSRGPB( 'PUT' , 'LSRGP_N' , N + N1 ) C1 End of LSRGP. RETURN END