* SDP4 3 NOV 80

	SUBROUTINE SDP4D(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, C1, SINIO, C4, XMDOT, OMGDOT, XNODOT, XNODCF
	SAVE T2COF, XLCOF, AYCOF, X7THM1

	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 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)
	    SING=dsin(OMEGAO)
	    COSG=dcos(OMEGAO)
	    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*AE*AE
	    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)))
	    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
	    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
	    X7THM1=7d0*THETA2-1d0
	    IFLAG=0
	    CALL DPINITD(EOSQ,SINIO,COSIO,BETAO,AODP,THETA2,SING,COSG,BETAO2,XMDOT,OMGDOT,XNODOT,XNODP)

	END IF

	! Update for secular gravity and atmospheric drag

	XMDF=XMO+XMDOT*TSINCE
 	OMGADF=OMEGAO+OMGDOT*TSINCE
	XNODDF=XNODEO+XNODOT*TSINCE
	TSQ=TSINCE*TSINCE
	XNODE=XNODDF+XNODCF*TSQ
	TEMPA=1d0-C1*TSINCE
	TEMPE=BSTAR*C4*TSINCE
	TEMPL=T2COF*TSQ
	XN=XNODP
	CALL DPSECD(XMDF,OMGADF,XNODE,EM,XINC,XN,OMGDOT,TSINCE)
	A=(XKE/XN)**TOTHRD*TEMPA*TEMPA
	E=EM-TEMPE
	XMAM=XMDF+XNODP*TEMPL
	CALL DPPERD(SINIO,COSIO,E,XINC,OMGADF,XNODE,XMAM,TSINCE)
	XL=XMAM+OMGADF+XNODE
	BETA=dsqrt(1d0-E*E)
	XN=XKE/A**1.5d0

	! Long period periodics

	AXN=E*dcos(OMGADF)
	TEMP=1d0/(A*BETA*BETA)
	XLL=TEMP*XLCOF*AXN
	AYNL=TEMP*AYCOF
	XLT=XL+XLL
	AYN=E*dsin(OMGADF)+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=XINC+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