C$Procedure PROP C SUBROUTINE PROP ( ORBIN , DATIN , & DT , jearth , & LTOP , LSFLAG , DRAG , & dsmadt , tsmaswitch, & sigma_dsma, HILO, & ORBOUT , BETAP ) 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 PROP C 13-Sep-1991 Bruce Shapiro jearth as a parameter C Dec-1992 additional along track force as table look-up C into array of da/dt values C C$ Purpose C C PROP propagates an orbit using the LSRGP Library and adjusts the C semi-major axis to reflect atmospheric drag. 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 DT DP 1 sec length of propagation time C jearth dp (2:29) - earth gravity field zonal coefficients 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 1 meters/day extra da/dt C HILO C*4 1 - 'HI' - use high density C 'LO' - use low density C 'TRUE' - use "true" density 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) after C time DT C C$ Restrictions C C 1] PROP 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. PROP does NOT check C for i=0, since execution time is a valuable commodity. C C$ Library_Links C C Entry Point Name Location C ----------------------------------------------------------------------------- C B2DRG TPXLIB C CRMGET GTARG C LSRGP LSRGP C LSRGPB LSRGP C ORBBP TPXORB C SUNORB TPXORB C double precision earth_rad ! in kilomters double precision WE ! radians / second double precision earth_rate ! meters / 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, WE, earth_rate, GMrth, & mu_moon, mu_sun, sid_day, deg_to_km C$ Parameters C C DOUBLE PRECISION GMRTH C PARAMETER ( GMRTH = 398600.44807345D0 ) logical trace, dsmatrace parameter ( trace = .false. ) parameter (dsmatrace = .false. ) DOUBLE PRECISION PI PARAMETER ( PI = 3. 14159 26535 89793 23846 D0 ) DOUBLE PRECISION D2R PARAMETER ( D2R = PI / 180D0 ) C DOUBLE PRECISION SPD C PARAMETER ( SPD = 86400D0 ) C DOUBLE PRECISION WE C PARAMETER ( WE = D2R * 360.9856473458400D0 / SPD ) C C$ Declarations_of_Input_and_Output_Arguments C CHARACTER*(*) DATIN LOGICAL DRAG DOUBLE PRECISION dsmadt(2), dsmaswitch, DT, tnow double precision sigma_dsma(2) CHARACTER * 4 HILO double precision jearth ( 2:29 ) LOGICAL LSFLAG INTEGER LTOP DOUBLE PRECISION ORBIN ( 6 ) DOUBLE PRECISION ORBOUT ( 6 ) common /boost/ ndsmadt_data, dsmadt_data, dsmadt_epoch, & dsmadt_data_sigma, plot_boost, & dsmadt_dates, xdsmadt_dates logical plot_boost double precision dsmadt_data(1000), xdsmadt_dates(1000) double precision dsmadt_data_sigma character*25 dsmadt_epoch, dsmadt_dates(1000) integer ndsmadt_data c C C$ Declarations_of_Local_Variables C DOUBLE PRECISION AREA DOUBLE PRECISION BETAP DOUBLE PRECISION CI DOUBLE PRECISION CRM, crmhi, crmlo DOUBLE PRECISION DADT, dadt1, dadt2, dadtnow, xt1,xt2 logical error DOUBLE PRECISION N DOUBLE PRECISION SUNSTT ( 6 ) DOUBLE PRECISION V2 double precision error_sign double precision tdsma C C$ External_Statements C DOUBLE PRECISION VMAREA EXTERNAL VMAREA DOUBLE PRECISION CH2SEC EXTERNAL CH2SEC integer finddate external finddate DOUBLE PRECISION ORBBP EXTERNAL ORBBP CHARACTER*25 SEC2CH EXTERNAL SEC2CH C C$ Method C-& C1 Propagate the orbit using LSRGP. This does not include drag effects. if (trace) then write(50,*) 'PROP (before LSRGP) ORBIN,DATIN,dt:',orbin,datin,dt end if CALL LSRGPj( ORBIN , DATIN , DT , JEARTH, LTOP , LSFLAG , ORBOUT ) if (trace) then write(50,*) 'PROP (after LSRGP) ORBOUT:',orbout end if C1 If desired, modify the semi-major axis to reflect the effects of C1 atmospheric drag. IF ( DRAG ) THEN C2 Compute square of satellite velocity. V2 = GMRTH / ORBOUT( 1 ) C2 Use ORBBP to determine the beta-prime angle of the Sun. CALL SUNORB( DATIN , SUNSTT ) BETAP = ORBBP( ORBIN , SUNSTT ) C2 Use B2DRG to determine the projected area of the TOPEX satellite C2 due to atmospheric drag in km**2. C AREA = B2DRG( BETAP ) * 1.0D-6 AREA = VMAREA( BETAP, 'DRAG' ) * 1.0D-6 C2 Compute cosine of inclination of satellite. CI = DCOS( D2R * ORBIN( 3 ) ) C2 Get the mean motion off of the LSRGP Bulletin Board. C C need to test at some point for invalid HILO value C also: may want to keep three separate bulletin boards ??? C for now, just use one and hope things don't get confused C if (hilo .eq. 'TRUE') then CALL LSRGPB( 'GET' , 'LSRGP_N' , N ) else if (hilo .eq. 'HI') then CALL LSRGPB( 'GET' , 'LSRGP_N' , N ) else if (hilo .eq. 'LO') then CALL LSRGPB( 'GET' , 'LSRGP_N' , N ) else write(6,*) 'PROP: Unknown HILO = '//hilo stop 'ERROR EXIT' end if C2 Determine the product of Cd * density / mass which is computed C2 at the time the PFLUX file is read. CALL CRMGET( DATIN , CRM, crmhi, crmlo ) C C use the 'TRUE' CRM unless otherwise requested C if (hilo .eq. 'HI') then CRM = CRMHI else if (hilo .eq. 'LO') then CRM = CRMLO end if if (trace) then write (50,*) 'PROP CRM:',crm end if C2 Compute da/dt in km/sec. DADT = - CRM & * AREA & * V2 & * ( 1D0 - WE * CI / N ) ** 2D0 & / N C2 Adjust semi-major axis. ORBOUT( 1 ) = ORBOUT( 1 ) + DADT * DT END IF C C C add extra (d/dt)sma C error_sign = 0.0d0 if (hilo .eq. 'HI') then error_sign = -1.0d0 else if (hilo .eq. 'LO') then error_sign = 1.0d0 end if tnow = ch2sec ( datin , .true., error) C C ***** C apply the selected boost model C if (ndsmadt_data .gt. 0 ) then if (dsmadt_epoch .eq. ' ') then C C use DSMADT_DATES and use the most recent value C supplied for da/dt C idate = finddate ( datin, .false., dsmadt_dates, & xdsmadt_dates, ndsmadt_data ) C C interpolate C if ( (idate .gt. 0) .and. & (idate .le. ndsmadt_data) ) then xt1 = xdsmadt_dates(idate) xt2 = xdsmadt_dates(idate+1) dadt1 = dsmadt_data(idate) dadt2 = dsmadt_data(idate+1) dadtnow = dadt1 + (tnow-xt1)*(dadt2-dadt1)/(xt2-xt1) C write(8,*) datin, ' ',dsmadt_dates(idate), ' ', C & dsmadt_dates(idate+1), ' ', C & dadt1, ' ',dadt2, ' ',dadtnow orbout(1) = orbout(1) + & dadtnow * 0.001d0 * (dt/86400.0d0) else write (8,*) 'Current Epoch = ', datin, & ' is out of range = ', idate, & ' of data, range = 1 to ', ndsmadt_data end if C***** else C C use DSMADT_EPOCH and assume data are C equally spaced at one point per orbit C tdsma = ch2sec ( dsmadt_epoch, .true., error) itnow = int ( (tnow - tdsma)/86400.0d0 ) + 1 if ( (itnow .ge. 1) .or. (itnow.le.ndsmadt_data)) then if (dsmatrace) * write(8,400) datin, itnow,dsmadt_data(itnow), * dsmadt_data_sigma 400 format(1x, a, 1x,i10, f10.4,1x, f10.4) orbout(1) = orbout(1) + & dsmadt_data(itnow) * 0.001d0 * (dt/86400.0d0) else write(8,*) 'Warning: tnow = ',datin, * ' is out of the range of dsmadt_data, i=',itnow end if end if else C C if this model is used, error model may be incorrect; should C use array of dsmadt_data instead of a pair of points C if ( tnow .le. tsmaswitch ) then if ( ( dsmadt(1) .ne. 0.0d0) .or. & (sigma_dsma(1) .ne. 0.0d0) ) then orbout(1) = orbout(1) + & ( dsmadt(1) + error_sign * sigma_dsma(1) ) & * 0.001d0 * (dt/86400.0d0) if ( dsmatrace) write(8,*) datin, ' ',dsmadt(1) end if else if ( (dsmadt(2) .ne. 0.0d0) .or. & (sigma_dsma(2) .ne. 0.0d0 ) ) then orbout(1) = orbout(1) + & ( dsmadt(2) + error_sign * sigma_dsma(2) ) & * 0.001d0 * (dt/86400.0d0) if (dsmatrace) write(8,*) datin, ' ',dsmadt(2) end if end if C1 End of PROP. RETURN END