$STORAGE:2
      PROGRAM OSP
C     NUMERICAL ORBIT INTEGRATION & SEARCH PROGRAM
C     UNIVERSAL VERSION FOR BOTH LOW AND HIGH ORBITTING SATELLITES
C     AUTOMATICALLY SEARCHES FOR THE PASS OF REQUIRED MIN.ELEVATION
C     GIVEN BY VALUES "ELMIN"
C
C     6x6 GRAVITY FIELD IS USED FOR HIGH SATELLITES (h > 4000km)
C     18x18 GRAVITY FIELD IS USED FOR OTHER SATELLITES
C
C     INTEGRATION WRITTEN BY A.T.SINCLAIR. (ROYAL GREENWICH OBSERVATORY)
C     NUMERICAL INTEGRATION BY GAUSS-JACKSON METHOD
C     REST BY I.PROCHAZKA, CVUT, PRAGUE
C
C
C      INPUT FILES
C            GEM10N.DAT      NORMALISED GRAVITY FIELD
C            SLRCOR.DAT      STATION COORD.FILE
C            IRV_****.DAT    STARTING VALUES AT 0 HOURS EACH DAY
C                            (FROM UNIVERSITY OF TEXAS)
C                             **** ARE THE FOUR SAT.ID.DIGITS
C      OUTPUT FILES
C            NFOUT.DAT      SUMMARY LISTING OF INTEGRATION OUTPUT
C            GEOC.DAT       UNFORMATTED INTEGRATION OUTPUT
C
C
      DOUBLE PRECISION BA(6),Y(3,9),YD(3,9),
     1 T,H,HMIN,YT(3),YDT(3),BAOUT(6),
     2      J2,GM,DAE,SMD,TSTART,STMJD,STMJDR,SATNO,
     3  RDSEC0,RADSEC,THD,CTH,STH,TEMP,DSEC,RMJD,FMJD,AZ,EL
      REAL TC(20,20),TS(20,20),ERR(3),ERRB(5),SAT(2),JJ(20),NMS(2),
     1     ELMIN(2)
      CHARACTER*12 FNAME
      CHARACTER *10 SNAM
      CHARACTER *2 SNAM2
      COMMON/TOM/TC,TS,SMD,TSTART,NZON,NMAX,MMAX,
     1  JJ,J2,GM,DAE,THD,RADSEC
      DATA RDSEC0/7.2921151463D-5/
C     MIN.ELEVATION FOR HIGH/LOW SATELLITES TO BE TABULATED
      DATA ELMIN /20.,0./
C
C      READ GIVING START OF PASS
C      CALCULATE STMJDR : 3 MINS EARLIER THAN START TIME.
C      CALCULATE STMJD  : 0H MJD AT LEAST 9 MINS EARLIER
C                     THAN STMJDR (TO ENSURE THAT STMJDR IS NOT IN
C                     THE STARTING ROUTINE)
C
      OPEN(UNIT=1,FILE='CON',STATUS='UNKNOWN')
      WRITE(1,950)
950   FORMAT(' ORBIT INTEGRATION PROGRAM WITH AUTOMATICAL PASS SEARCH '/
     #       ' ****************************************************** '/
     #   //  ' The starting values for the search :'/
     #       '       0 - current time '/
     #       '       1 - input values        >>> '\)
      READ(1,*) JUMP
      IF(JUMP.NE.0) GOTO 9502
      CALL GETDAT(IYS,IMS,IDS)
      IYS=IYS-1900
      CALL GETTIM(IHS,IMNS,ISEC,IS100)
      IF(IHS.EQ.0.AND.IMNS.LT.13) WRITE(1,9505)
9505  FORMAT(' MIDNIGHT PROBLEM, ENTER THE STARTING TIME MANUALLY')
      IF(IHS.EQ.0.AND.IMNS.LT.13) GOTO 9502
      GOTO 9504
9502  WRITE(1,9503)
9503  FORMAT('       date         (yy,mm,dd) >>> '\)
      READ(1,*) IYS,IMS,IDS
      WRITE(1,951)
951   FORMAT('       time  > 12 min  (hh,mm) >>> '\)
      READ(1,*) IHS,IMNS
C
9504  NPASS=1
      IOFF = 0
      WRITE(1,9501)
9501  FORMAT('       satellite name - 2 char >>> '\)
      READ(1,'(1A2)') SNAM2
      KODE=4
      CALL CATLG(SATNO,ISAT4,SNAM,SNAM2,KODE,NPTSEC)
      WRITE(1,9506) SNAM,SATNO, IYS,IMS,IDS,IHS,IMNS
9506  FORMAT( ' Satellite        :  ',1A10,2X,F9.0,/,
     #        ' Search starts on :  ',I4,1H/,I2,1H/,I2,' at ',I2,1H:,I2)
      ISAT=SATNO/1000.
      WRITE(FNAME,'(F8.0)') SATNO
      READ(FNAME,'(2A4)') NMS
      CALL MODJD(IYS,IMS,IDS,STMJDR)
C
C     ASSIGN THE IRV DATA FILE NAME
C
      WRITE(FNAME,4121) ISAT
4121  FORMAT(4HIRV_,I4,4H.DAT)
C
      STMJDR=STMJDR +(IHS*60 +IMNS -3)/1440.D0
      STMJD =IDINT(STMJDR -9.D0/1440.D0 -40000.D0) +40000.D0
C
      DO 300 L=2,18
      JJ(L)=0.0
      DO 300 M=1,L
      TC(L,M)=0.0
      TS(L,M)=0.0
 300  CONTINUE
      NZON=18
      NMAX=18
      MMAX=18
C
C      OPEN FILE CONTAINING NORMALISED GEM10N GEOPOTENTIAL
C      COEFFICIENTS.
C
      OPEN(UNIT=2,FILE='GEM10N.DAT',STATUS='OLD')
      READ (2,908)J2
 908      FORMAT(10X,D15.8)
      J2=-J2*DSQRT(5.D0)
      JJ(2)=SNGL(J2)
 303      READ(2,900,END=302) L,M,CC,SS
 900      FORMAT(6X,2I2,2E15.8)
      IF(L.GT.18) GO TO 303
      IF(M.GT.18) GO TO 303
      IF(M.EQ.0)GO TO 304
      TC(L,M)=CC
      TS(L,M)=SS
      GO TO 303
 304      JJ(L)=-CC*SQRT(2.0*FLOAT(L)+1.0)
      GO TO 303
 302      CONTINUE
      CLOSE(UNIT=2)
C
C     GM    : GM OF EARTH IN (METERS**3)/(SEC**2)
C     DAE   : EQUATORIAL RADIUS OF EARTH IN MEGAMETERS.
C     HMIN  : STEP LENGTH IN MINUTES
C     NST   : NUMBER OF STEPS (25 HOUR SPAN INTEGRATION)
C     SMD   : MJD OF BASE EPOCH
C     TSTART: START TIME OF INTEGRATION (MINUTES AFTER BASE EPOCH)
C     NA    : NUMBER OF STARTING ITERATIONS (DEFAULT 5)
C     NOCORR: SET TO 1 TO SUPPRESS G-J CORRECTOR.
C
      HMIN=1.D0
      NST=1510
      NA=4
      NOCORR=1
      GM=0.39860044D15
      DAE=6.378145D0
C
C
C      READ FILE OF STARTING VALUES AND LOCATE THOSE FOR
C      REQUIRED DAY.
C
C
C
C      TEXAS STARTING VALUES
C
C      REFERRED TO GREENWICH MERIDIAN OF UT1 AND TRUE EQUATOR
C      LEAP SECOND IS INCLUDED
C      DATA IN METERS, METERS/SEC
C
 305      CONTINUE
          OPEN(UNIT=2,FILE=FNAME,STATUS='OLD')
   5      READ(2,930) JUNK
 930      FORMAT(A1)
          READ(2,931) IY,IMN,ID,IH,IM,SEC,(BAOUT(K),K=1,3)
 931      FORMAT(1X,I4,4I3,F5.1,3F18.6)
          READ(2,932) IRVID,IRVSQ,IDYSQ,(BAOUT(K),K=4,6)
 932      FORMAT(I5,I6,I6,5X,3F18.9)
      READ(2,933) IXP,IYP,IRATE
 933      FORMAT(3(1X,I6))
      IY=IY -1900
      CALL MODJD(IY,IMN,ID,SMD)
      L=IDINT(STMJD -SMD)
      IF(L.LT.0) GO TO 97
      IF(L.EQ.0) GO TO 106
      DO 107 LL=1,L
      READ(2,930) JUNK
      READ(2,931) IY,IMN,ID,IH,IM,SEC,(BAOUT(K),K=1,3)
      READ(2,932) IRVID,IRVSQ,IDYSQ,(BAOUT(K),K=4,6)
107   READ(2,933) IXP,IYP,IRATE
      IY=IY -1900
      CALL MODJD(IY,IMN,ID,SMD)
      IF(SMD.NE.STMJD) GO TO 96
106   CONTINUE
      CLOSE(UNIT=2)
      TSTART=DBLE(FLOAT(IH*60 +IM))
C
C     CEHCK AND DISPLAY THE IRVs USED
C
      IF(ISAT4.NE.IRVID) WRITE(1,9510) ISAT4,IRVID
9510  FORMAT(/' *** ERROR IN IRV DATA SETS COMPATIBILITY *** ',2I6)
      IF(ISAT4.NE.IRVID) STOP
      WRITE(1,9511) IRVSQ,IDYSQ
9511  FORMAT(' IRV data set #',I6,'  day within the set #',I6,' used')
      WRITE(1,934)IY,IMN,ID,BAOUT
934   FORMAT(3I4,/,3F13.2,/,3F13.5)
C
C     COMPUTE THE SATELLITE HEIGHT ABOVE THE EARTH
C
      R2=BAOUT(1)**2+BAOUT(2)**2+BAOUT(3)**2
      R=SQRT(R2)
      H=(R-6378000.)/1000.
      LOW=1
      IF(H.GT.4000.) LOW=0
C
C     DECREASE THE NUMBER OF GRAVITY COEFFICIENTS USED FOR HIGH SATELLITES
C
      IF(LOW.EQ.0) NZON=7
      NMAX=NZON
      MMAX=NZON
C
      IF(SEC.NE.0.0) GO TO 98
      RADSEC=RDSEC0 +DBLE(FLOAT(IRATE))*1.D-14
      BAOUT(4)=BAOUT(4) -RADSEC*BAOUT(2)
      BAOUT(5)=BAOUT(5) +RADSEC*BAOUT(1)
      CALL STEQ(SMD,TSTART,THD,EQE)
      CTH=DCOS(THD)
      STH=DSIN(THD)
      BA(1)=(BAOUT(1)*CTH -BAOUT(2)*STH)*1.D-6
      BA(2)=(BAOUT(1)*STH +BAOUT(2)*CTH)*1.D-6
      BA(3)=BAOUT(3)*1.D-6
      BA(4)=(BAOUT(4)*CTH -BAOUT(5)*STH)*0.0864D0
      BA(5)=(BAOUT(4)*STH +BAOUT(5)*CTH)*0.0864D0
      BA(6)=BAOUT(6)*0.0864D0
 306      CONTINUE
      OPEN(UNIT=101,FILE='NFOUT.DAT',STATUS='UNKNOWN')
      WRITE(101,903)HMIN,NST,SMD,TSTART
 903      FORMAT(6H STEP=,F4.1,20H MINS   NO.OF STEPS=,I6/
     1   11H BASE EPOCH,F8.1/11H START TIME,F8.1,5H MINS/)
      WRITE(101,905) BA
 905      FORMAT(1H ,3D20.12/1H ,3D20.12/)
      WRITE(101,915)GM,NMAX,MMAX,NZON
 915      FORMAT(1H ,'GM = ',D16.9/' NMAX, MMAX, NZON =',3I4)
      IF(NOCORR.EQ.1)WRITE(101,916)
      IF(NOCORR.NE.1)WRITE(101,917)
 916      FORMAT(/' CORRECTOR OMITTED')
 917      FORMAT(/' CORRECTOR INCLUDED')
C
C      OPEN OUTPUT FILE GEOC.DAT
C      TEMP IS INTERVAL IN MINS OF INPUT START TIME FROM
C      EPOCH OF STARTING VALUES.
C
      OPEN(UNIT=3,FILE='GEOC.DAT',STATUS='UNKNOWN',
     1      FORM='UNFORMATTED')
      TEMP=3.D0 +AINT(0.5 +SNGL(STMJDR-SMD)*1440.0)
      WRITE(3) STMJDR,SMD,NMS,NPASS,IOFF,DUM,THD,RADSEC
C
C
      H=HMIN/1440.D0
      T=0.D0
      GM=GM*1.D-18
      GM=GM*86400.D0*86400.D0
C
C
C     STARTING PROCEDURE
C
      CALL STARTE(T,H,BA,NA,Y,YD,ERRB)
      WRITE(101,901) NA,(ERRB(I),I=1,3)
 901  FORMAT(' NUMBER OF STARTING ITERATIONS',I19/
     1' MAXIMUM TRUNCATION ERROR',22X,E9.2/
     2' MAXIMUM CONVERGENCE ERROR OF PENULT ITERATION',E10.2/
     3' PREDICTED CONVERGENCE ERROR OF LAST ITERATION',E10.2//)
C
C     MAIN INTEGRATION
C
      EMAX=0.0
      ISTEP=0
      ELM=ELMIN(LOW+1)
      T1=(STMJDR -SMD)*1440.0
      WRITE(1,971)
 971  FORMAT(/'   TIME     AZIM    ELEV '/
     1        '  hr min       (deg)     ')
      DO 103 I=1,NST
      CALL GAUSSJ(T,H,YT,YDT,NOCORR,ERR)
      IF(ABS(ERR(1)).GT.EMAX) EMAX=ERR(1)
      IF(ABS(ERR(2)).GT.EMAX) EMAX=ERR(2)
      IF(ABS(ERR(3)).GT.EMAX) EMAX=ERR(3)
      ST=AINT(0.5 +SNGL(T*1440.D0 +TSTART))
      DSEC=DBLE(ST)*60.D0
      RMJD=STMJD
      FMJD=ST/1440.D0
      CALL STPS(YT(1),YT(2),YT(3),RMJD,FMJD,AZ,EL,ISHAD,IDYTM)
      MIN=ST
      IHR=MIN/60
      MIN=MIN-IHR*60
961   FORMAT(2I4,2F8.1,4X,A2,2X,A2)
      IF(ST.LT.T1) GO TO 103
      IF(EL.LT.ELM.AND.ISTEP.EQ.0) GO TO 103
      IF(EL.LT.ELM.AND.ISTEP.GT.0) GO TO 104
      ISTEP=ISTEP+1
      IF(ISTEP.GT.50) GOTO 104
      WRITE(1,961) IHR,MIN,AZ,EL,IDYTM,ISHAD
      WRITE(101,907) ST,YT(1),YT(2),YT(3)
 907  FORMAT(1H ,F11.2,3F15.10)
      WRITE(3) DSEC,YT(1),YT(2),YT(3)
 103  CONTINUE
 104  CONTINUE
      IF(ISTEP.EQ.0)
     1 WRITE(1,*) ' *** No point above the minimal elevation *** '
C      IF(ISTEP.GT.0) WRITE(1,970) ISTEP
C 970  FORMAT(/ ' *** Totally ',I4,' minutes tabulated '//)
      WRITE(101,913) EMAX
 913  FORMAT(///' MAXIMUM GAUSS-JACKSON ERROR =',E10.2///)
      CLOSE(UNIT=101)
      CLOSE(UNIT=3)
      STOP
 99   WRITE(1,*)' STARTING VALUES FOR REQUIRED DATE NOT IN FILE'
      STOP
 98   WRITE(1,*) ' STARTING VALUES EPOCH NOT 0.0 SEC'
      STOP
 97   WRITE(1,*) ' START DATE TOO EARLY'
      STOP
 96   WRITE(1,*) ' STARTING VALUE FILE NOT AT 1 DAY INTERVAL'
      STOP
      END
C
C      *********************************************************
C
      SUBROUTINE GAUSSJ(T,H,Y,YD,NOCORR,ERR)
C
C     INTEGRATION BY GAUSS-JACKSON, OF POSITION ONLY
C     VELOCITY IS NOT USED, AND IS NOT RETURNED
C
      DOUBLE PRECISION T,H,Y(3),YD(3),D1(3),D2(3),H2,BA,YY(6)
     1 ,F(3,9),P(9),EE(9),FA(6),Q(9),R,S,FF,SE,ERRD(3)
      REAL ERR(3)
      COMMON/START/F,D1,D2
      DATA P/103798439.D0, -385853488.D0, 867424848.D0,-1274515624.D0,
     1 1258146350.D0, -831418464.D0, 354064088.D0, -88091848.D0,
     1  9751299.D0/
      DATA R,RA,RB/159667200.D0, 292524.0,9458775.0/
      DATA EE/-9.D0,36.D0,-84.D0,126.D0,-126.D0,84.D0,-36.D0,9.D0,-1.D0/
      DATA Q/9751299.D0, 16036748.D0, -34806724.D0, 48315732.D0,
     1 -45851950.D0,29482676.D0,-12309348.D0,3017324.D0,-330157.D0/
      H2=H*H
      FAC=H2*RA/R
      IF(NOCORR.EQ.1)FAC=H2*RB/R
C
C     PREDICTOR FOR Y
C
      DO 100 K=1,3
      S=0.D0
      SE=0.D0
      DO 101 J=1,9
      JJ=10 -J
      FF=F(K,J)
      S=S +P(JJ)*FF
      SE=SE +EE(JJ)*FF
 101  CONTINUE
      ERRD(K)=SE
      BA=H2*(D2(K) +S/R)
      Y(K)=BA
      YY(K)=BA
 100  CONTINUE
C
C     CORRECTOR FOR Y, AND ERROR ESTIMATE
C
      T=T +H
      IF(NOCORR.EQ.1)GO TO 110
      CALL FIT(T,YY,FA,0)
      DO 104 K=1,3
      S=Q(1)*FA(K+3)
      DO 105 J=2,9
      JJ=11 -J
      FF=F(K,J)
      S=S +Q(JJ)*FF
 105  CONTINUE
      BA=H2*(D2(K) +S/R)
      Y(K)=BA
      YY(K)=BA
 104  CONTINUE
C
C     ADVANCE TABLE BY ONE STEP
C
 110  CONTINUE
      CALL FIT(T,YY,FA,1)
      DO 102 K=1,3
      ERR(K)=FAC*SNGL(ERRD(K) +FA(K+3))
      DO 103 J=1,8
 103  F(K,J)=F(K,J+1)
      F(K,9)=FA(K+3)
      D1(K)=D1(K) +FA(K+3)
      D2(K)=D2(K) +D1(K)
 102  CONTINUE
      RETURN
      END
C
C      *********************************************************
C
      SUBROUTINE FIT(T,X,F,IND)
      DOUBLE PRECISION X(6),F(6),T,GM,DR,DA,DAE,J2,JJ,DPH,DR2,GMR2
     1   ,TS,JED,GMS,GMP(2),DEL,DEL3,RP,RP3,DGMP
     2   ,JEDP,XP1,XP2,XP3,SMD,TSTART,FAC,AUMM
     3   ,THDZ,RADSEC
      REAL J(20),C(20,20),S(20,20),CC(20),CS(20),PD(29),RT(57)
     1 ,XS(3,2),XM(3,2)
      COMMON/TOM/C,S,SMD,TSTART,NZON,NMAX,MMAX,J,J2,GM,DAE
     1 ,THDZ,RADSEC
      DATA NUM/0/,JEDP/0.0D0/
     1  ,GMS/9.906930412D11/,AUMM/149597.87D0/
      IF(NUM.EQ.1) GO TO 400
      AE=SNGL(DAE)
      DO 401 I=1,57
      RT(I)=SQRT(FLOAT(I))
 401  CONTINUE
      NUM=1
      GMP(1)=GMS
      GMP(2)=GM/81.30068D0
 400  CONTINUE
C
C
      DR2=X(1)**2 +X(2)**2 +X(3)**2
      GMR2=GM/DR2
      DR=DSQRT(DR2)
      R2=SNGL(DR2)
      R=SNGL(DR)
      DPH=X(3)/DR
      PH=SNGL(DPH)
C
C     COMPUTE ANGLE OF GREENWICH MERIDIAN FROM INERTIAL EQUINOX
C
      TH=SNGL(THDZ +RADSEC*T*86400.D0)
      SY = SNGL(X(2))
      SX = SNGL(X(1))
      CTH = COS(TH)
      STH = SIN(TH)
      CCL=(SX*CTH +SY*STH)/R
      CSL=(SY*CTH -SX*STH)/R
      CC(1)=CCL
      CS(1)=CSL
      DO 104 M=2,MMAX
      CC(M) = CC(M-1)*CCL -CS(M-1)*CSL
  104 CS(M) = CS(M-1)*CCL +CC(M-1)*CSL
      PD(1) = 1.0/RT(2)
      PD(2) = 0.0
      BET=0.0
      GAM=0.0
      AR=AE/R
      ARL=AR
      PB=1.0
      PC=PH
      PDD=3.0*PH
      A=0.0
      B=0.0
      DO 100 LZ=2,NMAX
      L=LZ
      PA=PB
      PB=PC
      PDC=PDD
      AL=FLOAT(L)
      PC=((2.0*AL-1.0)*PH*PB -(AL-1.0)*PA)/AL
      PDD=(AL+1.0)*PC +PH*PDC
      ARL=ARL*AR
      PLM1 = PDC
      LL=MIN0(L,MMAX)
      DO 101 M=1,LL
      AM = FLOAT(M)
      PLM = PLM1/(RT(L+M)*RT(L-M+1))
      IF(M.EQ.L)GO TO 102
      PLM1=RT(L+M)*RT(L-M)*PD(M)
      IF(M.EQ.L-1) GO TO 112
      PLM1=PLM1 +PH*RT(L-M)*RT(L-M-1)*PD(M+1)
 112  CONTINUE
      PD(M)=PLM
      IF(M.EQ.MMAX) PD(M+1)=PLM1/(RT(L-M)*RT(L+M+1))
      GO TO 103
  102 PLM1=0.0
      PD(M)=PLM
      PD(M+1)=0.0
  103 CONTINUE
      CLM=C(L,M)
      SLM=S(L,M)
      ARLL=ARL*RT(2)*RT(2*L +1)
      CCSS=CLM*CC(M)+SLM*CS(M)
      PL1M1=PH*PLM1 +(AL+AM+1.0)*PLM
      A=A -ARLL*PL1M1*CCSS
      B=B -ARLL*PLM1*CCSS
      IF(M.EQ.1)GO TO 105
      BET=BET +ARLL*AM*PLM*(CLM*CC(M-1) +SLM*CS(M-1))
      GAM=GAM +ARLL*AM*PLM*(CLM*CS(M-1) -SLM*CC(M-1))
      GO TO 106
  105 BET=BET +ARLL*PLM*CLM
      GAM=GAM -ARLL*PLM*SLM
  106 CONTINUE
  101 CONTINUE
      IF(L.GT.NZON)GO TO 100
      AJ=J(L)*ARL
      IF(IND.EQ.0) GO TO 107
 107  IF(L.EQ.2) GO TO 100
      A=A +AJ*PDD
      B=B +AJ*PDC
 100  CONTINUE
      JJ=J2*DAE*DAE/DR2
      DA=(1.D0 -JJ*(7.5D0*DPH*DPH -1.5D0)-A)*GMR2/DR
      B=GMR2*B
      BET = GMR2*BET
      GAM = GMR2*GAM
      F(1)=X(4)
      F(2)=X(5)
      F(3)=X(6)
      F(4)=-DA*X(1) +BET*CTH +GAM*STH
      F(5)=-DA*X(2) +BET*STH -GAM*CTH
      F(6)=-DA*X(3) -B -GMR2*JJ*3.D0*DPH
C
C
C      LUNAR & SOLAR PERTURBATIONS. ESTIMATE OF ET-UTC
C
      JED=SMD +TSTART/1440.D0 +T +2400000.5006D0
      IF(JED-JEDP.LT.0.05D0) GO TO 250
      CALL LSXYZ(JED,XS,XM)
      JEDP=JED
 250      CONTINUE
      Q=(JED -JEDP)/0.05
      DO 20 I=1,2
      DGMP=GMP(I)
      IF(I.EQ.2)GO TO 22
      XP1=DBLE(XS(1,1) +XS(1,2)*Q)
      XP2=DBLE(XS(2,1) +XS(2,2)*Q)
      XP3=DBLE(XS(3,1) +XS(3,2)*Q)
      GO TO 21
 22      XP1=DBLE(XM(1,1) +XM(1,2)*Q)
      XP2=DBLE(XM(2,1) +XM(2,2)*Q)
      XP3=DBLE(XM(3,1) +XM(3,2)*Q)
  21      DEL=DSQRT((X(1)-XP1)**2+(X(2)-XP2)**2+(X(3)-XP3)**2)
      DEL3=DEL*DEL*DEL
      RP=DSQRT(XP1*XP1 + XP2*XP2 + XP3*XP3)
      RP3=RP*RP*RP
      F(4)=F(4)+DGMP*((XP1-X(1))/DEL3 -XP1/RP3)
      F(5)=F(5)+DGMP*((XP2-X(2))/DEL3 -XP2/RP3)
      F(6)=F(6)+DGMP*((XP3-X(3))/DEL3 -XP3/RP3)
C
C      SOLAR RADIATION PRESSURE(FOR I=1).NO EARTH SHADDOW EFFECT.
C
      IF(I.EQ.2)GO TO 20
      FAC=3.70286D-9*7.46496D3*(AUMM**2)/DEL3
      F(4)=F(4)+FAC*(X(1)-XP1)
      F(5)=F(5)+FAC*(X(2)-XP2)
      F(6)=F(6)+FAC*(X(3)-XP3)
  20      CONTINUE
C
      IF(IND.EQ.0) GO TO 108
 108  RETURN
      END
C
C      ********************************************************
C
      SUBROUTINE STARTE(T,H,BA,NA,Y,YD,ERRB)
      DOUBLE PRECISION T,H,BA(6),YD(3,9),Y(3,9),F(3,9),D1(3),D2(3)
     1  ,CONST, FA(6),BB(6),H2,C(10,9),TA,BC,AA(9),AB(9),XP(3),XPP(3)
     2 ,CU,CS(40),CT(50)
      REAL ERRB(5)
      COMMON/START/F,D1,D2
      DATA AB/-1.D0,9.D0,-36.D0,84.D0,-126.D0,126.D0,-84.D0,36.D0,
     1   -9.0D0/
      DATA AA/ 9751299.D0,   -88091848.D0, 354064088.D0, -831418464.D0,
     1      1258146350.D0, -1274515624.D0,867424848.D0, -385853488.D0,
     2       103798439.D0/
      DATA CS/    9751299.D0, 16036748.D0,-34806724.D0, 48315732.D0,
     1  -45851950.D0, 29482676.D0,-12309348.D0,3017324.D0,
     2    -330157.D0,  -330157.D0, 12722712.D0,4151096.D0,
     3   -7073536.D0,  6715950.D0, -4252168.D0,1749488.D0,
     4    -423696.D0,    45911.D0,    45911.D0, -743356.D0,
     5   14375508.D0,   294572.D0, -1288750.D0,  931164.D0,
     6    -395644.D0,    96692.D0,  -10497.D0, -10497.D0,
     7     140384.D0, -1121248.D0, 15257256.D0, -1028050.D0,
     8      33872.D0,    49416.D0,  -17752.D0,     2219.D0,
     9       2219.D0,   -30468.D0,  220268.D0, -1307644.D0/
      DATA CT/   15536850.D0, -1307644.D0,  220268.D0,   -30468.D0,
     1       2219.D0,     2219.D0,  -17752.D0,    49416.D0,
     2              33872.D0, -1028050.D0,15257256.D0, -1121248.D0,
     3             140384.D0,   -10497.D0,  -10497.D0,    96692.D0,
     4            -395644.D0,   931164.D0,-1288750.D0,   294572.D0,
     5           14375508.D0,  -743356.D0,   45911.D0,    45911.D0,
     6            -423696.D0,  1749488.D0,-4252168.D0,  6715950.D0,
     7           -7073536.D0,  4151096.D0,12722712.D0,  -330157.D0,
     8            -330157.D0,  3017324.D0, -12309348.D0, 29482676.D0,
     9          -45851950.D0, 48315732.D0,-34806724.D0, 16036748.D0,
     1            9751299.D0,-45820566.D0,-94047140.D0,153921548.D0,
     2         -198129492.D0,182110720.D0, -115111084.D0, 47557620.D0,
     3          -11575388.D0,  1260182.D0/
      DATA CONST/ 159667200.D0/
      IF(NA.EQ.0) NA=5
      IF(NA.LT.3) NA=3
      IF(NA.GT.10) NA=10
      H2=H*H
      TA=T
      DO 101 I=1,10
      DO 120 J=1,9
      JJ=(I-1)*9 +J
      IF(JJ.LE.40) CU=CS(JJ)
      IF(JJ.GE.41) CU=CT(JJ-40)
      C(I,J)=(CU)/CONST
 120  CONTINUE
 101  CONTINUE
      CALL FIT(T,BA,FA,1)
      DO 100 I=1,3
      Y(I,1)=BA(I)
      YD(I,1)=BA(I+3)
      F(I,1)=FA(I+3)
      DO 100 J=1,3
 100  CONTINUE
C
C     APPROX INTEGRATION
C
      DO 102 I=1,3
      D2(I)=BA(I)/H2 -FA(I+3)/12.D0
 102  D1(I)=BA(I+3)/H +FA(I+3)/2.D0
      DO 103 J=2,9
      TA=TA +H
      DO 104 I=1,3
      D2(I)=D2(I) +D1(I)
 104  BB(I)=H2*D2(I)
      CALL FIT(TA,BB,FA,0)
      DO 105 I=1,3
      F(I,J)=FA(I+3)
      D1(I)=D1(I) +FA(I+3)
 105  CONTINUE
 103  CONTINUE
C
C    ITERATE TO IMPROVE ACCURACY
C
      DO 106 L=1,NA
      TA=T
C
C  COMPUTE INITIAL 2ND AND 1ST TERMS
      DO 107 I=1,3
      D2(I)=BA(I)/H2
      D1(I)=BA(I+3)/H
      DO 108 K=1,9
      D2(I)=D2(I) -C(1,K)*F(I,K)
      D1(I)=D1(I) -C(10,K)*F(I,K)
 108  CONTINUE
 107  CONTINUE
      DO 109 J=2,9
      TA=TA +H
C
C  INCREMENT 2ND SUMS AND FORM Y AT STEP J
      DO 110 I=1,3
      D2(I)=D2(I) +D1(I)
      BC=H2*D2(I)
      DO 111 K=1,9
 111  BC=BC +H2*C(J,K)*F(I,K)
 110  BB(I)=BC
C
C  COMPUTE FUNCTION VALUES AND INCREMENT 1ST SUMS
      N=0
      IF(L.EQ.NA) N=1
      CALL FIT(TA,BB,FA,N)
      DO 112 I=1,3
      D1(I)=D1(I) +FA(I+3)
      F(I,J)=FA(I+3)
      IF(L.NE.NA) GO TO 112
      Y(I,J)=BB(I)
      YD(I,J)=BB(I+3)
      DO 125 JJ=1,3
 125  CONTINUE
 112  CONTINUE
 109  CONTINUE
      IF(L.NE.NA-2) GO TO 124
      DO 123 I=1,3
 123  XPP(I)=BB(I)
 124  CONTINUE
      IF(L.NE.NA-1) GO TO 106
      DO 121 I=1,3
      XP(I)=BB(I)
 121  CONTINUE
 106  CONTINUE
      DIFF=0.0
      DIFFP=0.0
      DO 122 I=1,3
      DIFFA=ABS(SNGL(Y(I,9) -XP(I)))
      DIFFB=ABS(SNGL(XP(I) -XPP(I)))
      IF(DIFFB.GT.DIFFP) DIFFP=DIFFB
 122  IF(DIFFA.GT.DIFF) DIFF=DIFFA
C
C  COMPUTE COORDS AT STEP 10 USING GAUSS-JACKN PREDICTOR
      DO 114 I=1,3
      D2(I)=D2(I) +D1(I)
      BC=H2*D2(I)
      DO 115 K=1,9
 115  BC=BC +H2*AA(K)*F(I,K)/CONST
 114  BB(I)=BC
C  COMPUTE FUNCTION VALUES AT STEP 10
C  COMPUTE ERROR/STEP DUE TO TRUNCATION OF VELOCITY FORMULA
C
      CALL FIT(TA+H,BB,FA,0)
      ERRB(1)=0.0
      ERRB(4)=0.0
      DO 116 I=1,3
      BC=FA(I+3)
      DO 117 K=1,9
      BC=BC +AB(K)*F(I,K)
 117  CONTINUE
      ERR=ABS(SNGL(H2*BC)*0.007)
      IF(ERR.GT.ERRB(1)) ERRB(1)=ERR
 116  CONTINUE
      ERRB(2)=DIFF
      ERRB(3)=DIFF*DIFF/DIFFP
      T=TA
      ERRB(4)=0.0
      ERRB(5)=0.0
      RETURN
      END
C
C      ************************************************************
C
      SUBROUTINE LSXYZ(JED,XS,XM)
C
C      COMPUTES POLYNOMIALS FOR SUN & MOON COORDS.
C      MEAN EQUATOR AND EQUINOX OF DATE
C
      DOUBLE PRECISION JED,T,AJD
      REAL LP,MP,L,M,M2,XS(3,2),XM(3,2)
      DATA DGRD/0.0174532925/
C
      AJD=JED
      ITER=1
 100      CONTINUE
C
      T=(AJD-2415020.D0)/36525.D0
      EPS=(23.452-0.0130*T)*DGRD
      COSE=COS(EPS)
      SINE=SIN(EPS)
C
C      LUNAR COORDINATES.
C
      TH=DMOD(477198.849D0*T,360.D0)
      M=(296.10 +TH)*DGRD
      TH=DMOD(483202.025D0*T,360.D0)
      F =( 11.25 +TH)*DGRD
      TH=DMOD(445267.114D0*T,360.D0)
      D =(350.74 +TH)*DGRD
      TH=DMOD(481267.883D0*T,360.D0)
      L =(270.43 +TH)*DGRD
      TH=DMOD(35999.05D0*T,360.D0)
      MP=(358.48 +TH)*DGRD
      TH=DMOD(36000.769D0*T,360.D0)
      LP=(279.70 +TH)*DGRD
      D2=D +D
      M2=M +M
      F2=F +F
C
C      LUNAR COORDINATES
C
      DLM=L +(6.289*SIN(M) +0.214*SIN(M2))*DGRD
      DBM=(5.128*SIN(F) +0.281*SIN(M +F) +0.278*SIN(M -F))*DGRD
      DPI=0.95075 +0.05182*COS(M) +0.00283*COS(M2)
      R=6.378137/(DPI*DGRD)
      SINB=SIN(DBM)
      COSB=COS(DBM)
      SINLM=SIN(DLM)
      COSLM=COS(DLM)
      XM(1,ITER)=R*COSB*COSLM
      XM(2,ITER)=R*(COSB*SINLM*COSE -SINB*SINE)
      XM(3,ITER)=R*(COSB*SINLM*SINE +SINB*COSE)
C
C      SOLAR COORDINATES
C
      DLS=LP +0.0355*SIN(MP)
      CLS=COS(DLS)
      SLS=SIN(DLS)
      R=(1.00014 -0.01675*COS(MP) -0.00014*COS(MP+MP))*149597.87
      XS(1,ITER)=R*CLS
      XS(2,ITER)=R*COSE*SLS
      XS(3,ITER)=R*SINE*SLS
      IF(ITER.EQ.2) GO TO 101
      ITER=2
      AJD=AJD +0.05D0
      GO TO 100
C
C
 101      CONTINUE
      DO 102 I=1,3
      XM(I,2)=XM(I,2) -XM(I,1)
 102      XS(I,2)=XS(I,2) -XS(I,1)
      RETURN
      END
C
C      *************************************************************
C
      SUBROUTINE MODJD(I,J,K,DMJD)
C
C      COMPUTE MJD
C
      DOUBLE PRECISION DMJD
      L=(J -14)/12
      AMJD=AINT(365.25*FLOAT(I+L)) +(15048 +K +367*(J -2 -12*L)/12
     1  +(1202 -12*I -J)/1200 )
      DMJD=DBLE(AMJD)
      RETURN
      END
C
C      ****************************************************************
C
      SUBROUTINE STEQ(SMD,TSTART,THD,EQE)
C
C      COMPUTES NEWCOMBS MEAN SIDEREAL TIME AND EQUATION OF EQUINOX
C
      DOUBLE PRECISION SMD,TSTART,THD,T,UT,P2
      REAL N,NA,NB
      DATA P2/6.283185307179586477D0/
      DATA FA,FB/0.0774848,  0.23089572/
      DATA DA,DB/2.1771852,  0.21276871/
      DATA NA,NB/2.8719993, -0.00092422/
      DATA GA,GB/2.1234424,  0.01720197/
      DATA EA,EB/0.4091395, -6.21796E-9/
      DATA S1,S2,S3,S4,S5/-17.25,  0.21, -1.27, -0.20,  0.13/
C
C
      UT=SMD +TSTART/1440.D0
      T=(UT -15019.5D0)/36525.D0
      THD=UT +(23925.836D0 +8640184.542D0*T
     1      +0.0929D0*T*T)/86400.D0
      THD=DMOD(THD,1.D0)*P2
C
C
      D=(UT -44000.D0)
      N=NA +NB*D
      F2=(FA +FB*D)*2.0
      D2=(DA +DB*D)*2.0
      G=GA +GB*D
      EPS=EA +EB*D
      CE=COS(EPS)
      SA=SIN(N)
      SB=SIN(2.0*N)
      SC=SIN(F2 -D2 +2.0*N)
      SD=SIN(F2 +2.0*N)
      SG=SIN(G)
      DPSI=(SA*S1 +SB*S2 +SC*S3 +SD*S4 +SG*S5)*4.8481368E-6
      EQE=DPSI*CE
      RETURN
      END
C
      SUBROUTINE STPS(XS,YS,ZS,RMJD,FMJD,AZ,EL,ISHAD,IDYTM)
C     **************************************
C
C     CONVERTS THE SATELLITE X,Y,Z INERTIAL COORDINATES INTO TOPOCENTRIC
C     AZIMUTH,ELEVATION; ROUGH VALUES, ONLY.
C
C     INPUT      XS,YS,ZS .........    INERTIAL COORD.IN MEGAMETERS
C                RMJD,FMJD ........    INTEGER,FRACT. DAY OF EPOCH
C                SLRCOR.DAT .......    FILE WITH THE STATION COORD.
C
C     OUTPUT     AZ,EL ............    AZIMUTH,ELEVATION IN DEGREES
C                ISHAD,IDYTM ......    SHADOW AND DAYTIME FLAGS
C
C     REV 030592    IP
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION STAT(3)
      INTEGER NMST(3)
      DATA PI,TWPI,PI2
     #       /3.1415926535897D0, 6.283185307179D0, 1.5707963268D0/
      DATA IPASS /0/
C
C     GET STATION COORDINATES IN THE FIRST CALL
C
      IPASS=IPASS+1
      IF(IPASS.GT.1) GOTO 10
      OPEN(UNIT=22,FILE='SLRCOR.DAT',STATUS='OLD')
      READ(22,3) NMST
3     FORMAT(3A2)
      READ(22,*) DLOND
      READ(22,*) DP
      READ(22,*) HT
      READ(22,*) DAE
      READ(22,*) DF
      READ(22,*) (STAT(J),J=1,3)
      CLOSE(UNIT=22)
      PHI=DP/180.D0*PI
      DLAM=DLOND/180.D0*PI
      CPHI=DCOS(PHI)
      SPHI=DSIN(PHI)
      SLAM=DSIN(DLAM)
      CLAM=DCOS(DLAM)
10    CONTINUE
C
      CALL SIDTM(RMJD,FMJD,TETA)
      CT3=DCOS(TETA)
      ST3=DSIN(TETA)
C
C     EXAMINE THE SKY AND ILLUMINATION CONDITIONS
C
      EPOCH=RMJD+FMJD
      CALL SKY(EPOCH,TETA,CT3,ST3,XS,YS,ZS,STAT,ISHAD,IDYTM)
C
C     TRANSFORMS TO THE ROTATING EARTH
      XSF=CT3*XS+ST3*YS
      YSF=-ST3*XS+CT3*YS
      ZSF=ZS
C     TRANSFORM TO THE LOCAL X,Y,Z
      SX=XSF-STAT(1)
      SY=YSF-STAT(2)
      SZ=ZSF-STAT(3)
C     COMPUTE AZ,EL
      HX=SPHI*CLAM*SX+SPHI*SLAM*SY-CPHI*SZ
      HY=-SLAM*SX+CLAM*SY
      HZ=CPHI*CLAM*SX+CPHI*SLAM*SY+SPHI*SZ
C
      S=DSQRT(HX*HX+HY*HY+HZ*HZ)
      HX=HX/S
      HY=HY/S
      HZ=HZ/S
C
      SXX=DSQRT(1.D0-HZ*HZ)
      CXX=HZ
      SZD=SXX
      CALL QUADR(SXX,CXX,ZD)
      SXX=-HY/SZD
      CXX=HX/SZD
      CALL QUADR(SXX,CXX,AZ)
      ALT=PI2-ZD
      IF(ALT.GT.PI2) ALT=ALT-PI
C
C     AZIMUTH IS NOW EQUAL ZERO FOR THE SOUTH DIRECTION
      AZ=AZ+PI
      CALL DDMOD(AZ,TWPI)
      AZ=AZ*180.D0/PI
      EL=ALT*180.D0/PI
      RETURN
      END
C
      SUBROUTINE QUADR(SXX,CXX,XX)
C     ****************
C
C     MODIFIED FUNCTION DATAN2
C
      DOUBLE PRECISION SXX,CXX,XX,PI,PI2,TWPI,X
      DATA PI,TWPI/3.1415926535897D0,6.283185307179586D0/
      DATA PI2/1.5707963268D0/
      DFNCT(X)=DATAN(DSQRT(1.D0-X*X)/X)
C
      IF(SXX) 1,2,2
2     IF(CXX)3,5,4
4     XX=DFNCT(CXX)
      RETURN
3     XX=DFNCT(-CXX)
      XX=PI-XX
      RETURN
5     XX=PI2
      RETURN
1     IF(CXX) 6,7,8
8     XX=DFNCT(CXX)
      XX=TWPI-XX
      RETURN
6     XX=DFNCT(-CXX)
      XX=PI+XX
      RETURN
7     XX=3.D0*PI2
      RETURN
      END
C
C
      SUBROUTINE DDMOD(A,B)
C     ****************
C
C     MODIFEID FUNCTION DMOD
C
      DOUBLE PRECISION A,B,C
      C=DMOD(A,B)
      IF(C.LT.0.D0) C=C+B
      A=C
      RETURN
      END
C
C
      SUBROUTINE MJD(IYR,MO,IDY,RMJD)
C     **************
C     TRANSFERS THE YEAR:MONTH:DAY INTO MODIFIED JULIAN DATE
C
C                                        R.N.1979
C
      DOUBLE PRECISION RMJD
      IF (MO.GT.2) GO TO 1
      MO=MO+9
      IYR=IYR-1
      GO TO 2
1     MO=MO-3
2     RHLP=15078.+AINT(365.25*IYR)+AINT(30.6*MO+0.4)+IDY
      RMJD=RHLP
      RETURN
      END
C
C
      SUBROUTINE SIDTM(RMJD,FMJD,DTETA)
C     ****************
C
C     CALCULATES THE VALUE OF SID.TIME
C
C     INPUT : RMJD  ... INTEGER PART OF EPOCH
C             FMJD  ... FRACT.PART OF EPOCH   (DOUBLE)
C             C     OUTPUT: DTETA ... SID. ANGLE (RADIANS)
C
C                                        R.N.1979
      DOUBLE PRECISION FMJD,ST,CT,TWPI,RIMJD,T,DTETA
     #                 ,RMJD
      RIMJD=RMJD
      RIMJD=RIMJD-42412.D0
      DTETA=1.728501408D0+0.0172027915D0*RIMJD
      DTETA=DTETA+6.3003881D0*FMJD
      T=RIMJD+FMJD
      DTETA=DTETA-8.360947D-5*DSIN(-1.943525D0-9.2415D-4*T)
      CALL DDMOD(DTETA,6.2831853071796D0)
      RETURN
      END
C
      SUBROUTINE SKY(TIME,STG,CSTG,SSTG,SX,SY,SZ,STAT,ISHAD,IDYTM)
C     **************
C
C     COMPUTES THE SKY AND SHADOW CONDITION
C     SCHEME FROM THE SAO AIMLASER
C     INPUT TIME     ... EPOCH
C           STG,CSTG,SSTG SIDEREAL TIME AND ITS COS AND SIN
C           SX,SY,SZ ... SATELLITE COORDINATES IN THE IN.FRAME
C           STAT     ... STATION COORDINATES
C     OUTP. IDYTM,ISHAD
C
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION SUN(3),STAT(3),XS(3)
C
C     SUN 1 DEG. BELOW HORIZONT FOR NIGHT
      DATA DUSK /0.02D0/
C
C
C
      ISHAD=2HIL
      IDYTM=2HDY
      RS=0.D0
      DO 53 J=1,3
53    RS=RS+STAT(J)**2
      RS=DSQRT(RS)
C
      X=STAT(1)
      Y=STAT(2)
      XS(1)=X*CSTG-Y*SSTG
      XS(2)=Y*CSTG+X*SSTG
      XS(3)=STAT(3)
C
C     GET POSITION OF SUN
C
225   CALL SUNVC(TIME,SUN,RSUN)
      DO 226 K=1,3
226   SUN(K)=SUN(K)/RSUN
      CALL PRECS(SUN,SUN,TIME,1)
C
C     EXAMINE SKY CONDITION
C
      X=0.D0
      DO 228 K=1,3
228   X=X+SUN(K)*XS(K)
      X=X/RS
      IF(X+DUSK)229,229,230
230   RETURN
229   IDYTM=2HNI
C
C     SIMPLE SHADOW TEST FOR SPHERICAL EARTH (NIGHT ONLY)
C
      D2=((SX*SUN(2)-SY*SUN(1))**2+(SY*SUN(3)
     #    -SZ*SUN(2))**2+(SZ*SUN(1)-SX*SUN(3))**2)
C     EARTH RADIUS SQUARED
      IF(D2.LT.40.68D0) ISHAD=2HSH
C
      RETURN
      END
C
C
      SUBROUTINE PRECS (XZ,X,T,JAY)
C     ******************
C
C              XZ IS VECTOR REFERRED TO EQUINOX OF 1950.0
C              X  IS VECTOR REFERRED TO EQUINOX OF SMITHSONIAN DAY = T
C
C              JAY=1 MEANS XZ IS INPUT VECTOR AND X IS OUTPUT
C              JAY=2 MEANS X IS INPUT VECTOR AND XZ IS OUTPUT
C              JAY=3 MEANS REPEAT THE LAST CALL TO PRECS USING NEW
C                    VALUES OF COMPONENTS IN THE INPUT VECTOR.
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION V(3),W(3),X(3),XZ(3)
      IF(JAY-2)1,1,2
    1 KAY=JAY
      D=(T-33281.923D0)/36524.22D0
      D2=D*D
      D3=D2*D
      XX=1.-.00029696D0*D2-.00000013D0*D3
      YX=-.02234941D0*D-.00000676D0*D2+.00000221D0*D3
      ZX=-.00971690D0*D+.00000207D0*D2+.00000096D0*D3
      YY=1.-.00024975D0*D2-.00000015D0*D3
      ZY=-.00010858D0*D2
      ZZ=1.-.00004721D0*D2
    2 IF(KAY-1)3,3,6
C                             XZ IS INPUT
    3 DO 4 J=1,3
    4 W(J)=XZ(J)
      V(1)=ZX*W(3)+YX*W(2)+XX*W(1)
      V(2)=ZY*W(3)+YY*W(2)-YX*W(1)
      V(3)=ZZ*W(3)+ZY*W(2)-ZX*W(1)
      DO 5 J=1,3
    5 X(J)=V(J)
      RETURN
C                             X IS INPUT
    6 DO 7 J=1,3
    7 V(J)=X(J)
      W(1)=XX*V(1)-YX*V(2)-ZX*V(3)
      W(2)=YX*V(1)+YY*V(2)+ZY*V(3)
      W(3)=ZX*V(1)+ZY*V(2)+ZZ*V(3)
      DO 8 J=1,3
    8 XZ(J)=W(J)
      RETURN
C
      END
C
      SUBROUTINE SUNVC (T,X,R)
C     *****************
C
C
C            T= TIME IN SMITHSONIAN DAYS
C              X=VECTOR TO SUN REFERRED TO EQUINOX AND EQUATOR OF 1950.0
C              R=MAGNITUDE OF X VECTOR    (LENGTHS ARE IN A.U.)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(3)
C
      A=6.2291402D0+.0172021242D0*(T-39125.D0)
      E=A
      DO 1 J=1,6
    1 E=A+.016728643D0*DSIN(E)
      SE=DSIN(E)
      CE=DCOS(E)
      R=1.0000002D0-.016728643D0*CE
      X(1)=.20947919D0*CE+.97767652D0*SE-.003504301D0
      X(2)=-.89708519D0*CE+.19215772D0*SE+.015007012D0
      X(3)=.43367832D0*X(2)
      RETURN
      END
C
      SUBROUTINE CATLG(SATNM7, ISAT4, SNAM, SNAM2, KODE, NPTSEC)
C     *****************
C
C     ROUTINE TO INTERFACE THE CATALOG FILE OF THE SATELLITES
C
C     INPUT/OUTPUT       SATNM7 ...   7 DIGITS , DOUBLE PRECISION
C                                     SATELLITE 7 DIGIT INDICATOR
C                        ISAT4  ...   INTEGER  SATELLITE 4 DIGIT IND.
C                        SNAM   ...   CHARACTE*10, SATELLITE NAME
C                        SNAM2  ...   CHARACTER*2, ABBREV.SATEL.NAME
C      INPUT             KODE   ...   DESCRIBING, WHICH OF THE FORMER
C                                     FOUR IS INPUT
C      OUTPUT            KODE   ...   NEGATIVE IF SEARCH FAILED
C                        NPTSEC ...   INTEGER, NORMAL POINTS IN SECONDS
C
C      INPUT FILE      SATCAT.DAT IN FIXED FORMAT, STARTING AT LINE
C                      NUMBER 7, SATELLITE CATALOGUE FILE
C
      DOUBLE PRECISION SATNM7, DHELP
      CHARACTER*10 SNAM  ,SNAMF
      CHARACTER*2  SNAM2 ,SNAM2F, DUMMY
C
      OPEN(9,FILE='SATCAT.DAT',STATUS='OLD')
      DO 2 I=1,6
2     READ(9,1) DUMMY
1     FORMAT(1A2)
C
C     READ AND MATCH THE SATELLITE ID
C
10    READ(9,11,END=100) SNAMF,SNAM2F,DHELP,ISAT4F,NPTSEC
11    FORMAT(1A10,2X,1A2,2X,F7.0,9X,I4,9X,I4)
      IF(KODE.EQ.3.AND.SNAMF .EQ.SNAM ) GOTO 20
      IF(KODE.EQ.4.AND.SNAM2F.EQ.SNAM2) GOTO 20
      IF(KODE.EQ.1.AND.DHELP.EQ.SATNM7) GOTO 20
      IF(KODE.EQ.2.AND.ISAT4F.EQ.ISAT4) GOTO 20
      GOTO 10
C
100   WRITE(*,101)
101   FORMAT(' **** ERROR **** the followIng satellite ID not found ')
      WRITE(*,11) SNAM,SNAM2,SATNM7,ISAT4,KODE
      KODE=-1*KODE
      CLOSE(9)
      STOP 102
C
20    CONTINUE
      SNAM=SNAMF
      SNAM2=SNAM2F
      SATNM7=DHELP
      ISAT4=ISAT4F
C
      CLOSE(9)
      RETURN
      END

