* SGP8 14 NOV 80

	SUBROUTINE SGP8D(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

	DATA RHO/0.15696615d0/

	SAVE COSI, THETA2, TTHMUN, XNODP, ISIMP, SINI, SINIO2, COSIO2 
	SAVE UNM5TH, UNMTH2, A3COF, XMDT1, XGDT1, XHDT1, XLLDOT, OMGDT
	SAVE XNODOT, XNDT, EDOT, PP, GAMMA, XND, QQ, ED, OVGPP

	IF (IFLAG .NE. 0) THEN

	    ! Recover original mean motion (xnodp) and semimajor axis (aodp)
	    ! from input elements --------- calculate ballistic coefficient
	    ! (B term) from input B* drag term

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

	    ! Initialization

	    ISIMP=0
	    PO=AODP*BETAO2
	    POM2=1d0/(PO*PO)
	    SINI=dsin(XINCL)
	    SING=dsin(OMEGAO)
	    COSG=dcos(OMEGAO)
	    TEMP=0.5d0*XINCL
	    SINIO2=dsin(TEMP)
	    COSIO2=dcos(TEMP)
	    THETA4=THETA2*THETA2
	    UNM5TH=1d0-5d0*THETA2
	    UNMTH2=1d0-THETA2
	    A3COF=-XJ3/CK2*AE*AE*AE
	    PARDT1=3d0*CK2*POM2*XNODP
	    PARDT2=PARDT1*CK2*POM2
	    PARDT4=1.25d0*CK4*POM2*POM2*XNODP
	    XMDT1=0.5d0*PARDT1*BETAO*TTHMUN
	    XGDT1=-0.5d0*PARDT1*UNM5TH
	    XHDT1=-PARDT1*COSI
	    XLLDOT=XNODP+XMDT1+0.0625d0*PARDT2*BETAO*(13d0-78d0*THETA2+137d0*THETA4)
	    OMGDT=XGDT1+0.0625d0*PARDT2*(7d0-114d0*THETA2+395d0*THETA4)+PARDT4*(3d0-36d0*THETA2+49d0*THETA4)
	    XNODOT=XHDT1+(0.5d0*PARDT2*(4d0-19d0*THETA2)+2d0*PARDT4*(3d0-7d0*THETA2))*COSI
	    TSI=1d0/(PO-S)
	    ETA=EO*S*TSI
	    ETA2=ETA*ETA
	    PSIM2=dabs(1d0/(1d0-ETA2))
	    ALPHA2=1d0+EOSQ
	    EETA=EO*ETA
	    COS2G=2d0*COSG*COSG-1d0
	    D5=TSI*PSIM2
	    D1=D5/PO
	    D2=12d0+ETA2*(36d0+4.5d0*ETA2)
	    D3=ETA2*(15d0+2.5d0*ETA2)
	    D4=ETA*(5d0+3.75d0*ETA2)
	    B1=CK2*TTHMUN
	    B2=-CK2*UNMTH2
	    B3=A3COF*SINI
	    C0=0.5d0*B*RHO*QOMS2T*XNODP*AODP*TSI**4d0*PSIM2**3.5d0/dsqrt(ALPHA2)
	    C1=1.5d0*XNODP*ALPHA2*ALPHA2*C0
	    C4=D1*D3*B2
	    C5=D5*D4*B3
	    XNDT=C1*((2d0+ETA2*(3d0+34d0*EOSQ)+5d0*EETA*(4d0+ETA2)+8.5d0*EOSQ)+D1*D2*B1+C4*COS2G+C5*SING)
	    XNDTN=XNDT/XNODP

	    ! If drag is very small, the isimp flag is set and the
	    ! equations are truncated to linear variation in mean
	    ! motion and quadratic variation in mean anomaly

	    IF(dabs(XNDTN*XMNPDA) .GE. 2.16d-3) THEN

		D6=ETA*(30d0+22.5d0*ETA2)
		D7=ETA*(5d0+12.5d0*ETA2)
		D8=1d0+ETA2*(6.75d0+ETA2)
		C8=D1*D7*B2
		C9=D5*D8*B3
		EDOT=-C0*(ETA*(4d0+ETA2+EOSQ*(15.5d0+7d0*ETA2))+EO*(5d0+15d0*ETA2)+D1*D6*B1 +C8*COS2G+C9*SING)
		D20=0.5d0*TOTHRD*XNDTN
		ALDTAL=EO*EDOT/ALPHA2
		TSDTTS=2d0*AODP*TSI*(D20*BETAO2+EO*EDOT)
		ETDT=(EDOT+EO*TSDTTS)*TSI*S
		PSDTPS=-ETA*ETDT*PSIM2
		SIN2G=2d0*SING*COSG
		C0DTC0=D20+4d0*TSDTTS-ALDTAL-7d0*PSDTPS
		C1DTC1=XNDTN+4d0*ALDTAL+C0DTC0
		D9=ETA*(6d0+68d0*EOSQ)+EO*(20d0+15d0*ETA2)
		D10=5d0*ETA*(4d0+ETA2)+EO*(17d0+68d0*ETA2)
		D11=ETA*(72d0+18d0*ETA2)
		D12=ETA*(30d0+10d0*ETA2)
		D13=5d0+11.25d0*ETA2
		D14=TSDTTS-2d0*PSDTPS
		D15=2d0*(D20+EO*EDOT/BETAO2)
		D1DT=D1*(D14+D15)
		D2DT=ETDT*D11
		D3DT=ETDT*D12
		D4DT=ETDT*D13
		D5DT=D5*D14
		C4DT=B2*(D1DT*D3+D1*D3DT)
		C5DT=B3*(D5DT*D4+D5*D4DT)
		D16= D9*ETDT+D10*EDOT +B1*(D1DT*D2+D1*D2DT) +C4DT*COS2G+C5DT*SING+XGDT1*(C5*COSG-2d0*C4*SIN2G)
		XNDDT=C1DTC1*XNDT+C1*D16
		EDDOT=C0DTC0*EDOT-C0*(
     &		    (4d0+3d0*ETA2+30d0*EETA+EOSQ*(15.5d0+21d0*ETA2))*ETDT+(5d0+15d0*ETA2+
     &		    EETA*(31d0+14d0*ETA2))*EDOT + B1*(D1DT*D6+D1*ETDT*(30d0+67.5d0*ETA2)) +
     &		    B2*(D1DT*D7+D1*ETDT*(5d0+37.5d0*ETA2))*COS2G+
     &		    B3*(D5DT*D8+D5*ETDT*ETA*(13.5d0+4d0*ETA2))*SING+XGDT1*(C9*COSG-2d0*C8*SIN2G))
		D25=EDOT*EDOT
		D17=XNDDT/XNODP-XNDTN*XNDTN
		TSDDTS=2d0*TSDTTS*(TSDTTS-D20)+AODP*TSI*(TOTHRD*BETAO2*D17-4d0*D20*EO*EDOT+2d0*(D25+EO*EDDOT))
		ETDDT =(EDDOT+2d0*EDOT*TSDTTS)*TSI*S+TSDDTS*ETA
		D18=TSDDTS-TSDTTS*TSDTTS
		D19=-PSDTPS*PSDTPS/ETA2-ETA*ETDDT*PSIM2-PSDTPS*PSDTPS
		D23=ETDT*ETDT
		D1DDT=D1DT*(D14+D15)+D1*(D18-2d0*D19+TOTHRD*D17+2d0*(ALPHA2*D25/BETAO2+EO*EDDOT)/BETAO2)
		XNTRDT=XNDT*(2d0*TOTHRD*D17+3d0*(D25+EO*EDDOT)/ALPHA2-6d0*ALDTAL*ALDTAL +4d0*D18-7d0*D19 ) +
     &		    C1DTC1*XNDDT+C1*(C1DTC1*D16+D9*ETDDT+D10*EDDOT+D23*(6d0+30d0*EETA+68d0*EOSQ)+
     &		    ETDT*EDOT*(40d0+30d0*ETA2+272d0*EETA)+D25*(17d0+68d0*ETA2) +
     &		    B1*(D1DDT*D2+2d0*D1DT*D2DT+D1*(ETDDT*D11+D23*(72d0+54d0*ETA2))) +
     &		    B2*(D1DDT*D3+2d0*D1DT*D3DT+D1*(ETDDT*D12+D23*(30d0+30d0*ETA2))) * COS2G+
     &		    B3*((D5DT*D14+D5*(D18-2d0*D19)) *D4+2d0*D4DT*D5DT+D5*(ETDDT*D13+22.5d0*ETA*D23)) *SING+XGDT1*
     &		    ((7d0*D20+4d0*EO*EDOT/BETAO2)*(C5*COSG-2d0*C4*SIN2G)+
     &		    ((2d0*C5DT*COSG-4d0*C4DT*SIN2G)-XGDT1*(C5*SING+4d0*C4*COS2G))))
		TMNDDT=XNDDT*1d9
		TEMP=TMNDDT*TMNDDT-XNDT*1d18*XNTRDT
		PP=(TEMP+TMNDDT*TMNDDT)/TEMP
		GAMMA=-XNTRDT/(XNDDT*(PP-2d0))
		XND=XNDT/(PP*GAMMA)
		QQ=1d0-EDDOT/(EDOT*GAMMA)
		ED=EDOT/(QQ*GAMMA)
		OVGPP=1d0/(GAMMA*(PP+1d0))

	    ELSE

		ISIMP=1
		EDOT=-TOTHRD*XNDTN*(1d0-EO)

	    END IF

	    IFLAG=0

	END IF

	! Update for secular gravity and atmospheric drag

	XMAM=dmod(dmod(XMO+XLLDOT*TSINCE,TWOPI)+TWOPI,TWOPI)
 	OMGASM=OMEGAO+OMGDT*TSINCE
	XNODES=XNODEO+XNODOT*TSINCE
	IF (ISIMP .NE. 1) THEN
	    TEMP=1d0-GAMMA*TSINCE
	    TEMP1=TEMP**PP
	    XN=XNODP+XND*(1d0-TEMP1)
	    EM=EO+ED*(1d0-TEMP**QQ)
	    Z1=XND*(TSINCE+OVGPP*(TEMP*TEMP1-1d0))
	ELSE
 	    XN=XNODP+XNDT*TSINCE
 	    EM=EO+EDOT*TSINCE
	    Z1=0.5d0*XNDT*TSINCE*TSINCE
	END IF
 	Z7=3.5d0*TOTHRD*Z1/XNODP
 	XMAM=dmod(dmod(XMAM+Z1+Z7*XMDT1,TWOPI)+TWOPI,TWOPI)
	OMGASM=OMGASM+Z7*XGDT1
	XNODES=XNODES+Z7*XHDT1

	! Solve keplers equation

	ZC2=XMAM+EM*dsin(XMAM)*(1d0+EM*dcos(XMAM))
	DO I=1,10
	    SINE=dsin(ZC2)
	    COSE=dcos(ZC2)
	    ZC5=1d0/(1d0-EM*COSE)
	    CAPE=(XMAM+EM*SINE-ZC2)*ZC5+ZC2
	    IF(dabs(CAPE-ZC2) .LE. E6A) GO TO 140
	    ZC2=CAPE
 	END DO

	! Short period preliminary quantities

 140	AM=(XKE/XN)**TOTHRD
 	BETA2M=1d0-EM*EM
	SINOS=dsin(OMGASM)
	COSOS=dcos(OMGASM)
	AXNM=EM*COSOS
	AYNM=EM*SINOS
	PM=AM*BETA2M
	G1=1d0/PM
	G2=0.5d0*CK2*G1
	G3=G2*G1
	BETA=dsqrt(BETA2M)
	G4=0.25d0*A3COF*SINI
	G5=0.25d0*A3COF*G1
	SNF=BETA*SINE*ZC5
	CSF=(COSE-EM)*ZC5
	FM=dmod(datan2(SNF,CSF)+TWOPI,TWOPI)
	SNFG=SNF*COSOS+CSF*SINOS
	CSFG=CSF*COSOS-SNF*SINOS
	SN2F2G=2d0*SNFG*CSFG
	CS2F2G=2d0*CSFG*CSFG-1d0
	ECOSF=EM*CSF
	G10=FM-XMAM+EM*SNF
	RM=PM/(1d0+ECOSF)
	AOVR=AM/RM
	G13=XN*AOVR
	G14=-G13*AOVR
	DR=G2*(UNMTH2*CS2F2G-3d0*TTHMUN)-G4*SNFG
	DIWC=3d0*G3*SINI*CS2F2G-G5*AYNM
	DI=DIWC*COSI

	! Update for short period periodics

	SNI2DU=SINIO2*(G3*(0.5d0*(1d0-7d0*THETA2)*SN2F2G-3d0*UNM5TH*G10)-G5*SINI*CSFG*(2d0+
     &	    ECOSF))-0.5d0*G5*THETA2*AXNM/COSIO2
	XLAMB=FM+OMGASM+XNODES+G3*(0.5d0*(1d0+6d0*COSI-7d0*THETA2)*SN2F2G-3d0*
     &	    (UNM5TH+2d0*COSI)*G10)+G5*SINI*(COSI*AXNM/(1d0+COSI)-(2d0+ECOSF)*CSFG)
	Y4=SINIO2*SNFG+CSFG*SNI2DU+0.5d0*SNFG*COSIO2*DI
	Y5=SINIO2*CSFG-SNFG*SNI2DU+0.5d0*CSFG*COSIO2*DI
	R=RM+DR
	RDOT=XN*AM*EM*SNF/BETA+G14*(2d0*G2*UNMTH2*SN2F2G+G4*CSFG)
	RVDOT=XN*AM*AM*BETA/RM+G14*DR+AM*G13*SINI*DIWC

	! Orientation vectors

	SNLAMB=dsin(XLAMB)
	CSLAMB=dcos(XLAMB)
	TEMP=2d0*(Y5*SNLAMB-Y4*CSLAMB)
	UX=Y4*TEMP+CSLAMB
	VX=Y5*TEMP-SNLAMB
	TEMP=2d0*(Y5*CSLAMB+Y4*SNLAMB)
	UY=-Y4*TEMP+SNLAMB
	VY=-Y5*TEMP+CSLAMB
	TEMP=2d0*dsqrt(1d0-Y4*Y4-Y5*Y5)
	UZ=Y4*TEMP
	VZ=Y5*TEMP

	! Position and velocity

	X=R*UX
	Y=R*UY
	Z=R*UZ
	XDOT=RDOT*UX+RVDOT*VX
	YDOT=RDOT*UY+RVDOT*VY
	ZDOT=RDOT*UZ+RVDOT*VZ

	RETURN
	END