* SGP4 3 NOV 80

	SUBROUTINE SGP4D(IFLAG,TSINCE)

	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

	SAVE COSIO, X3THM1, XNODP, AODP
	SAVE ISIMP, ETA, C1, SINIO, X1MTH2, C4, C5, XMDOT, OMGDOT, XNODOT, OMGCOF
	SAVE XMCOF, XNODCF, T2COF, XLCOF, AYCOF, DELMO, SINMO, X7THM1
	SAVE D2, D3, D4, T3COF, T4COF, T5COF

	IF (IFLAG .NE. 0) THEN

	    ! Recover original mean motion (XNODP) and semimajor axis (AODP) from input elements

	    A1	   = (XKE/XNO)**TOTHRD
	    COSIO  = dcos(XINCL)
	    THETA2 = COSIO*COSIO
	    X3THM1 = 3d0*THETA2-1d0
	    EOSQ   = EO*EO
	    BETAO2 = 1d0-EOSQ
	    BETAO  = dsqrt(BETAO2)
	    DEL1   = 1.5d0*CK2*X3THM1/(A1*A1*BETAO*BETAO2)
	    AO	   = A1*(1d0-DEL1*(0.5d0*TOTHRD+DEL1*(1d0+134d0/81d0*DEL1)))
	    DELO   = 1.5d0*CK2*X3THM1/(AO*AO*BETAO*BETAO2)
	    XNODP  = XNO/(1d0+DELO)
	    AODP   = AO/(1d0-DELO)

	    ! Initialization
	    ! For perigee less than 220 kilometers, the isimp flag is set and
	    ! the equations are truncated to linear variation in sqrt a and
	    ! quadratic variation in mean anomaly. also, the c3 term, the
	    ! delta omega term, and the delta m term are dropped.

	    ISIMP = 0
	    IF ((AODP*(1d0-EO)/AE) .LT. (220d0/XKMPER+AE)) ISIMP=1

	    ! For perigee below 156 km, the values of S and QOMS2T are altered

	    S4		= S
	    QOMS24	= QOMS2T
	    PERIGE	= (AODP*(1d0-EO)-AE)*XKMPER
	    IF (PERIGE .LT. 156d0) THEN
		S4	= PERIGE-78d0
		IF(PERIGE .LE. 98d0) S4=20d0
		QOMS24	= ((120d0-S4)*AE/XKMPER)**4d0
		S4	= S4/XKMPER+AE
	    END IF
 	    PINVSQ	= 1d0/(AODP*AODP*BETAO2*BETAO2)
	    TSI		= 1d0/(AODP-S4)
	    ETA		= AODP*EO*TSI
	    ETASQ	= ETA*ETA
	    EETA	= EO*ETA
	    PSISQ	= dabs(1d0-ETASQ)
	    COEF	= QOMS24*TSI**4d0
	    COEF1	= COEF/PSISQ**3.5d0
	    C2=COEF1*XNODP*(AODP*(1d0+1.5d0*ETASQ+EETA*(4d0+ETASQ))+0.75d0*CK2*TSI/PSISQ*X3THM1*(8d0+3d0*ETASQ*(8d0+ETASQ)))
	    C1		= BSTAR*C2
	    SINIO	= dsin(XINCL)
	    A3OVK2	= -XJ3/CK2*AE**3d0
	    C3		= COEF*TSI*A3OVK2*XNODP*AE*SINIO/EO
	    X1MTH2	= 1d0-THETA2
	    C4=2d0*XNODP*COEF1*AODP*BETAO2*(ETA*(2d0+0.5d0*ETASQ)+EO*(0.5d0+2d0*ETASQ)-2d0*CK2*TSI/
     &		(AODP*PSISQ)*(-3d0*X3THM1*(1d0-2d0*EETA+ETASQ*(1.5d0-0.5d0*EETA))+0.75d0*X1MTH2*(2d0*ETASQ-EETA*
     &		(1d0+ETASQ))*dcos(2d0*OMEGAO)))
	    C5=2d0*COEF1*AODP*BETAO2*(1d0+2.75d0*(ETASQ+EETA)+EETA*ETASQ)
	    THETA4=THETA2*THETA2
	    TEMP1=3d0*CK2*PINVSQ*XNODP
	    TEMP2=TEMP1*CK2*PINVSQ
	    TEMP3=1.25d0*CK4*PINVSQ*PINVSQ*XNODP
	    XMDOT=XNODP+0.5d0*TEMP1*BETAO*X3THM1+0.0625d0*TEMP2*BETAO*(13d0-78d0*THETA2+137d0*THETA4)
	    X1M5TH=1d0-5d0*THETA2
	    OMGDOT=-0.5d0*TEMP1*X1M5TH+0.0625d0*TEMP2*(7d0-114d0*THETA2+395d0*THETA4)+TEMP3*(3d0-36d0*THETA2+49d0*THETA4)
	    XHDOT1=-TEMP1*COSIO
	    XNODOT=XHDOT1+(0.5d0*TEMP2*(4d0-19d0*THETA2)+2d0*TEMP3*(3d0-7d0*THETA2))*COSIO
	    OMGCOF=BSTAR*C3*dcos(OMEGAO)
	    XMCOF=-TOTHRD*COEF*BSTAR*AE/EETA
	    XNODCF=3.5d0*BETAO2*XHDOT1*C1
	    T2COF=1.5d0*C1
	    XLCOF=0.125d0*A3OVK2*SINIO*(3d0+5d0*COSIO)/(1d0+COSIO)
	    AYCOF=0.25d0*A3OVK2*SINIO
	    DELMO=(1d0+ETA*dcos(XMO))**3d0
	    SINMO=dsin(XMO)
	    X7THM1=7d0*THETA2-1d0
	    IF(ISIMP .NE. 1) THEN
		C1SQ=C1*C1
		D2=4d0*AODP*TSI*C1SQ
		TEMP=D2*TSI*C1/3d0
		D3=(17d0*AODP+S4)*TEMP
		D4=0.5d0*TEMP*AODP*TSI*(221d0*AODP+31d0*S4)*C1
		T3COF=D2+2d0*C1SQ
		T4COF=0.25d0*(3d0*D3+C1*(12d0*D2+10d0*C1SQ))
		T5COF=0.2d0*(3d0*D4+12d0*C1*D3+6d0*D2*D2+15d0*C1SQ*(2d0*D2+C1SQ))
	    END IF

	    IFLAG=0

	END IF

	! Update for secular gravity and atmospheric drag

 	XMDF=XMO+XMDOT*TSINCE
	OMGADF=OMEGAO+OMGDOT*TSINCE
	XNODDF=XNODEO+XNODOT*TSINCE
	OMEGA=OMGADF
	XMP=XMDF
	TSQ=TSINCE*TSINCE
	XNODE=XNODDF+XNODCF*TSQ
	TEMPA=1d0-C1*TSINCE
	TEMPE=BSTAR*C4*TSINCE
	TEMPL=T2COF*TSQ
	IF(ISIMP .NE. 1) THEN
	    DELOMG=OMGCOF*TSINCE
	    DELM=XMCOF*((1d0+ETA*dcos(XMDF))**3d0-DELMO)
	    TEMP=DELOMG+DELM
	    XMP=XMDF+TEMP
	    OMEGA=OMGADF-TEMP
	    TCUBE=TSQ*TSINCE
	    TFOUR=TSINCE*TCUBE
	    TEMPA=TEMPA-D2*TSQ-D3*TCUBE-D4*TFOUR
	    TEMPE=TEMPE+BSTAR*C5*(dsin(XMP)-SINMO)
	    TEMPL=TEMPL+T3COF*TCUBE+TFOUR*(T4COF+TSINCE*T5COF)
	END IF
	A=AODP*TEMPA*TEMPA
	E=EO-TEMPE
	XL=XMP+OMEGA+XNODE+XNODP*TEMPL
	BETA=dsqrt(1d0-E*E)
	XN=XKE/A**1.5d0

	! Long period periodics

	AXN=E*dcos(OMEGA)
	TEMP=1d0/(A*BETA*BETA)
	XLL=TEMP*XLCOF*AXN
	AYNL=TEMP*AYCOF
	XLT=XL+XLL
	AYN=E*dsin(OMEGA)+AYNL

	! Solve keplers equation

	CAPU=dmod(dmod(XLT-XNODE,TWOPI)+TWOPI,TWOPI)
	TEMP2=CAPU
	DO I=1,10
	    SINEPW=dsin(TEMP2)
	    COSEPW=dcos(TEMP2)
	    TEMP3=AXN*SINEPW
	    TEMP4=AYN*COSEPW
	    TEMP5=AXN*COSEPW
	    TEMP6=AYN*SINEPW
	    EPW=(CAPU-TEMP4+TEMP3-TEMP2)/(1d0-TEMP5-TEMP6)+TEMP2
	    IF (dabs(EPW-TEMP2) .LE. E6A) GO TO 140
	    TEMP2=EPW
	END DO

	! Short period preliminary quantities

 140	ECOSE=TEMP5+TEMP6
	ESINE=TEMP3-TEMP4
	ELSQ=AXN*AXN+AYN*AYN
	TEMP=1d0-ELSQ
	PL=A*TEMP
	R=A*(1d0-ECOSE)
	TEMP1=1d0/R
	RDOT=XKE*dsqrt(A)*ESINE*TEMP1
	RFDOT=XKE*dsqrt(PL)*TEMP1
	TEMP2=A*TEMP1
	BETAL=dsqrt(TEMP)
	TEMP3=1d0/(1d0+BETAL)
	COSU=TEMP2*(COSEPW-AXN+AYN*ESINE*TEMP3)
	SINU=TEMP2*(SINEPW-AYN-AXN*ESINE*TEMP3)
	U=dmod(datan2(SINU,COSU)+TWOPI,TWOPI)
	SIN2U=2d0*SINU*COSU
	COS2U=2d0*COSU*COSU-1d0
	TEMP=1d0/PL
	TEMP1=CK2*TEMP
	TEMP2=TEMP1*TEMP

	! Update for short periodics

	RK=R*(1d0-1.5d0*TEMP2*BETAL*X3THM1)+0.5d0*TEMP1*X1MTH2*COS2U
	UK=U-0.25d0*TEMP2*X7THM1*SIN2U
	XNODEK=XNODE+1.5d0*TEMP2*COSIO*SIN2U
	XINCK=XINCL+1.5d0*TEMP2*COSIO*SINIO*COS2U
	RDOTK=RDOT-XN*TEMP1*X1MTH2*SIN2U
	RFDOTK=RFDOT+XN*TEMP1*(X1MTH2*COS2U+1.5d0*X3THM1)

	! Orientation vectors

	SINUK=dsin(UK)
	COSUK=dcos(UK)
	SINIK=dsin(XINCK)
	COSIK=dcos(XINCK)
	SINNOK=dsin(XNODEK)
	COSNOK=dcos(XNODEK)
	XMX=-SINNOK*COSIK
	XMY=COSNOK*COSIK
	UX=XMX*SINUK+COSNOK*COSUK
	UY=XMY*SINUK+SINNOK*COSUK
	UZ=SINIK*SINUK
	VX=XMX*COSUK-COSNOK*SINUK
	VY=XMY*COSUK-SINNOK*SINUK
	VZ=SINIK*COSUK

	! Position and velocity

	X=RK*UX
	Y=RK*UY
	Z=RK*UZ
	XDOT=RDOTK*UX+RFDOTK*VX
	YDOT=RDOTK*UY+RFDOTK*VY
	ZDOT=RDOTK*UZ+RFDOTK*VZ

	RETURN
	END
