Driver and Function Subroutines

Document Pages

Page 72

Page 73

Page 74

Page 75

Page 76

Page 77

Page 78

Transcription

The DRIVER controls the input and output function and the selection of the model. The input consists of a program card which specifies the model to be used and the output times and either a G-card or T-card element set.

The DRIVER reads and converts the input elements to units of radians and minutes. These are communicated to the prediction model through the COMMON E1. Values of the physical and mathematical constants are set and communicated through the COMMONs C1 and C2, respectively.

The program card indicates the mathematical model to be used and the start and stop time of prediction as well as the increment of time for output. These times are in minutes since epoch.

In the interest of efficiency the DRIVER sets a flag (IFLAG) the first time the model is called. This flag tells the model to calculate all initialized (time independent) quantities. After initialization, the model subroutine turns off the flag so that all subsequent calls only access the time dependent part of the model. This mode continues until another input case is encountered.

The DRIVER takes the output from the mathematical model (communicated through the COMMON E1) and converts it to units of kilometers and seconds for printout.

The function subroutine ACTAN is passed the values of sine and cosine in that order and it returns the angle in radians within the range of 0 to $2\pi$ . The function subroutine FMOD2P is passed an angle in radians and returns the angle in radians within the range of 0 to $2\pi$ . The function subroutine THETAG is passed the epoch time exactly as it appears on the input element cards.¹ The routine converts this time to days since 1950 Jan 0.0 UTC, stores this in the COMMON E1, and returns the right ascension of Greenwich at epoch (in radians).

FORTRAN IV computer code listings of the routines DRIVER, ACTAN, FMOD2P, and THETAG are given below.

¹ If only one year digit is given (as on standard G-cards) the program assumes the 80 decade. This may be overridden by putting a 2 digit year in columns 30-31 of the first G-card.

FORTRAN Implementation

DRIVER

*     DRIVER                                             3 NOV 80
*
*     WGS-72 PHYSICAL AND GEOPOTENTIAL CONSTANTS
*         CK2= .5*J2*AE**2        CK4=-.375*J4*AE**4
*
      DOUBLE PRECISION EPOCH,DS50
      COMMON/E1/XMO,XNODEO,OMEGAO,EO,XINCL,XNO,XNDT2O,XNDD6O,BSTAR,
     1          X,Y,Z,XDOT,YDOT,ZDOT,EPOCH,DS50
      COMMON/C1/CK2,CK4,E6A,QOMS2T,S,TOTHRD,
     1          XJ3,XKE,XKMPER,XMNPDA,AE
      COMMON/C2/DE2RA,PI,PIO2,TWOPI,X3PIO2
      DATA IHG/1HG/
      DATA DE2RA,E6A,PI,PIO2,QO,SO,TOTHRD,TWOPI,X3PIO2,XJ2,XJ3,
     1     XJ4,XKE,XKMPER,XMNPDA,AE/.174532925E-1,1.E-6,
     2     3.14159265,1.57079633,120.0,78.0,.66666667,
     4     6.2831853,4.71238898,1.082616E-3,-.253881E-5,
     5     -1.65597E-6,.743669161E-1,6378.135,1440.,1./
      DIMENSION ISET(5)
      CHARACTER ABUF*80(2)
      DATA (ISET(I),I=1,5)/3HSGP,4HSGP4,4HSDP4,4HSGP8,4HSDP8/
*
*     SELECT EPHEMERIS TYPE AND OUTPUT TIMES
*
      CK2=.5*XJ2*AE**2
      CK4=-.375*XJ4*AE**4
      QOMS2T=((QO-SO)*AE/XKMPER)**4
      S=AE*(1.+SO/XKMPER)
    2 READ (5,700) IEPT, TS,TF,DELT
      IF(IEPT.LE.0) STOP
      IDEEP=0
*
*     READ IN MEAN ELEMENTS FROM 2 CARD T(TRANS) OR G(INTERN) FORMAT
*
      READ (5,706) ABUF
      DECODE(ABUF(1),707)    ITYPE
      IF(ITYPE.EQ.IHG) GO TO 5
      DECODE (ABUF,702) EPOCH,XNDT2O,XNDD6O,IEXP,BSTAR,IBEXP,XINCL,
     1       XNODEO,EO,OMEGAO,XMO,XNO
      GO TO 7
    5 DECODE(ABUF,701) EPOCH,XMO,XNODEO,OMEGAO,EO,XINCL,XNO,XNDT2O,
     1       XNDD6O,IEXP,BSTAR,IBEXP
    7 IF(XNO.LE.0.) STOP
      WRITE(6,704) ABUF,ISET(IEPT)
      IF(IEPT.GT.5) GO TO 900
      XNDD6O=XNDD6O*(10.**IEXP)
      XNODEO=XNODEO*DE2RA
      OMEGAO=OMEGAO*DE2RA
      XMO=XMO*DE2RA
      XINCL=XINCL*DE2RA
      TEMP=TWOPI/XMNPDA/XMNPDA
      XNO=XNO*TEMP*XMNPDA
      XNDT2O=XNDT2O*TEMP
      XNDD6O=XNDD6O*TEMP/XMNPDA
*
*     INPUT CHECK FOR PERIOD VS EPHEMERIS SELECTED
*     PERIOD GE 225 MINUTES   IS DEEP SPACE
*
      A1=(XKE/XNO)**TOTHRD
      TEMP=1.5*CK2*(3.*COS(XINCL)**2-1.)/(1.-EO*EO)**1.5
      DEL1=TEMP/(A1*A1)
      AO=A1*(1.-DEL1*(.5*TOTHRD+DEL1*(1.+134./81.*DEL1)))
      DELO=TEMP/(AO*AO)
      XNODP=XNO/(1.+DELO)
      IF((TWOPI/XNODP/XMNPDA) .GE. .15625) IDEEP=1
      BSTAR=BSTAR*(10.**IBEXP)/AE
      TSINCE=TS
      IFLAG=1
      IF(IDEEP .EQ. 1 .AND. (IEPT .EQ. 1 .OR. IEPT .EQ. 2
     1             .OR. IEPT .EQ. 4)) GO TO 800
    9 IF(IDEEP .EQ. 0 .AND. (IEPT .EQ. 3 .OR. IEPT .EQ. 5))
     1             GO TO 850
   10 GO TO (21,22,23,24,25), IEPT
   21 CALL SGP(IFLAG,TSINCE)
      GO TO 60
   22 CALL SGP4(IFLAG,TSINCE)
      GO TO 60
   23 CALL SDP4(IFLAG,TSINCE)
      GO TO 60
   24 CALL SGP8(IFLAG,TSINCE)
      GO TO 60
   25 CALL SDP8(IFLAG,TSINCE)
   60 X=X*XKMPER/AE
      Y=Y*XKMPER/AE
      Z=Z*XKMPER/AE
      XDOT=XDOT*XKMPER/AE*XMNPDA/86400.
      YDOT=YDOT*XKMPER/AE*XMNPDA/86400.
      ZDOT=ZDOT*XKMPER/AE*XMNPDA/86400.
      WRITE(6,705) TSINCE,X,Y,Z,XDOT,YDOT,ZDOT
      TSINCE=TSINCE+DELT
      IF(ABS(TSINCE) .GT. ABS(TF)) GO TO 2
      GO TO 10
  700 FORMAT(I1,3F10.0)
  701 FORMAT(29X,D14.8,1X,3F8.4,/,6X,F8.7,F8.4,1X,2F11.9,1X,F6.5,I2,
     1       4X,F8.7,I2)
  702 FORMAT(18X,D14.8,1X,F10.8,2(1X,F6.5,I2),/,7X,2(1X,F8.4),1X,
     1       F7.7,2(1X,F8.4),1X,F11.8)
  703 FORMAT(79X,A1)
  704 FORMAT(1H1,A80,/,1X,A80,//,1X,A4,7H TSINCE,
     1 14X,1HX,16X,1HY,16X,1HZ,14X,
     1 4HXDOT,13X,4HYDOT,13X,4HZDOT,//)
  705 FORMAT(7F17.8)
  706 FORMAT(A80)
  707 FORMAT(79X,A1)
  930 FORMAT("SHOULD USE DEEP SPACE EPHEMERIS")
  940 FORMAT("SHOULD USE NEAR EARTH EPHEMERIS")
  950 FORMAT("EPHEMERIS NUMBER",I2," NOT LEGAL, WILL SKIP THIS CASE")
  800 WRITE(6,930)
      GO TO 9
  850 WRITE(6,940)
      GO TO 10
  900 WRITE(6,950) IEPT
      GO TO 2
      END

ACTAN

      FUNCTION  ACTAN(SINX,COSX)
      COMMON/C2/DE2RA,PI,PIO2,TWOPI,X3PIO2
      ACTAN=0.
      IF (COSX.EQ.0. ) GO TO 5
      IF (COSX.GT.0. ) GO TO 1
      ACTAN=PI
      GO TO 7
    1 IF (SINX.EQ.0. ) GO TO 8
      IF (SINX.GT.0. ) GO TO 7
      ACTAN=TWOPI
      GO TO 7
    5 IF (SINX.EQ.0. ) GO TO 8
      IF (SINX.GT.0. ) GO TO 6
      ACTAN=X3PIO2
      GO TO 8
    6 ACTAN=PIO2
      GO TO 8
    7 TEMP=SINX/COSX
      ACTAN=ACTAN+ATAN(TEMP)
    8 RETURN
      END

FMOD2P

      FUNCTION FMOD2P(X)
      COMMON/C2/DE2RA,PI,PIO2,TWOPI,X3PIO2
      FMOD2P=X
      I=FMOD2P/TWOPI
      FMOD2P=FMOD2P-I*TWOPI
      IF(FMOD2P.LT.0) FMOD2P=FMOD2P+TWOPI
      RETURN
      END

THETAG

      FUNCTION THETAG(EP)
      COMMON /E1/XMO,XNODEO,OMEGAO,EO,XINCL,XNO,XNDT2O,XNDD6O,BSTAR,
     1 X,Y,Z,XDOT,YDOT,ZDOT,EPOCH,DS50
      DOUBLE PRECISION EPOCH,D,THETA,TWOPI,YR,TEMP,EP,DS50
      TWOPI=6.28318530717959D0
      YR=(EP+2.D-7)*1.D-3
      JY=YR
      YR=JY
      D=EP-YR*1.D3
      IF(JY.LT.10) JY=JY+80
      N=(JY-69)/4
      IF(JY.LT.70) N=(JY-72)/4
      DS50=7305.D0 + 365.D0*(JY-70) +N + D
      THETA=1.72944494D0 + 6.3003880987D0*DS50
      TEMP=THETA/TWOPI
      I=TEMP
      TEMP=I
      THETAG=THETA-TEMP*TWOPI
      IF(THETAG.LT.0.D0) THETAG=THETAG+TWOPI
      RETURN
      END