* DRIVER 3 NOV 80
* WGS-72 PHYSICAL AND GEOPOTENTIAL CONSTANTS
* CK2= .5*J2*AE**2 CK4=-.375*J4*AE**4

	program tle_test

	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 DE2RA,E6A,PI,PIO2,QO,SO,TOTHRD,TWOPI,X3PIO2,XJ2,XJ3,XJ4,XKE,XKMPER,XMNPDA,AE
     &	    /.174532925d-1,1d-6,3.14159265d0,1.57079633d0,120d0,78d0,0.66666667d0,6.2831853d0,
     &	     4.71238898d0,1.082616d-3,-0.253881d-5,-1.65597d-6,0.743669161d-1,6378.135d0,1440d0,1d0/

	CHARACTER ABUF(2)*80
	CHARACTER CSET(5)*4	/'SGP ','SGP4','SDP4','SGP8','SDP8'/

	include		'openfile.h'
	include		'planet.h'
	include		'math.h'

	logical		bOpenFile
	character	cStr*80
	character	cFile*80
	parameter	(nVar=1)
	character	cVar(nVar)*20
	character	cArg*20
	integer		time(2)
	logical		bReset
	logical		ForeignArgSet

	double precision JD

	XKMPER = PLANET__EARTH_RADIUS
	PI   = MATH__PI
	PIO2 = PI/2d0
	TWOPI = 2d0*PI

	! SELECT EPHEMERIS TYPE AND OUTPUT TIMES

	CK2=0.5d0*XJ2*AE**2d0
	CK4=-0.375d0*XJ4*AE**4d0
	QOMS2T=((QO-SO)*AE/XKMPER)**4d0
	S=AE*(1d0+SO/XKMPER)

	iVar = nVar
	call ForeignArg(' ',iVar,cVar,cArg)
	call ForeignI4Arg(cArg,'ephemeris',2,IEPT0)
	bReset = ForeignArgSet(cArg,'reset_flag')

	TBeg = 2453104.5d0 !2452541.5d0
	TEnd = TBeg!+3
	DelT = 1.0d0

	! READ IN MEAN ELEMENTS FROM 2 CARD T(TRANS) OR G(INTERN) FORMAT

	cFile = 'tle.txt'
	iRecl = 0
	if (.not. bOpenFile(OPN__TEXT+OPN__REOPEN,iU,cFile,iRecl)) continue

	WRITE (*,'(10X,5H TIME,12X,1HX,14X,1HY,14X,1HZ,9X,4HXDOT,6X,4HYDOT,6X,4HZDOT)')

	read (iU,'(A)',iostat=I) cStr

	do while (I .eq. 0)

	    READ (iU,'(A)') cStr
	    READ (cStr,'(18X,D14.8,1X,F10.8,2(1X,F6.5,I2))') EPOCH,XNDT2O,XNDD6O,IEXP,BSTAR,IBEXP
	    READ (iU,'(7X,2(1X,F8.4),1X,F7.7,2(1X,F8.4),1X,F11.8)') XINCL,XNODEO,EO,OMEGAO,XMO,XNO

	    READ (cStr(19:20),'(I2.2)' ) iYr
	    READ (cStr(21:32),'(F12.8)') JD 

	    iYr = 2000+iYr
	    call Time2Day8(1,time,JD)
	    call Time2YDoy(1,time,iYr,time)
	    call Time2JD(0,time,time)
	    call Time2Day8(0,time,JD)

	    XNDD6O=XNDD6O*(10d0**IEXP)
	    XNODEO=XNODEO*DE2RA
	    OMEGAO=OMEGAO*DE2RA
	    XMO=XMO*DE2RA
	    XINCL=XINCL*DE2RA
	    TEMP=TWOPI/XMNPDA/XMNPDA
	    XNO=XNO*TEMP*XMNPDA
	    XNDT2O=XNDT2O*TEMP
	    XNDD6O=XNDD6O*TEMP/XMNPDA

	    ! INPUT CHECK FOR PERIOD VS EPHEMERIS SELECTED
	    ! PERIOD GE 225 MINUTES IS DEEP SPACE

	    A1=(XKE/XNO)**TOTHRD
	    TEMP=1.5d0*CK2*(3d0*dcos(XINCL)**2d0-1d0)/(1d0-EO*EO)**1.5d0
	    DEL1=TEMP/(A1*A1)
	    AO=A1*(1d0-DEL1*(0.5d0*TOTHRD+DEL1*(1d0+134d0/81d0*DEL1)))
	    DELO=TEMP/(AO*AO)
	    XNODP=XNO/(1d0+DELO)
	    IDEEP = 0
	    IF((TWOPI/XNODP/XMNPDA) .GE. 0.15625d0) IDEEP=1
	    BSTAR = BSTAR*(10d0**IBEXP)/AE

	    TT = TBeg

	    IEPT = IEPT0
	    IF (IEPT .NE. 0) THEN
		IF(IDEEP .EQ. 1 .AND. (IEPT .EQ. 1 .OR. IEPT .EQ. 2 .OR. IEPT .EQ. 4)) THEN
		    !IEPT = IEPT+1
		    WRITE(6,'(2A)') ' SWITCHING TO DEEP SPACE EPHEMERIS ',CSET(IEPT)
		ELSE IF(IDEEP .EQ. 0 .AND. (IEPT .EQ. 3 .OR. IEPT .EQ. 5)) THEN
		    !IEPT = IEPT-1
		    WRITE(6,'(2A)') ' SWITCHING TO NEAR EARTH EPHEMERIS ',CSET(IEPT)
 		END IF
	    END IF

	    IFLAG_SGP  = 1
	    IFLAG_SGP4 = 1
	    IFLAG_SDP4 = 1
	    IFLAG_SGP8 = 1
	    IFLAG_SDP8 = 1


	    DO WHILE (TT .LE. TEnd)

		TSINCE = (TT-JD)*1440d00

		IF (IEPT .EQ. 0) THEN

		    CALL SGPD(IFLAG_SGP ,TSINCE)
		    CALL TLE_SHOW(CSET(IEPT),JD,TT,X,Y,Z,XDOT,YDOT,ZDOT,XKMPER,AE,XMNPDA)
		    CALL SGP4D(IFLAG_SGP4,TSINCE)
		    CALL TLE_SHOW(CSET(IEPT),JD,TT,X,Y,Z,XDOT,YDOT,ZDOT,XKMPER,AE,XMNPDA)
		    !CALL SDP4D(IFLAG_SDP4,TSINCE)
		    !CALL TLE_SHOW(CSET(IEPT),JD,TT,X,Y,Z,XDOT,YDOT,ZDOT,XKMPER,AE,XMNPDA)
		    CALL SGP8D(IFLAG_SGP8,TSINCE)
		    CALL TLE_SHOW(CSET(IEPT),JD,TT,X,Y,Z,XDOT,YDOT,ZDOT,XKMPER,AE,XMNPDA)
		    !CALL SDP8D(IFLAG_SDP8,TSINCE)
		    !CALL TLE_SHOW(CSET(IEPT),JD,TT,X,Y,Z,XDOT,YDOT,ZDOT,XKMPER,AE,XMNPDA)

		ELSE IF (IEPT .EQ. 1) THEN
		    CALL SGPD(IFLAG_SGP,TSINCE)
		    CALL TLE_SHOW(CSET(IEPT),JD,TT,X,Y,Z,XDOT,YDOT,ZDOT,XKMPER,AE,XMNPDA)

		ELSE IF (IEPT .EQ. 2) THEN
		    CALL SGP4D(IFLAG_SGP4,TSINCE)
		    CALL TLE_SHOW(CSET(IEPT),JD,TT,X,Y,Z,XDOT,YDOT,ZDOT,XKMPER,AE,XMNPDA)

		ELSE IF (IEPT .EQ. 3) THEN
		    CALL SDP4D(IFLAG_SDP4,TSINCE)
		    CALL TLE_SHOW(CSET(IEPT),JD,TT,X,Y,Z,XDOT,YDOT,ZDOT,XKMPER,AE,XMNPDA)

		ELSE IF (IEPT .EQ. 4) THEN
		    CALL SGP8D(IFLAG_SGP8,TSINCE)
		    CALL TLE_SHOW(CSET(IEPT),JD,TT,X,Y,Z,XDOT,YDOT,ZDOT,XKMPER,AE,XMNPDA)

		ELSE IF (IEPT .EQ. 5) THEN
		    CALL SDP8D(IFLAG_SDP8,TSINCE)
		    CALL TLE_SHOW(CSET(IEPT),JD,TT,X,Y,Z,XDOT,YDOT,ZDOT,XKMPER,AE,XMNPDA)

		END IF

		IF (bReset) THEN
		    IFLAG_SGP  = 1
		    IFLAG_SGP4 = 1
		    IFLAG_SDP4 = 1
		    IFLAG_SGP8 = 1
		    IFLAG_SDP8 = 1
		END IF

		TT=TT+DelT

	    END DO

	    read (iU,'(A)',iostat=I) cStr

	end do

	END

	subroutine TLE_SHOW(CEPH,JD,TT,X,Y,Z,XDOT,YDOT,ZDOT,XKMPER,AE,XMNPDA)

	implicit double precision (A-H,O-Z)
	double precision JD

	X   =X   *XKMPER/AE
	Y   =Y   *XKMPER/AE
	Z   =Z   *XKMPER/AE
	XDOT=XDOT*XKMPER/AE*XMNPDA/86400.
	YDOT=YDOT*XKMPER/AE*XMNPDA/86400.
	ZDOT=ZDOT*XKMPER/AE*XMNPDA/86400.

	WRITE (*,'(A4,F11.5,F11.5,3F15.5,3F10.5)') CEPH,TT-245d4,TT-JD,X,Y,Z,XDOT,YDOT,ZDOT

	RETURN
	END