$STORAGE:2
C
C
C
      PROGRAM N06M
C
C
C     EVALUATION OF TESSERAL AMPLITUDES
C     EVALUATION OF LUNAR AND SOLAR PERTURBATIONS
C     BASED ON N06M CODE BY A.T.SINCLAIR,RGO
C
C
C     INPUT FILES
C            N6MDAT.DAT      PASS START AND END TIMES
C            TESSER.DAT      NORMALISED GRAVITY FIELD COEFFICIENTS
C            SAOELM.DAT      ORBITAL ELEMENTS
C                        (FILE IS OPENED IN SBR READEL)
C     OUTPUT FILES
C            TESSPU.DAT      UNFORMATTED OUTPUT FILE
C            TESSPR.DAT      SUMMARY LISTING OF COMPUTATION
C
      REAL I,N,N1,K,LP,LPQ,AMP(6),
     1  J,N0,I0,I1,B(5,3),M0,M2,M3,ND,KAP,M1,NMS(2)
      DOUBLE PRECISION DEPOCH,DTIME,DKAP,DND,DDI,DEC,DES,DA,DDN,DW0
      DOUBLE PRECISION DW1,DN0,DN1,DI0,DI1,DE0,DE1,DM0,DM1,DRAD,DGM
     1   ,DDJ,TSTART,TEND,DT
      INTEGER P,Q,TDT(6)
      REAL GDS(24,9),GS(24,9),DE(2),DI(2),DM(2),DW(2),DN(2)
      COMMON/CBB/DW1,DN1,NMS
      COMMON /CBC/DRAD,DGM,DDJ,DT
C
      OPEN(UNIT=4,FILE='TESSPU.DAT',STATUS='UNKNOWN',
     #     FORM='UNFORMATTED')
      OPEN(UNIT=3,FILE='TESSPR.DAT',STATUS='UNKNOWN')
      OPEN(UNIT=2,FILE='N6MDAT.DAT',STATUS='OLD')
      OPEN(UNIT=1,FILE='CON',STATUS='UNKNOWN')
C
      DRAD=0.01745329252D0
      RAD=SNGL(DRAD)
C
C     UNITS   MEGAMETRES AND REVOLUTIONS/DAY
C
      DGM=75371.38D0
      GM=SNGL(DGM)
      AE=6.378155
      THD=1.002737812
      DDJ=0.0660644D0
      J=SNGL(DDJ)
C
C     READ PASS INFORMATION FILE
C
      READ(2,900)IYS,IMS,IDS,IHS,IMNS,IYF,IMF,IDF,IHF,IMNF,
     1   NMS,NPASS,IOFF,DM2
 900  FORMAT(2(5I2,1X),2A4,1X,A2,I5,E10.2)
      CLOSE(UNIT=2)
C
C
      CALL MODJD(IYS,IMS,IDS,IHS,IMNS,TSTART)
      CALL MODJD(IYF,IMF,IDF,IHF,IMNF,TEND)
      WRITE(4)TSTART,TEND,NMS,NPASS,IOFF,DM2
C
C     READ ELEMENTS FILE
C
      CALL READEL(DEPOCH)
C
      DTIME=TSTART
      NQ=3
      NQ=2*NQ +1
      LUN=0
      TOL=1.0
C
C     NQ IS RANGE OF Q VALUES USED FOR ECC. FUNCTIONS (3 OR 4)
C     LUNAR AND SOLAR PERTS ARE OMITTED IF LUN=1
C      TESSERALS ARE EVALUATED AT MJD = DTIME
C     TOL IS ACCEPTANCE LEVEL IN MICRO RADIANS
C     LUNAR AND SOLAR PERTS ARE ALSO EVALUATED AT MJD
C
      IYON=0
      IF(DABS(DEPOCH-DTIME).LT.20.D0)GO TO 2
      WRITE(1,*)'REQUIRED EPOCH MORE THAN 20 DAYS FROM ELEM.EPOCH'
      WRITE(1,*)'NO TESSERALS EVALUATED. LUNAR + SOLAR PERTS ONLY'
      IYON=1
  2   CONTINUE
      CALL LONGP(DTIME,DEC,DES,DDI,DND,DKAP,DA,DDN,TH)
      EC=DEC
      ES=DES
      KAP=DKAP
      I=DDI
      ND=DND
      A=DA
      N=DDN
      W1=DW1
      N1=DN1
      E2=EC*EC +ES*ES
      E=SQRT(E2)
      DEGI=I/RAD
      W1=W1/360.0
      N1=N1/360.0
C
C     EVALUATE LUNAR AND SOLAR PERTURBATIONS
C
      DET=0.0
      DIT=0.0
      DMT=0.0
      DWT=0.0
      DNT=0.0
      IF(LUN.EQ.1) GO TO 201
      W=ATAN2(ES,EC)
      SJD=SNGL(DTIME)
      CALL LUNAR(SJD,E,I,W,ND,N,W1,N1,DE,DI,DM,DW,DN)
      ANOM=KAP -W
      WRITE(3,914) SJD,E,DE,I,DI,ANOM,DM,W,DW,ND,DN
 914  FORMAT(21H0LUNAR+SOLAR PERTS AT,F12.3//4H E :,3F14.8/4H I :,
     1 3F14.8/4H M :,3F14.8/4H W :,3F14.8/4H ND:,3F14.8///)
      DET=DE(1)
      DIT=DI(1)
      DMT=DM(1)
      DWT=DW(1)
      DNT=DN(1)
 201  CONTINUE
      WRITE(4)DET,DIT,DMT,DWT,DNT
C
C     COMPUTE TESSERAL PERTURBATIONS
C
C      SKIP IF EPOCH MORE THAN 20 DAYS FROM ELEMENTS EPOCH
C
      IF(IYON.EQ.0)GO TO 800
      IWU=1
      RD=0.0
      qWRITE(4)IWU,IWU,IWU,RD,RD,RD,RD,RD,RD,RD
      GO TO 200
 800  GTOL=TOL*(A -AE)/A
      WRITE(3,909) DTIME,A,E,DEGI,N,W1,N1
 909  FORMAT(1H0/20H0ELEMENTS AT EPOCH =,D21.11/4H A =,F15.10/
     1 4H E =,F15.10/4H I =,F15.10/4H N =,F15.10/4H W1=,F15.10/
     2 4H N1=,F15.10////)
      WRITE(3,910) TOL
 910  FORMAT(12H TOLERANCE =,F6.2,14H MICRO RADIANS////)
      WRITE(3,911)
 911  FORMAT(13H   L  M  P  Q,11H   W  M N-T,5X,4HDA/A,9X,2HDE,10X,2HDI,
     1 10X,2HDM,10X,2HDW,10X,2HDN,9X,5HPHASE//)
      GMA=GM/(N*A*A*A)
      AEA=AE/A
      ZE=1.0 -E2
      RE=SQRT(ZE)
      SI=SIN(I)
      CI=COS(I)
C
C     READ TESSERAL DATA FILE
C
      OPEN(UNIT=2,FILE='TESSER.DAT',STATUS='OLD')
      LP=0
 100  READ(2,920,END=200) L,M,C,S
 920  FORMAT(2I3,2E13.5)
      AAL=AEA**L
      AL=FLOAT(L)
      AM=FLOAT(M)
      K=1
      AA=C
      BB=S
      IF(MOD(L-M,2).EQ.0) GO TO 101
      AA=-S
      BB=C
 101  CONTINUE
C
C     CALCULATE ECCENTRICITY FUNCTIONS
C
      IF(L.EQ.LP) GO TO 104
      DO 102 IP=1,L+1
      P=IP -1
      DO 103 IQ=1,NQ
      Q=IQ -NQ/2  -1
      CALL ECCFUN(L,P,Q,E,GSS,GDSS)
      GS(IP,IQ)=GSS
      GDS(IP,IQ)=GDSS
 103  CONTINUE
 102  CONTINUE
C
C     CALCULATE INCLINATION FUNCTIONS
C     READ ECCENTRICITY FUNCTIONS FROM MEMORY
C
 104  CONTINUE
      DO 105 IP=1,L+1
      P=IP -1
      CALL INCFUN(L,M,P,DDI,F,FD)
      DO 106 IQ=1,NQ
      Q=IQ -NQ/2    -1
      G=GS(IP,IQ)
      GD=GDS(IP,IQ)
      AP=FLOAT(P)
      AQ=FLOAT(Q)
      LP=AL -2.0*AP
      LPQ=LP +AQ
      SD=LP*W1 +LPQ*N +AM*(N1 -THD)
      Z=AAL*GMA/SD
C
C     AMPLITUDES OF PERTURBATIONS
C
      R=Z*K*SQRT(AA*AA +BB*BB)
      RFG=R*F*G
      AMP(1)=2.0*LPQ*RFG
      AMP(2)=RFG*(AQ*ZE/E-LP*E*RE/(1.0+RE))
      AMP(3)=RFG*(LP*CI -AM)/(RE*SI)
      AMP(4)=R*F*(-3.0*N*LPQ*G/SD +(2.0*AL +2.0)*G -ZE*GD/E)
      AMP(5)=R*(RE*F*GD/E -CI*FD*G/(SI*RE))
      AMP(6)=R*FD*G/(RE*SI)
      IPRINT=0
      DO 107 II=1,6
      AMP(II)=AMP(II)*1.0E6
 107  CONTINUE
      DO 108 II=1,3
 108  IF(ABS(AMP(II)).GT.GTOL) IPRINT=1
      IF(ABS(AMP(6)*SI).GT.GTOL) IPRINT=1
      IF(ABS(AMP(6)*CI +AMP(5)).GT.GTOL/E) IPRINT=1
      IF(ABS(AMP(6)*CI +AMP(5) +AMP(4)).GT.GTOL) IPRINT=1
      IF(IPRINT.EQ.0) GO TO 109
      EPS=ATAN2(BB,AA)
      ILPQ=LPQ
      ILP=LP
      WRITE(3,912) L,M,P,Q,ILP,ILPQ,M,AMP,EPS
 912  FORMAT(1H ,4I3,2X,3I3,6F12.4,F15.6)
      WRITE(4)ILP,ILPQ,M,AMP,EPS
 109  CONTINUE
 106  CONTINUE
 105  CONTINUE
      LP=L
      GO TO 100
 200  CONTINUE
C
      CLOSE (UNIT=1)
      CLOSE (UNIT=2)
      CLOSE (UNIT=3)
      CLOSE (UNIT=4)
      CLOSE (UNIT=100)
      END
C
C      ************************************************************
C
      SUBROUTINE LUNAR(MJD,E,I,W,ND,N,W1,N1,DE,DI,DM,DW,DN)
C
C     COMPUTES LUNAR PERTURBATIONS OF 14-DAY PERIOD AT MOD.JUL.DATE MJD
C     ALSO SOLAR PERTS OF 6-MONTH PERIOD
C     INPUT ANGLES IN RADIANS:  N,W1,N1 IN REVS/DAY
C     OUTPUT PERTURBATIONS IN RADIANS
C
      REAL MJD,I,ND,N,N1,NDD,LI,LND,NM,MU,NDM,IM,LM,NB,NDH
      REAL D(3),CA(3),SA(3),CP(15),SP(15),Q(15),QD(15),
     1      DE(2),DI(2),DM(2),DW(2),DN(2)
      RAD=0.01745329252
      TD=MJD-44000.0
      EPS=23.441*RAD
      DO 200 K=1,2
      IF(K.EQ.2)GO TO 201
C
C      LUNAR DATA
C
      NDD=(164.552 -0.0529539*TD)*RAD
      LI=5.1454*RAD
      LND=(4.442 +13.2293505*TD)*RAD
      NM=0.036601
      MU=0.0121506
      SI=SIN(LI)
      CI=COS(LI)
      SE=SIN(EPS)
      CE=COS(EPS)
      SN=SIN(NDD)
      CN=COS(NDD)
      NDM=ATAN2(SI*SN, CI*SE +SI*CE*CN)
      ETA=ATAN2(SE*SN, SI*CE +CI*SE*CN)
      IM=ATAN2(SE*SN/SIN(ETA),CI*CE-SI*SE*CN)
      LM=LND +ETA
      GO TO 202
C
C      SOLAR DATA
C
 201      CONTINUE
      NM=0.273802E-2
      MU=0.999997
      NDM=0.0
      IM=EPS
      LM=(44.249+0.9856473*TD)*RAD
C
C      COMMON CALCULATION FOR LUNAR AND SOLAR PERTS.
C
 202      CONTINUE
      NB=N1/NM
      WB=W1/NM
      F=MU*NM/N
      PREC=(MJD -33281.923)*3.07364*RAD/(365.2422*240.0)
      NDH=ND +PREC -NDM
C
C     CALCULATE TRIG FUNCTIONS
C
      D(1)=2.0*(1.0 +WB -NB)
      D(2)=2.0*(1.0 -WB -NB)
      D(3)=2.0*(1.0 -NB)
      ARG=2.0*(LM +W -NDH)
      CA(1)=COS(ARG)
      SA(1)=SIN(ARG)
      ARG=2.0*(LM -W -NDH)
      CA(2)=COS(ARG)
      SA(2)=SIN(ARG)
      ARG=2.0*(LM -NDH)
      CA(3)=COS(ARG)
      SA(3)=SIN(ARG)
      CN=COS(NDH)
      SN=SIN(NDH)
      DO 100 IA=1,3
      DEN=D(IA)
      CB=CA(IA)
      SB=SA(IA)
      IC=5*IA -5
      DO 101 IB=1,5
      ID=IC +IB
      CP(ID)=CB/DEN
      SP(ID)=SB/DEN
      IF(IB.EQ.5) GO TO 101
      DEN=DEN +NB
      Z= CB*CN -SB*SN
      SB=SB*CN +CB*SN
      CB=Z
 101  CONTINUE
 100  CONTINUE
C
C     CALCULATE COEFFICIENTS
C
      S=SIN(I/2.0)
      G=COS(I/2.0)
      SM=SIN(IM/2.0)
      GM=COS(IM/2.0)
      X1=S**4
      FAC=G/S
      X2=X1*FAC
      X3=X2*FAC
      X4=X3*FAC
      X5=X4*FAC
      Y1=SM**4
      FAC=GM/SM
      Y2=Y1*FAC
      Y3=Y2*FAC
      Y4=Y3*FAC
      Y5=Y4*FAC
      Q(1)=X1*Y5
      Q(2)=-4.0*X2*Y4
      Q(3)=6.0*X3*Y3
      Q(4)=-4.0*X4*Y2
      Q(5)=X5*Y1
      Q(6)=X5*Y5
      Q(7)=4.0*X4*Y4
      Q(8)=Q(3)
      Q(9)=4.0*X2*Y2
      Q(10)=X1*Y1
      Q(11)=X3*Y5
      Q(12)=2.0*(X2 -X4)*Y4
      Q(13)=(X1 -4.0*X3 +X5)*Y3
      Q(14)=2.0*(X4 -X2)*Y2
      Q(15)=X3*Y1
      QD(1)=2.0*X2*Y5
      QD(2)=2.0*(X1 -3.0*X3)*Y4
      QD(3)=6.0*(X4 -X2)*Y3
      QD(4)=2.0*(3.0*X3 -X5)*Y2
      QD(5)=-2.0*X4*Y1
      QD(6)=-2.0*X4*Y5
      QD(7)=2.0*(X5 -3.0*X3)*Y4
      QD(8)=QD(3)
      QD(9)=2.0*(3.0*X3 -X1)*Y2
      QD(10)=2.0*X2*Y1
      QD(11)=(X4 -X2)*Y5
      QD(12)=(-X1 +6.0*X3 -X5)*Y4
      QD(13)=6.0*(X2 -X4)*Y3
      QD(14)=(X1 -6.0*X3 +X5)*Y2
      QD(15)=(X4 -X2)*Y1
C
C     SUMMATION OF PERIODIC TERMS
C
      SUM1=0.0
      SUM2=0.0
      SUM3=0.0
      SUM4=0.0
      SGN=1.0
      DO 102 IA=1,10
      IF(IA.GE.6) SGN=-1.0
      AM=MOD(IA-1,5) -2.0
      SUM1=SUM1 +Q(IA)*SP(IA)
      SUM2=SUM2 +Q(IA)*CP(IA)*SGN
      SUM3=SUM3 +Q(IA)*CP(IA)*AM
      SUM4=SUM4 +QD(IA)*SP(IA)
 102  CONTINUE
        E2=E*E
        ZA=1.875*E2
      ZB=1.5+2.25*E2
      R=ZA*SUM1
      RW=2.0*ZA*SUM2
      RN=ZA*SUM3
      RE=3.75*E*SUM1
      RI=ZA*SUM4
      SUM1=0.0
      SUM3=0.0
      SUM4=0.0
      DO 103 IA=11,15
      AM=MOD(IA-1,5)-2.0
      SUM1=SUM1+Q(IA)*SP(IA)
      SUM3=SUM3 +Q(IA)*CP(IA)*AM
      SUM4=SUM4 +QD(IA)*SP(IA)
 103  CONTINUE
      R=R +ZB*SUM1
      RN=RN +ZB*SUM3
      RE=RE +4.5*E*SUM1
      RI=RI +ZB*SUM4
C
C     CALCULATE PERTURBATIONS
C
      ZC=1.0 -E2
      ZD=SQRT(ZC)
      CI=COS(I)
      SI=SIN(I)
      DE(K)=-F*ZD*RW/E
      DW(K)=F*(-CI*RI/(ZD*SI) +ZD*RE/E)
      DI(K)=F*(CI*RW -RN)/(ZD*SI)
      DN(K)=F*RI/(ZD*SI)
      DM(K)=-F*(ZC*RE/E +4.0*R)
 200      CONTINUE
      RETURN
      END
C
C      ***********************************************************
C
      SUBROUTINE ECCFUN(L,P,Q,E,G,GD)
C
C     COMPUTES ECCENTRICITY FUNCTION AND DERIVATIVE
C
      INTEGER P,Q,PD,QD,H,R,P2,PL,RB
      MQ=IABS(Q)
      RE=SQRT(1.0 -E*E)
      GA=1.0 +RE
      BE=E/GA
      BE2=BE*BE
      BED=1.0/GA +BE2/RE
      DE=BE/RE
      PD=L -P
      QD=-Q
      IF(P.GT.L/2) GO TO 100
      PD=P
      QD=Q
 100  CONTINUE
      P2=-2*PD
      PL=2*PD -2*L
      GLPQ=GA*FLOAT(L +P2 +QD)/2.0
      BKK=1.0
      BK=1.0
      DO 101 KK=1,5
      K=KK -1
      H=K +QD
      N=K
      IF(QD.GT.0) GO TO 102
      H=K
      N=K -QD
 102  CONTINUE
      SP=BIN(PL,H)
      SQ=BIN(P2,N)
      SQD=0.0
      SPD=0.0
      RB=K +MQ
      IF(RB.EQ.0) GO TO 110
      SGN=1.0
      A=1.0
      DO 103 R=1,RB
      FR=FLOAT(R)
      A=A*GLPQ/FR
      IF(H.LT.R) GO TO 104
      SGN=-SGN
      B=BIN(PL,H-R)*SGN*A
      SP=SP +B
      SPD=SPD -B*FR*DE
 104  CONTINUE
      IF(N.LT.R) GO TO 103
      B=BIN(P2,N-R)*A
      SQ=SQ +B
      SQD=SQD -B*FR*DE
 103  CONTINUE
 110  CONTINUE
      IF(K.NE.0) GO TO 105
      C=SP*SQ
      D=SPD +SQD
      GO TO 101
 105  BK=BK*BE2
      BKK=BKK*BE
      IF(K.NE.1) BKK=BKK*BE
      C=C +SP*SQ*BK
      D=D +2.0*FLOAT(K)*SP*SQ*BED*BKK +(SPD*SQ +SP*SQD)*BK
 101  CONTINUE
      Y=1.0
      Z=1.0
      X=(1.0 +BE2)**L
      IF(MQ.EQ.0) GO TO 107
      DO 106 I=1,MQ
      X=-X
      Y=Y*BE
      IF(I.EQ.1) GO TO 106
      Z=Z*BE
 106  CONTINUE
 107  Z=Z*FLOAT(MQ)
      G=X*Y*C
      GD=X*BED*(2.0*FLOAT(L)*BE*Y/(1.0+BE*BE) +Z)*C +X*Y*D
      RETURN
      END
C
C      **********************************************************
C
      SUBROUTINE INCFUN(L,M,P,I,SF,SFD)
C
C  COMPUTES INCLINATION FUNCTION WITH SINGLE SUMMATION
C
      DOUBLE PRECISION F,FD,CI,SI,S,I,
     1  SGN,ZA,ZB,SX,C,CX,X,DBIN,DI
      INTEGER P,P2
      DI=I/2.D0
      CI=DCOS(DI)
      SI=DSIN(DI)
      P2=2*P
      LP2=2*L -P2
      LM=L -M
      KA=MAX0(0, LM-P2) +1
      KB=MIN0(LP2,LM) +1
      KK=(LM +1)/2
      SS=SQRT(4.0*FLOAT(L) +2.0)/(2.0**L)*(-1.0)**KK
      A1=RFACL(L+M)
      A2=RFACL(L-M)
      A3=RFACL(P)
      A4=RFACL(L-P)
C
C      USE LOGS TO AVOIB EXCEEDING MAX EXPONENT OF COMPUTER
C
      A=A1+A2-2.0*A3-2.0*A4-76.0
      A=10.0**A
      A=A*SS
      SGN=(-1.D0)**KA
      F=0.D0
      FD=0.D0
      DO 100 IK=KA,KB
      SGN=-SGN
      K=IK -1
      N=P2 +2*K -LM
      LN=2*L -N
      ZA=DBIN(LP2,K)
      ZB=DBIN(P2,LM-K)
      S=1.D0
      SX=1.D0
      C=1.D0
      CX=1.D0
      IF(N.EQ.0) GO TO 101
      NN=N -1
      IF(NN.NE.0) SX=SI**NN
      S=SX*SI
 101      IF(LN.EQ.0) GO TO 102
      LNN=LN -1
      IF(LNN.NE.0) CX=CI**LNN
      C=CX*CI
 102      CONTINUE
      X=0.5D0*(DBLE(FLOAT(N))*C*CI*SX -DBLE(FLOAT(LN))*S*SI*CX)
      F=F +SGN*ZA*ZB*S*C
      FD=FD +SGN*ZA*ZB*X
 100      CONTINUE
      SF=A*F
      SFD=A*FD
      RETURN
      END
C
      FUNCTION RFACL(N)
C     **************
C
C      RETURNS ALOG SQRT(N FACTORIAL)*1.0E-38 FOR N=0,......,97
C
      REAL F(98)
        DATA F/
     # -.38000000E+02, -.38000000E+02, -.37849480E+02, -.37610920E+02,
     # -.37309890E+02, -.36960410E+02, -.36571330E+02, -.36148780E+02,
     # -.35697240E+02, -.35220120E+02, -.34720120E+02, -.34199420E+02,
     # -.33659830E+02, -.33102860E+02, -.32529800E+02, -.31941750E+02,
     # -.31339690E+02, -.30724470E+02, -.30096830E+02, -.29457450E+02,
     # -.28806940E+02, -.28145830E+02, -.27474620E+02, -.26793750E+02,
     # -.26103650E+02, -.25404680E+02, -.24697190E+02, -.23981510E+02,
     # -.23257930E+02, -.22526730E+02, -.21788170E+02, -.21042490E+02,
     # -.20289910E+02, -.19530660E+02, -.18764920E+02, -.17992880E+02,
     # -.17214730E+02, -.16430630E+02, -.15640740E+02, -.14845210E+02,
     # -.14044180E+02, -.13237790E+02, -.12426160E+02, -.11609430E+02,
     # -.10787700E+02, -.99610940E+01, -.91297150E+01, -.82936660E+01,
     # -.74530450E+01, -.66079470E+01, -.57584620E+01, -.49046770E+01,
     # -.40466760E+01, -.31845380E+01, -.23183410E+01, -.14481600E+01,
     # -.57406560E+00,  .30387180E+00,  .11855860E+01,  .20710120E+01,
     #  .29600880E+01,  .38527520E+01,  .47489480E+01,  .56486190E+01,
     #  .65517080E+01,  .74581650E+01,  .83679370E+01,  .92809740E+01,
     #  .10197230E+02,  .11116650E+02,  .12039200E+02,  .12964830E+02,
     #  .13893500E+02,  .14825160E+02,  .15759780E+02,  .16697310E+02,
     #  .17637710E+02,  .18580960E+02,  .19527010E+02,  .20475820E+02,
     #  .21427360E+02,  .22381610E+02,  .23338510E+02,  .24298050E+02,
     #  .25260190E+02,  .26224900E+02,  .27192150E+02,  .28161910E+02,
     #  .29134150E+02,  .30108850E+02,  .31085970E+02,  .32065490E+02,
     #  .33047380E+02,  .34031620E+02,  .35018190E+02,  .36007050E+02,
     #  .36998180E+02,  .37991570E+02/
      IF(N.LT.0) GO TO 100
      IF(N.GT.97) GO TO 101
      M=N +1
      RFACL=F(M)
      RETURN
 100  WRITE(1,105) N
 105  FORMAT('NEGATIVE ARGUMENT OF ROOT-FACTORIAL',I4)
      STOP
 101  WRITE(1,106) N
106   FORMAT('ARGUMENT OF ROOT-FACTORIAL TOO LARGE',I4)
      STOP
      END
C
C      *********************************************************
C
      FUNCTION DBIN(N,M)
C
C      COMPUTES BINOMIAL COEFFICIENT
C
      DOUBLE PRECISION DBIN,AN,AM
      IF(M) 103,102,99
 99      IF(N.EQ.0) GO TO 101
      AN=DBLE(FLOAT(N))
      AM=DBLE(FLOAT(M))
      DBIN=1.D0
      DO 100 I=1,M
      DBIN=DBIN*AN/AM
      AN=AN -1.D0
      AM=AM -1.D0
 100      CONTINUE
      RETURN
 101      DBIN=0.D0
      RETURN
 102      DBIN=1.D0
      RETURN
 103      DBIN=0.D0
      WRITE(1,104) N,M
 104  FORMAT('NEGATIVE LOWER ARGUMENT OF BINOMIAL',2I5)
      RETURN
      END
C
C      ********************************************************
C
      FUNCTION BIN(N,M)
C
C     COMPUTES BINOMIAL COEFFICIENT
C
      IF(M) 103,102,90
 90   IF(N.EQ.0) GO TO 101
      AN=FLOAT(N)
      AM=FLOAT(M)
      BIN=1.0
      DO 100 I=1,M
      BIN=BIN*AN/AM
      AN=AN -1.0
      AM=AM -1.0
 100  CONTINUE
      RETURN
 101  BIN=0.0
      RETURN
 102  BIN=1.0
      RETURN
 103  BIN=0.0
      WRITE(6,900) N,M,BIN
 900  FORMAT(1H0,'NEGATIVE LOWER ARGUMENT OF BINOMIAL'/
     13H N=,I4,3H M=,I4,5H BIN=,F3.1//)
      RETURN
      END
C
C      *******************************************************
C
      FUNCTION FACT(N)
C
C     COMPUTES FACTORIAL
C
      IF(N) 102,101,90
 90   FACT=1.0
      DO 100 I=1,N
      FACT=FACT*FLOAT(I)
 100  CONTINUE
      RETURN
 101  FACT=1.0
      RETURN
 102  FACT=0.0
      WRITE(6,900) N,FACT
 900  FORMAT(1H0,'NEGATIVE ARGUMENT OF FACTORIAL'/
     13H N=,I4,6H FACT=,F3.1//)
      RETURN
      END
C
C      ******************************************************
C
      SUBROUTINE LONGP(MJD,EC,ES,I,ND,KAP,A,N,TH)
C
C      CALCULATES SECULAR & LONG PERIOD PERTURBATIONS
C
      REAL M2,M3,B(5,3)
      DOUBLE PRECISION EPOCH,MJD,DT,M1,DM,W0,W1,N0,N1,I0,I1,DTT,
     1  E0,E1,M0,RAD,I,ND,EC,KAP,ES,W,E,M,CW,SW,ZA,N,RE,ZB,SI,P,A,GM
       DOUBLE PRECISION J,AJ,DELA
      COMMON/CBA/EPOCH,W0,W1,N0,N1,I0,I1,E0,E1,M0,M1,M2,M3,
     1    ETA,B,DELA
      COMMON/CBC/RAD,GM,J,DT
      DT=MJD-EPOCH
      T=DT
      W=(W0+W1*DT)*RAD
      ND=(N0+N1*DT)*RAD
      I=(I0+I1*DT)*RAD
      E=E0+E1*DT
        DM=M0+M1*DT+DBLE(M2)*DT*DT+DBLE(M3)*DT*DT*DT
      DM=DM-IDINT(DM)
      M=DM*360.D0*RAD
      CW=DCOS(W)
      SW=DSIN(W)
        C2W=SNGL(CW*CW-SW*SW)
        S2W=SNGL(2.0D0*SW*CW)
        C3W=C2W*SNGL(CW)-S2W*SNGL(SW)
        S3W=S2W*SNGL(CW)+C2W*SNGL(SW)
      EC=E*CW
      ES=E*SW
      KAP=M+W
          SCW=SNGL(CW)
         SSW=SNGL(SW)
      DEC=B(1,1)*SCW+B(1,2)*S2W+B(1,3)*C3W
      DND=(B(2,1)*SCW+B(2,2)*S2W+B(2,3)*C3W)*RAD
      DI=(B(3,1)*SSW+B(3,2)*C2W+B(3,3)*S3W)*RAD
      DES=B(4,1)*SSW+B(4,2)*C2W+B(4,3)*S3W
      DES=DES+ETA
      DKAP=B(5,1)*SCW+B(5,2)*S2W+B(5,3)*C3W
C
       EC=EC+DBLE(DEC)
       ND=ND+DBLE(DND)
       I=I+DBLE(DI)
       KAP=KAP+DBLE(DKAP)
       ES=ES+DBLE(DES)
C
C      CALC P & A  SEE SAO SR200,PAGE 101
C
  10    CONTINUE
      ZA=1.0D0-EC*EC-ES*ES
      RE=DSQRT(ZA)
       N=M1+2.0D0*DBLE(M2)*DT+3.0D0*DBLE(M3)*DT*DT
      ZB=(GM/(N*N))**(0.33333333333D0)
      SI=DSIN(I)
      P=ZB*ZA
      A=ZB*(1.0D0+DBLE(J)*RE*(1.5D0*SI*SI-1.0D0)/(3.0D0*P*P))
C
      TS=T*86400.0
      SEPOCH=EPOCH
      TH=1.7466477+0.01720218*(SEPOCH-33282.0)+7.292115E-5*TS
      TH =AMOD(TH,6.2831853)
      RETURN
      END
C
C      *****************************************************
C
      SUBROUTINE READEL(DEPOCH)
       IMPLICIT DOUBLE PRECISION (D)
       REAL B(5,3),M2,M3,NMS(2)
      COMMON/CBA/DEP,DW0,DW1,DN0,DN1,DI0,DI1,DE0,DE1,DM0,DM1,
     1  M2,M3,ETA,B,DELA
      COMMON/CBB/DDW1,DDN1,NMS
C
      OPEN(UNIT=100,FILE='SAOELM.DAT',STATUS='OLD')
C
  12  READ(100,919,END=20)RN1,RN2,EPOCH,FRACT
 919  FORMAT(2A4,1X,F6.0,F9.6,F7.2)
      IF(RN1.NE.NMS(1))GO TO 15
      IF(RN2.NE.NMS(2))GO TO 15
        DEPOCH=DBLE(EPOCH)+DBLE(FRACT)
      TIME=DTIME
      DEP=DEPOCH
      EPOCH=EPOCH+FRACT
      READ(100,*)DW0,DW1
      READ(100,*)DN0,DN1
      READ(100,*)DI0,DI1
      READ(100,*)DE0,DE1
      READ(100,*)DM0,DM1,M2,M3
      READ(100,*)((B(K,JJ),JJ=1,3),K=1,5)
      READ(100,*)ETA
      DDW1=DW1
      DDN1=DN1
C
      CLOSE (UNIT=100)
      RETURN
   15   DO 16 I=1,11
      READ(100,916)
 916  FORMAT(1X)
  16  CONTINUE
      GO TO 12
20    WRITE(1,*)'SATELLITE ELEMENTS NOT FOUND'
      CLOSE (UNIT=100)
      RETURN
      END
C
C      *************************************************************
C
      SUBROUTINE MODJD(I,J,K,IH,MIN,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) +DBLE(FLOAT(IH*60 +MIN))/1440.D0
      RETURN
      END
