C$Procedure MNODES C SUBROUTINE MNODES ( ORBIN , DATIN , & M , jearth , & LTOP , LSFLAG , & DRAG , & srpflag, cr, mass,dsmadt, tsmaswitch, & sigma_dsma, HILO, & ORBOUT , DATOUT , SUNalpha, & BETAPrime ) C C C******************************************************************************* C C Copyright (C) 1993, California Institute of Technology. U.S. C Government Sponsorhip under NASA Contract NAS7-918 is C acknowledged. C C******************************************************************************* C C$ Log C C Date Name Description C ----------------------------------------------------------------------------- C 11-JUL-1990 Eric Cannell creation of MNODES C 13-SEp-1991 Bruce Shapiro add JEARTH as parameter C C$ Purpose C C MNODES determines the classical elements and epoch of the Mth ascending C node. C C$ Input_Arguments C C Name Type Dim Units Description C ----------------------------------------------------------------------------- C ORBIN DP 6 km,deg classical elements (a,e,i,LAN,w,M) C DATIN C*(*) 1 - epoch of ORBIN, TIMETRANS format C 'dd-mmm-yyyy hh:mm:ss.ffff' C M I 1 - index of desired ascending node. 0 < M. C JEARTH DP (2:29) - earth zonals C LTOP I 1 - LTOP passed to LSRGP C LSFLAG L 1 - if true, luni-solar effects are turned C on in LSRGP C DRAG L 1 - if true, atmoshperic drag is turned on C dsmadt r*8 m/day C tsmaswitch r*8 C sigma_dsma C HILO c*4 1 - which drag to use, HI, LO, TRUE C C$ Output_Arguments C C Name Type Dim Units Description C ----------------------------------------------------------------------------- C ORBOUT DP 6 km,deg classical elements (a,e,i,LAN,w,M) at C the Mth ascending node C DATOUT C*(*) 1 - epoch of Mth ascending node, TIMETRANS C format 'dd-mmm-yyyy hh:mm:ss.ffff' C BETAPrime DP 1 deg beta prime angle at START OF RUN C C$ Restrictions C C 1] MNODES is strictly for Earth orbiters only, since the LSRGP library C is used to propagate the orbit. C C 2] ORBIN must not have a inclination of zero. MNODES does NOT check C for i=0, since execution time is a valuable commodity. C C 3] Index of desired ascending node (M) must be greater than zero. MNODES C does NOT check for M>0, since execution time is a valuable commodity. C C 4] If the initial orbit is at an ascending node (u0=0), then that node C is the zeroth node, not the first node. C C 5] Often a programmer will use MNODES to iterate on an orbit, feeding C the orbit at the ascending node back into MNODES. Unfortunately, C the orbit returned for the ascending node MIGHT have an exceedingly C small -Z component. In such a case, simply feeding the orbit back to C MNODES would have caused either the same ascending node to be counted C twice OR crash LSRGP when the propagation time is on the order of C 10**-15 (the time needed to go from, say, Z=-10**-13 to Z=0). C C =====> Because GTARG runs out a ground track by propagating an orbit C M ascending nodes from the current ascending node (M>0), C I have arbitrarily decided that if the initial argument of C latitude (U0) is in the range of [359.999,360.000] then C MNODES will set U0 to zero. In which case, Restriction #4 C is invoked. C C The only possible loss is that MNODES may skip the very first C ascending node of a ground track if the initial orbit provided C to GTARG happens to have an exceedingly small Z-component. Such C a loss does not affect the results of GTARG. C C$ Library_Links C C Entry Point Name Location C ----------------------------------------------------------------------------- C CH2SEC TIMETRANS C DVMOVE TPXUTIL C ORB2U TPXORB C PROP GTARG C SEC2CH TIMETRANS C C common area for physical constants C double precision earth_rad ! in kilomters double precision earth_freq ! radians / second double precision earth_rate ! kilomters / day double precision GMrth ! km**3/sec double precision mu_moon ! km**3/sec double precision mu_sun ! km**3/sec double precision sid_day ! seconds double precision deg_to_km ! kilometers/deg common / physical_constants / & earth_rad, earth_freq, earth_rate, GMrth, & mu_moon, mu_sun, sid_day, deg_to_km C$ Parameters C DOUBLE PRECISION EPS PARAMETER ( EPS = 1.0D-5 ) C DOUBLE PRECISION GMRTH C PARAMETER ( GMRTH = 398600.44807345D0 ) INTEGER KMAX PARAMETER ( KMAX = 20 ) DOUBLE PRECISION PI PARAMETER ( PI = 3. 14159 26535 89793 23846 D0 ) DOUBLE PRECISION D2R PARAMETER ( D2R = PI / 180D0 ) DOUBLE PRECISION TWOPI PARAMETER ( TWOPI = 2D0 * PI ) C C$ Declarations_of_Input_and_Output_Arguments C double precision betap, betaprime, alpha, sunalpha CHARACTER*(*) DATIN CHARACTER*(*) DATOUT CHARACTER*4 HILO double precision jearth ( 2:29 ) LOGICAL DRAG, srpflag double precision cr, mass double precision dsmadt(2), tsmaswitch double precision sigma_dsma(2) LOGICAL LSFLAG INTEGER LTOP INTEGER M DOUBLE PRECISION ORBIN ( 6 ) DOUBLE PRECISION ORBOUT ( 6 ) C C$ Declarations_of_Local_Variables C CHARACTER*25 DATSHY LOGICAL ERROR DOUBLE PRECISION JSEC INTEGER KFLAG DOUBLE PRECISION O1 ( 6 ) DOUBLE PRECISION O2 ( 6 ) DOUBLE PRECISION ORBSHY ( 6 ) DOUBLE PRECISION P DOUBLE PRECISION SU1 DOUBLE PRECISION SU2 DOUBLE PRECISION T1 DOUBLE PRECISION T2 DOUBLE PRECISION T3 DOUBLE PRECISION T4MORE DOUBLE PRECISION TSHY DOUBLE PRECISION U0 DOUBLE PRECISION U1 DOUBLE PRECISION U2 DOUBLE PRECISION U2GO C C$ External_Statements C DOUBLE PRECISION CH2SEC EXTERNAL CH2SEC CHARACTER*25 SEC2CH EXTERNAL SEC2CH C C$ Method C C In general- C C 1- Propagate orbit forward in one step to a point about 2 degrees C shy of the Mth ascending node. C C 2- Propagate about 4 degrees beyond that point to a second point C in the +z hemisphere. C C 3- Given two orbits that bound the Mth ascending node, use the C secant method to search for the orbit and epoch where the C sin(u) is zero. This is equivalent to finding the root C of: C C sin( u(t) ) = 0. C C Since the function sin(u(t)) is a well-behaved function of t, the C secant method will converge quickly to the desired root. C C 4- Note that sin(u(t)) is zero at both the ascending and descending C nodes. However, the initial step taken by MNODES will place the C orbit just shy of the ascending node. C C Specifically- C-& C1 Estimate the orbital period (in seconds) during this time to be C1 roughly a constant 2pi * sqrt( a**3 / gm ). P = TWOPI * DSQRT( ORBIN(1)**3 / GMRTH ) C1 Determine the initial U0. CALL ORB2U( ORBIN , U0 ) C1 Invoke Restriction #5, if necessary. IF ( U0 .GT. 359.999 ) U0 = 0D0 C1 Compute the remaining U-to-go to the next ascending node. This is C1 the first ascending node. U2GO = 360D0 - U0 C1 Propagate initial orbit up to about 2 degrees shy of Mth ascending C1 node. Note that this places the orbit in the -z hemisphere near an C1 ascending node, not a descending node. TSHY = ( DFLOAT( M - 1 ) & + ( U2GO - 2D0 ) & / 360D0 & ) & * P C write(8,*) 'MNODES1 hilo = ',hilo CALL PROP( ORBIN , DATIN , TSHY , JEARTH, LTOP , & LSFLAG , DRAG , & srpflag, cr, mass,dsmadt, tsmaswitch, & sigma_dsma, HILO, O1, sunalpha, BETAPRIME ) C write(8,*)'MNODES BETAP=',betap C1 Compute the epoch at this point. JSEC = CH2SEC( DATIN , .TRUE. , ERROR ) + TSHY DATSHY = SEC2CH( JSEC ) C1 Propagate orbit about 4 degrees more, placing orbit in +z hemisphere. T4MORE = ( 4D0 / 360D0 ) * P C write(8,*) 'MNODES2 hilo = ',hilo CALL PROP( O1 , DATSHY , T4MORE , jearth , & LTOP , LSFLAG , DRAG , & srpflag, cr, mass, dsmadt, tsmaswitch, & sigma_dsma, HILO, O2, alpha, BETAP ) C write(8,*)'MNODES BETAP=',betap C1 Setup for secant search for ascending node. Let the -z orbit be C1 designated (O1,T1,SU1) and the +z orbit be designated (O2,T2,SU2), C1 where SUn = sin( Un ). The node is somewhere in between orbit O1 and C1 orbit O2. T1 = 0D0 T2 = T4MORE CALL ORB2U( O1 , U1 ) CALL ORB2U( O2 , U2 ) SU1 = DSIN( U1 * D2R ) SU2 = DSIN( U2 * D2R ) C1 To avoid propagating an orbit backwards in time (although LSRGP should C1 be able to do it), save O1 and always propagate from there. CALL DVMOVE( 6 , O1 , ORBSHY ) C1 Iterate with the secant method until the delta time satisfies C1 the stopping criteria, until the change in sin( U ) is zero, or C1 until the maximum number of iterations has been reached. One hopes C1 that if delta(sin(U)) is zero, then U is zero. KFLAG = 0 901 CONTINUE C2 Increment the iteration counter. KFLAG = KFLAG + 1 C2 Determine the next time iterate. T3 = T1 - SU1 * ( T2 - T1 ) / ( SU2 - SU1 ) C2 Transfer the current orbit's time and sin(U). Do not bother to C2 save the orbital elements, since MNODES always propagate from ORBSHY. T1 = T2 SU1 = SU2 C2 Generate the next orbit. T2 = T3 C write(8,*) 'MNODES3 hilo = ',hilo CALL PROP( ORBSHY , DATSHY , T2 , jearth, & LTOP , LSFLAG , DRAG , & srpflag, cr, mass, dsmadt, tsmaswitch, & sigma_dsma, HILO, O2, alpha, BETAP ) c write(8,*)'MNODES BETAP=',betap CALL ORB2U( O2 , U2 ) SU2 = DSIN( U2 * D2R ) C2 If a stopping criteria was reached, exit the secant method. IF ( DABS( T2 - T1 ) .LE. EPS & .OR. SU1 .EQ. SU2 & .OR. KFLAG .EQ. KMAX ) GOTO 902 C2 If not, begin another iteration. GOTO 901 902 CONTINUE C1 Once a stopping criteria was reached, return the last iterate C1 as the final orbit. CALL DVMOVE( 6 , O2 , ORBOUT ) C1 Update the nodal epoch. JSEC = CH2SEC( DATSHY , .TRUE. , ERROR ) + T2 DATOUT = SEC2CH( JSEC ) C1 End of MNODES. RETURN END