* DEEP SPACE 31 OCT 80

	SUBROUTINE DEEPD

	IMPLICIT DOUBLE PRECISION (A-H,O-Z)

	COMMON/E1/XMO,XNODEO,OMEGAO,EO,XINCL,XNO,XNDT2O,XNDD6O,BSTAR,X,Y,Z,XDOT,YDOT,ZDOT,EPOCH,DS50
	COMMON/C1/CK2,CK4,E6A,QOMS2T,S,TOTHRD,XJ3,XKE,XKMPER,XMNPDA,AE
	COMMON/C2/DE2RA,PI,PIO2,TWOPI,X3PIO2

	DATA ZNS, C1SS, ZES/1.19459d-5, 2.9864797d-6, 0.01675d0/
	DATA ZNL, C1L, ZEL/1.5835218d-4, 4.7968065d-7, 0.05490d0/
	DATA ZCOSIS, ZSINIS, ZSINGS/0.91744867d0, 0.39785416d0, -0.98088458d0/
	DATA ZCOSGS, ZCOSHS, ZSINHS/0.1945905d0, 1d0, 0d0/
	DATA Q22,Q31,Q33/1.7891679d-6,2.1460748d-6,2.2123015d-7/
	DATA G22,G32/5.7686396d0,0.95240898d0/
	DATA G44,G52/1.8014998d0,1.0508330d0/
	DATA G54/4.4108898d0/
	DATA ROOT22,ROOT32/1.7891679d-6,3.7393792d-7/
	DATA ROOT44,ROOT52/7.3636953d-9,1.1428639d-7/
	DATA ROOT54/2.1765803d-9/
	DATA THDT/4.3752691d-3/

	SAVE PREEP,ZCOSIL, ZSINIL, ZSINHL, ZCOSHL ZCOSGL, ZSINGL
	SAVE THGR, XNQ, XQNCL, OMEGAQ, ZMOL, ZMOS, SAVTSN
	SAVE SGH2, SGH3, SGH4, XGH2, XGH3, XGH4, SSL, SSE, SSI, SSG, SSH, SE2, SI2, SL2, SH2
	SAVE SE3, SI3, SL3, SH3, SL4
	SAVE IRESFL, ISYNFL, XLAMO, DEL1, DEL2, DEL3, FASX2, FASX4, FASX6, XFACT
	SAVE D2201, D2211, D3210, D3222, D4410, D4422, D5220, D5232, D5421, D5433
	SAVE XLI, XNI, ATIME, STEPP, STEPN, STEP2
	SAVE SGHS, SHS, SGHL, SHL, PE, PINC, PL

	! ENTRANCE FOR DEEP SPACE INITIALIZATION

	ENTRY DPINITD(EQSQ,SINIQ,COSIQ,RTEQSQ,AO,COSQ2,SINOMO,COSOMO,BSQ,XLLDOT,OMGDT,XNODOT,XNODP)

	THGR=THETAGD(EPOCH)
	EQ = EO
	XNQ = XNODP
	AQNV = 1d0/AO
	XQNCL = XINCL
	XMAO=XMO
	XPIDOT=OMGDT+XNODOT
	SINQ = dsin(XNODEO)
	COSQ = dcos(XNODEO)
	OMEGAQ = OMEGAO

	! INITIALIZE LUNAR SOLAR TERMS

   5	DAY=DS50+18261.5d0
	IF (DAY .NE. PREEP) THEN
	    PREEP = DAY
	    XNODCE=4.5236020d0-9.2422029d-4*DAY
	    STEM=dsin (XNODCE)
	    CTEM=dcos (XNODCE)
	    ZCOSIL=0.91375164d0-0.03568096d0*CTEM
	    ZSINIL=dsqrt (1d0-ZCOSIL*ZCOSIL)
	    ZSINHL=0.089683511d0*STEM/ZSINIL
	    ZCOSHL=dsqrt (1d0-ZSINHL*ZSINHL)
	    C=4.7199672d0+0.22997150d0*DAY
	    GAM=5.8351514d0+0.0019443680d0*DAY
	    ZMOL = dmod(dmod(C-GAM,TWOPI)+TWOPI,TWOPI)
	    ZX= 0.39785416d0*STEM/ZSINIL
	    ZY= ZCOSHL*CTEM+0.91744867d0*ZSINHL*STEM
	    ZX=dmod(datan2(ZX,ZY)+TWOPI,TWOPI)
	    ZX=GAM+ZX-XNODCE
	    ZCOSGL=dcos (ZX)
	    ZSINGL=dsin (ZX)
	    ZMOS=6.2565837d0+0.017201977d0*DAY
	    ZMOS=dmod(dmod(ZMOS,TWOPI)+TWOPI,TWOPI)
	END IF

	! DO SOLAR TERMS

	LS = 0
	SAVTSN=1.d20
	ZCOSG=ZCOSGS
	ZSING=ZSINGS
	ZCOSI=ZCOSIS
	ZSINI=ZSINIS
	ZCOSH=COSQ
	ZSINH=SINQ
	CC=C1SS
	ZN=ZNS
	ZE=ZES
	ZMO=ZMOS	! NOT USED?
	XNOI=1d0/XNQ
	ASSIGN 30 TO LS

  20	A1=ZCOSG*ZCOSH+ZSING*ZCOSI*ZSINH
	A3=-ZSING*ZCOSH+ZCOSG*ZCOSI*ZSINH
	A7=-ZCOSG*ZSINH+ZSING*ZCOSI*ZCOSH
	A8=ZSING*ZSINI
	A9=ZSING*ZSINH+ZCOSG*ZCOSI*ZCOSH
	A10=ZCOSG*ZSINI
	A2= COSIQ*A7+ SINIQ*A8
	A4= COSIQ*A9+ SINIQ*A10
	A5=- SINIQ*A7+ COSIQ*A8
	A6=- SINIQ*A9+ COSIQ*A10

	X1=A1*COSOMO+A2*SINOMO
	X2=A3*COSOMO+A4*SINOMO
	X3=-A1*SINOMO+A2*COSOMO

	X4=-A3*SINOMO+A4*COSOMO
	X5=A5*SINOMO
	X6=A6*SINOMO
	X7=A5*COSOMO
	X8=A6*COSOMO

	Z31=12d0*X1*X1-3d0*X3*X3
	Z32=24d0*X1*X2-6d0*X3*X4
	Z33=12d0*X2*X2-3d0*X4*X4
	Z1=3d0*(A1*A1+A2*A2)+Z31*EQSQ
	Z2=6d0*(A1*A3+A2*A4)+Z32*EQSQ
	Z3=3d0*(A3*A3+A4*A4)+Z33*EQSQ
	Z11=-6d0*A1*A5+EQSQ *(-24d0*X1*X7-6d0*X3*X5)
	Z12=-6d0*(A1*A6+A3*A5)+EQSQ *(-24d0*(X2*X7+X1*X8)-6d0*(X3*X6+X4*X5))
	Z13=-6d0*A3*A6+EQSQ *(-24d0*X2*X8-6d0*X4*X6)
	Z21=6d0*A2*A5+EQSQ *(24d0*X1*X5-6d0*X3*X7)
	Z22=6d0*(A4*A5+A2*A6)+EQSQ *(24d0*(X2*X5+X1*X6)-6d0*(X4*X7+X3*X8))
	Z23=6d0*A4*A6+EQSQ *(24d0*X2*X6-6d0*X4*X8)
	Z1=Z1+Z1+BSQ*Z31
	Z2=Z2+Z2+BSQ*Z32
	Z3=Z3+Z3+BSQ*Z33
	S3=CC*XNOI
	S2=-0.5d0*S3/RTEQSQ
	S4=S3*RTEQSQ
	S1=-15d0*EQ*S4
	S5=X1*X3+X2*X4
	S6=X2*X3+X1*X4
	S7=X2*X4-X1*X3
	SE=S1*ZN*S5
	SI=S2*ZN*(Z11+Z13)
	SL=-ZN*S3*(Z1+Z3-14d0-6d0*EQSQ)
	SGH=S4*ZN*(Z31+Z33-6d0)
	SH=-ZN*S2*(Z21+Z23)
	IF(XQNCL.LT.5.2359877d-2) SH=0d0
	EE2=2d0*S1*S6
	E3=2d0*S1*S7
	XI2=2d0*S2*Z12
	XI3=2d0*S2*(Z13-Z11)
	XL2=-2d0*S3*Z2
	XL3=-2d0*S3*(Z3-Z1)
	XL4=-2d0*S3*(-21d0-9d0*EQSQ)*ZE
	XGH2=2d0*S4*Z32
	XGH3=2d0*S4*(Z33-Z31)
	XGH4=-18d0*S4*ZE
	XH2=-2d0*S2*Z22
	XH3=-2d0*S2*(Z23-Z21)
	GO TO LS

	! DO LUNAR TERMS

  30	SSE = SE
	SSI=SI
	SSL=SL
	SSH=SH/SINIQ
	SSG=SGH-COSIQ*SSH
	SE2=EE2
	SI2=XI2
	SL2=XL2
	SGH2=XGH2
	SH2=XH2
	SE3=E3
	SI3=XI3
	SL3=XL3
	SGH3=XGH3
	SH3=XH3
	SL4=XL4
	SGH4=XGH4
	LS=1
	ZCOSG=ZCOSGL
	ZSING=ZSINGL
	ZCOSI=ZCOSIL
	ZSINI=ZSINIL
	ZCOSH=ZCOSHL*COSQ+ZSINHL*SINQ
	ZSINH=SINQ*ZCOSHL-COSQ*ZSINHL
	ZN=ZNL
	CC=C1L
	ZE=ZEL
	ZMO=ZMOL	! NOT USED ?
	ASSIGN 40 TO LS
	GO TO 20

  40	SSE = SSE+SE
	SSI=SSI+SI
	SSL=SSL+SL
	SSG=SSG+SGH-COSIQ/SINIQ*SH
	SSH=SSH+SH/SINIQ

	! GEOPOTENTIAL RESONANCE INITIALIZATION FOR 12 HOUR ORBITS

	IRESFL=0
	ISYNFL=0
	IF (XNQ .GE. 0.0052359877d0 .OR. XNQ .LE. 0.0034906585d0) THEN
	    IF (XNQ .LT. 8.26d-3 .OR. XNQ .GT. 9.24d-3) RETURN
	    IF (EQ .LT. 0.5d0) RETURN
	    IRESFL =1
	    EOC=EQ*EQSQ
	    G201=-0.306d0-(EQ-0.64d0)*0.440d0
	    IF (EQ .LE. 0.65d0) THEN
		G211=3.616d0-13.247d0*EQ+16.290d0*EQSQ
		G310=-19.302d0+117.390d0*EQ-228.419d0*EQSQ+156.591d0*EOC
		G322=-18.9068d0+109.7927d0*EQ-214.6334d0*EQSQ+146.5816d0*EOC
		G410=-41.122d0+242.694d0*EQ-471.094d0*EQSQ+313.953d0*EOC
		G422=-146.407d0+841.880d0*EQ-1629.014d0*EQSQ+1083.435d0*EOC
		G520=-532.114d0+3017.977d0*EQ-5740d0*EQSQ+3708.276d0*EOC
	    ELSE
		G211=-72.099d0+331.819d0*EQ-508.738d0*EQSQ+266.724d0*EOC
		G310=-346.844d0+1582.851d0*EQ-2415.925d0*EQSQ+1246.113d0*EOC
		G322=-342.585d0+1554.908d0*EQ-2366.899d0*EQSQ+1215.972d0*EOC
		G410=-1052.797d0+4758.686d0*EQ-7193.992d0*EQSQ+3651.957d0*EOC
		G422=-3581.69d0+16178.11d0*EQ-24462.77d0*EQSQ+12422.52d0*EOC
		IF (EQ .LE. 0.715d0) THEN
		    G520=1464.74d0-4664.75d0*EQ+3763.64d0*EQSQ
		ELSE
		    G520=-5149.66d0+29936.92d0*EQ-54087.36d0*EQSQ+31324.56d0*EOC
		END IF
	    END IF
	    IF (EQ .LT. 0.7d0) THEN
		G533=-919.2277d0+4988.61d0*EQ-9064.77d0*EQSQ+5542.21d0*EOC
		G521 = -822.71072d0+4568.6173d0*EQ-8491.4146d0*EQSQ+5337.524d0*EOC
		G532 = -853.666d0+4690.25d0*EQ-8624.77d0*EQSQ+5341.4d0*EOC
	    ELSE
		G533=-37995.78d0+161616.52d0*EQ-229838.2d0*EQSQ+109377.94d0*EOC
		G521 = -51752.104d0+218913.95d0*EQ-309468.16d0*EQSQ+146349.42d0*EOC
		G532 = -40023.88d0+170470.89d0*EQ-242699.48d0*EQSQ+115605.82d0*EOC
	    END IF
	    SINI2=SINIQ*SINIQ
	    F220=0.75d0*(1d0+2d0*COSIQ+COSQ2)
	    F221=1.5d0*SINI2
	    F321=1.875d0*SINIQ*(1d0-2d0*COSIQ-3d0*COSQ2)
	    F322=-1.875d0*SINIQ*(1d0+2d0*COSIQ-3d0*COSQ2)
	    F441=35d0*SINI2*F220
	    F442=39.3750d0*SINI2*SINI2
	    F522=9.84375d0*SINIQ*(SINI2*(1d0-2d0*COSIQ-5d0*COSQ2)+0.33333333d0*(-2d0+4d0*COSIQ+6d0*COSQ2))
	    F523 = SINIQ*(4.92187512d0*SINI2*(-2d0-4d0*COSIQ+10d0*COSQ2)+6.56250012d0*(1d0+2d0*COSIQ-3d0*COSQ2))
	    F542 = 29.53125d0*SINIQ*(2d0-8d0*COSIQ+COSQ2*(-12d0+8d0*COSIQ+10d0*COSQ2))
	    F543=29.53125d0*SINIQ*(-2d0-8d0*COSIQ+COSQ2*(12d0+8d0*COSIQ-10d0*COSQ2))
	    XNO2=XNQ*XNQ
	    AINV2=AQNV*AQNV
	    TEMP1 = 3d0*XNO2*AINV2
	    TEMP = TEMP1*ROOT22
	    D2201 = TEMP*F220*G201
	    D2211 = TEMP*F221*G211
	    TEMP1 = TEMP1*AQNV
	    TEMP = TEMP1*ROOT32
	    D3210 = TEMP*F321*G310
	    D3222 = TEMP*F322*G322
	    TEMP1 = TEMP1*AQNV
	    TEMP = 2d0*TEMP1*ROOT44
	    D4410 = TEMP*F441*G410
	    D4422 = TEMP*F442*G422
	    TEMP1 = TEMP1*AQNV
	    TEMP = TEMP1*ROOT52
	    D5220 = TEMP*F522*G520
	    D5232 = TEMP*F523*G532
	    TEMP = 2d0*TEMP1*ROOT54
	    D5421 = TEMP*F542*G521
	    D5433 = TEMP*F543*G533
	    XLAMO = XMAO+XNODEO+XNODEO-THGR-THGR
	    BFACT = XLLDOT+XNODOT+XNODOT-THDT-THDT
	    BFACT=BFACT+SSL+SSH+SSH
	ELSE			! SYNCHRONOUS RESONANCE TERMS INITIALIZATION
	    IRESFL=1
	    ISYNFL=1
	    G200=1d0+EQSQ*(-2.5d0+0.8125d0*EQSQ)
	    G310=1d0+2d0*EQSQ
	    G300=1d0+EQSQ*(-6d0+6.60937d0*EQSQ)
	    F220=0.75d0*(1d0+COSIQ)*(1d0+COSIQ)
	    F311=0.9375*SINIQ*SINIQ*(1d0+3d0*COSIQ)-0.75*(1d0+COSIQ)
	    F330=1d0+COSIQ
	    F330=1.875d0*F330*F330*F330
	    DEL1=3d0*XNQ*XNQ*AQNV*AQNV
	    DEL2=2d0*DEL1*F220*G200*Q22
	    DEL3=3d0*DEL1*F330*G300*Q33*AQNV
	    DEL1=DEL1*F311*G310*Q31*AQNV
	    FASX2=0.13130908d0
	    FASX4=2.8843198d0
	    FASX6=0.37448087d0
	    XLAMO=XMAO+XNODEO+OMEGAO-THGR
	    BFACT = XLLDOT+XPIDOT-THDT
	    BFACT=BFACT+SSL+SSG+SSH
	END IF

	XFACT=BFACT-XNQ

	! INITIALIZE INTEGRATOR

	XLI=XLAMO
	XNI=XNQ
	ATIME=0.d0
	STEPP=720.d0
	STEPN=-720.d0
	STEP2 = 259200.d0

	RETURN

	! ENTRANCE FOR DEEP SPACE SECULAR EFFECTS

	ENTRY DPSECD(XLL,OMGASM,XNODES,EM,XINC,XN,OMGDT,T)
	XLL=XLL+SSL*T
	OMGASM=OMGASM+SSG*T
	XNODES=XNODES+SSH*T
	EM=EO+SSE*T
	XINC=XINCL+SSI*T
	IF (XINC .LT. 0d0) THEN
	    XINC = -XINC
	    XNODES = XNODES + PI
	    OMGASM = OMGASM - PI
	END IF
	IF(IRESFL .EQ. 0) RETURN
 100	IF (ATIME.EQ.0d0) GO TO 170
	IF(T.GE.0d0.AND.ATIME.LT.0d0) GO TO 170
	IF(T.LT.0d0.AND.ATIME.GE.0d0) GO TO 170
 105	IF(dabs(T).GE.dabs(ATIME)) GO TO 120
	DELT=STEPP
	IF (T.GE.0d0) DELT = STEPN
 110	ASSIGN 100 TO IRET
	GO TO 160
 120	DELT=STEPN
	IF (T.GT.0d0) DELT = STEPP
 125	IF (dabs(T-ATIME).LT.STEPP) GO TO 130
	ASSIGN 125 TO IRET
	GO TO 160
 130	FT = T-ATIME
	ASSIGN 140 TO IRETN
	GO TO 150
 140	XN = XNI+XNDOT*FT+XNDDT*FT*FT*0.5d0
	XL = XLI+XLDOT*FT+XNDOT*FT*FT*0.5d0
	TEMP = -XNODES+THGR+T*THDT
	XLL = XL-OMGASM+TEMP
	IF (ISYNFL.EQ.0) XLL = XL+TEMP+TEMP
	RETURN

	! DOT TERMS CALCULATED
 
 150	IF (ISYNFL .NE. 0) THEN
	    XNDOT=DEL1*dsin (XLI-FASX2)+DEL2*dsin (2d0*(XLI-FASX4))+DEL3*dsin (3d0*(XLI-FASX6))
	    XNDDT = DEL1*dcos(XLI-FASX2)+2d0*DEL2*dcos(2d0*(XLI-FASX4))+3d0*DEL3*dcos(3d0*(XLI-FASX6))
	ELSE
	    XOMI = OMEGAQ+OMGDT*ATIME
	    X2OMI = XOMI+XOMI
	    X2LI = XLI+XLI
	    XNDOT =  D2201*dsin(X2OMI+XLI-G22)+D2211*dsin(XLI-G22)+D3210*dsin(XOMI+XLI-G32)
     &		+D3222*dsin(-XOMI+XLI-G32)+D4410*dsin(X2OMI+X2LI-G44)+D4422*dsin(X2LI-G44)
     &		+D5220*dsin(XOMI+XLI-G52)+D5232*dsin(-XOMI+XLI-G52)+D5421*dsin(XOMI+X2LI-G54)
     &		+D5433*dsin(-XOMI+X2LI-G54)
	    XNDDT =  D2201*dcos(X2OMI+XLI-G22)+D2211*dcos(XLI-G22)+D3210*dcos(XOMI+XLI-G32)
     &		+D3222*dcos(-XOMI+XLI-G32)+D5220*dcos(XOMI+XLI-G52)+D5232*dcos(-XOMI+XLI-G52)
     &		+2.*(D4410*dcos(X2OMI+X2LI-G44)+D4422*COS(X2LI-G44)+D5421*dcos(XOMI+X2LI-G54)
     &		+D5433*dcos(-XOMI+X2LI-G54))
	END IF
	XLDOT=XNI+XFACT
	XNDDT = XNDDT*XLDOT
	GO TO IRETN

	! INTEGRATOR

  160	ASSIGN 165 TO IRETN
	GO TO 150
  165	XLI = XLI+XLDOT*DELT+XNDOT*STEP2
	XNI = XNI+XNDOT*DELT+XNDDT*STEP2
	ATIME=ATIME+DELT
	GO TO IRET

	! EPOCH RESTART

  170	IF (T .LT. 0d0) THEN
	    DELT=STEPN
	ELSE
	    DELT = STEPP
	END IF
	ATIME = 0d0
	XNI=XNQ
	XLI=XLAMO
	GO TO 125

	! ENTRANCES FOR LUNAR-SOLAR PERIODICS

	ENTRY DPPERD(SINIQ,COSIQ,EM,XINC,OMGASM,XNODES,XLL,T)
	SINIS = dsin(XINC)
	COSIS = dcos(XINC)
	IF (dabs(SAVTSN-T) .GE. 30d0) THEN
	    SAVTSN=T
	    ZM=ZMOS+ZNS*T
	    ZF=ZM+2d0*ZES*dsin (ZM)
	    SINZF=dsin (ZF)
	    F2= 0.5d0*SINZF*SINZF-0.25d0
	    F3=-0.5d0*SINZF*dcos (ZF)
	    SES=SE2*F2+SE3*F3
	    SIS=SI2*F2+SI3*F3
	    SLS=SL2*F2+SL3*F3+SL4*SINZF
	    SGHS=SGH2*F2+SGH3*F3+SGH4*SINZF
	    SHS=SH2*F2+SH3*F3
	    ZM=ZMOL+ZNL*T
	    ZF=ZM+2d0*ZEL*dsin (ZM)
	    SINZF=dsin (ZF)
	    F2= 0.5d0*SINZF*SINZF-0.25d0
	    F3=-0.5d0*SINZF*dcos (ZF)
	    SEL=EE2*F2+E3*F3
	    SIL=XI2*F2+XI3*F3
	    SLL=XL2*F2+XL3*F3+XL4*SINZF
	    SGHL=XGH2*F2+XGH3*F3+XGH4*SINZF
	    SHL=XH2*F2+XH3*F3
	    PE=SES+SEL
	    PINC=SIS+SIL
	    PL=SLS+SLL
	END IF
	PGH=SGHS+SGHL
	PH=SHS+SHL
	XINC = XINC+PINC
	EM = EM+PE
	IF (XQNCL .GE. 0.2d0) THEN	! APPLY PERIODICS DIRECTLY
	    PH=PH/SINIQ
	    PGH=PGH-COSIQ*PH
	    OMGASM=OMGASM+PGH
	    XNODES=XNODES+PH
	    XLL = XLL+PL
	ELSE				! APPLY PERIODICS WITH LYDDANE MODIFICATION
	    SINOK=dsin(XNODES)
	    COSOK=dcos(XNODES)
	    ALFDP=SINIS*SINOK
	    BETDP=SINIS*COSOK
	    DALF=PH*COSOK+PINC*COSIS*SINOK
	    DBET=-PH*SINOK+PINC*COSIS*COSOK
	    ALFDP=ALFDP+DALF
	    BETDP=BETDP+DBET
	    XLS = XLL+OMGASM+COSIS*XNODES
	    DLS=PL+PGH-PINC*XNODES*SINIS
	    XLS=XLS+DLS
	    XNODES=dmod(datan2(ALFDP,BETDP)+TWOPI,TWOPI)
	    XLL = XLL+PL
	    OMGASM = XLS-XLL-COS(XINC)*XNODES
	END IF

	RETURN
	END

	DOUBLE PRECISION FUNCTION THETAGD(EP)
	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
	COMMON /E1/XMO,XNODEO,OMEGAO,EO,XINCL,XNO,XNDT2O,XNDD6O,BSTAR,X,Y,Z,XDOT,YDOT,ZDOT,EPOCH,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
	THETAGD=THETA-TEMP*TWOPI
	IF(THETAGD.LT.0.d0) THETAGD=THETAGD+TWOPI

	RETURN
	END