* SGP 31 OCT 80

	SUBROUTINE SGPD(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 AO, QO, XLO, D1O, D2O, D3O, D4O, OMGDT, XNODOT, C5, C6

	IF (IFLAG .NE. 0) THEN

	    ! Initialization

	    C1= CK2*1.5d0
	    C2= CK2/4d0
	    C3= CK2/2d0
	    C4= XJ3*AE**3d0/(4d0*CK2)
	    COSIO=dcos(XINCL)
	    SINIO=dsin(XINCL)
	    A1=(XKE/XNO)**TOTHRD
	    D1= C1/A1/A1*(3d0*COSIO*COSIO-1d0)/(1d0-EO*EO)**1.5d0
	    AO=A1*(1d0-1d0/3d0*D1-D1*D1-134d0/81d0*D1*D1*D1)
	    PO=AO*(1d0-EO*EO)
	    QO=AO*(1d0-EO)
	    XLO=XMO+OMEGAO+XNODEO
	    D1O= C3 *SINIO*SINIO
	    D2O= C2 *(7d0*COSIO*COSIO-1d0)
	    D3O=C1*COSIO
	    D4O=D3O*SINIO
	    PO2NO=XNO/(PO*PO)
	    OMGDT=C1*PO2NO*(5d0*COSIO*COSIO-1d0)
	    XNODOT=-2d0*D3O*PO2NO
	    C5=0.5d0*C4*SINIO*(3d0+5d0*COSIO)/(1d0+COSIO)
	    C6=C4*SINIO
	    IFLAG=0

	END IF

	! Update for secular gravity and atmospheric drag
 
	A=XNO+(2d0*XNDT2O+3d0*XNDD6O*TSINCE)*TSINCE
  	A=AO*(XNO/A)**TOTHRD
	E=E6A
	IF(A.GT.QO) E=1d0-QO/A
	P=A*(1d0-E*E)
	XNODES= XNODEO+XNODOT*TSINCE
	OMGAS= OMEGAO+OMGDT*TSINCE
	XLS=dmod(dmod(XLO+(XNO+OMGDT+XNODOT+(XNDT2O+XNDD6O*TSINCE)*TSINCE)*TSINCE,TWOPI)+TWOPI,TWOPI)

	! Long period periodics

	AXNSL=E*dcos(OMGAS)
 	AYNSL=E*dsin(OMGAS)-C6/P
	XL=dmod(dmod(XLS-C5/P*AXNSL,TWOPI)+TWOPI,TWOPI)

	! Solve keplers equation

	U=dmod(dmod(XL-XNODES,TWOPI)+TWOPI,TWOPI)
	ITEM3=0
	EO1=U
	TEM5=1d0
 20	SINEO1=dsin(EO1)
 	COSEO1=dcos(EO1)
	IF (dabs(TEM5) .LT. E6A .OR. ITEM3 .GE. 10) GO TO 30
	ITEM3=ITEM3+1
	TEM5=1d0-COSEO1*AXNSL-SINEO1*AYNSL
	TEM5=(U-AYNSL*COSEO1+AXNSL*SINEO1-EO1)/TEM5
	TEM2=dabs(TEM5)
	IF(TEM2.GT.1d0) TEM5=TEM2/TEM5
	EO1=EO1+TEM5
	GO TO 20

	! Short period preliminary quantities

 30	ECOSE=AXNSL*COSEO1+AYNSL*SINEO1
 	ESINE=AXNSL*SINEO1-AYNSL*COSEO1
	EL2=AXNSL*AXNSL+AYNSL*AYNSL
	PL=A*(1d0-EL2)
	PL2=PL*PL
	R=A*(1d0-ECOSE)
	RDOT=XKE*dsqrt(A)/R*ESINE
	RVDOT=XKE*dsqrt(PL)/R
	TEMP=ESINE/(1d0+dsqrt(1d0-EL2))
	SINU=A/R*(SINEO1-AYNSL-AXNSL*TEMP)
	COSU=A/R*(COSEO1-AXNSL+AYNSL*TEMP)
	SU=dmod(datan2(SINU,COSU)+TWOPI,TWOPI)

	! Update for short periodics

	SIN2U=(COSU+COSU)*SINU
	COS2U=1d0-2d0*SINU*SINU
	RK=R+D1O/PL*COS2U
	UK=SU-D2O/PL2*SIN2U
	XNODEK=XNODES+D3O*SIN2U/PL2
	XINCK =XINCL+D4O/PL2*COS2U

	! Orientation vectors

	SINUK=dsin(UK)
	COSUK=dcos(UK)
	SINNOK=dsin(XNODEK)
	COSNOK=dcos(XNODEK)
	SINIK=dsin(XINCK)
	COSIK=dcos(XINCK)
	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=RDOT*UX
	YDOT=RDOT*UY
	ZDOT=RDOT*UZ
	XDOT=RVDOT*VX+XDOT
	YDOT=RVDOT*VY+YDOT
	ZDOT=RVDOT*VZ+ZDOT

	RETURN
	END