c*******************************************************************************
c$ Component_name:
c	ORBIT
c$ Abstract:
c	Sample program to tabulate satellite positions as a function of time.
c$ Keywords:
c	KEPLER, ORBITAL_MOTION
c	FORTRAN, SAMPLE_PROGRAM
c$ Declarations:
c	none
c$ Inputs:
c	none
c$ Outputs:
c	none
c$ Returns:
c	none
c$ Detailed_description:
c	This simple program uses the KEPLER library to write out the coordinates
c	of a Saturnian satellite at uniformly-spaced time intervals.  Output may
c	be directed either to the terminal or a to file.
c
c	First, the user is asked to enter the orbital elements of a satellite.
c	Then the user is prompted repeatedly for a start time, stop time and
c	time step, all assumed to be in units of the orbital period.  In each
c	case, the user may send the tabulated coordinates to either the terminal
c	or a file.  Enter <EOF> at the terminal when done.
c$ External_references:
c	FKEP_SETPLANET(), FKEP_SETORBIT(), FKEP_MEANANOM(), FKEP_OMEGA(), 
c	FKEP_LOCATE()
c$ Examples:
c	none
c$ Error_handling:
c	none
c$ Limitations:
c	none
c$ Author_and_institution:
c	Mark R. Showalter
c	NASA/Ames Research Center
c$ Version_and_date:
c	1991 June 18
c$ Change_history:
c	none
c*******************************************************************************

	program ORBIT

	real*8		a, e, i, peri, node, trueLon, meanLon, period
	real*8		t, t1, t2, dt, xyz(3), vel(3), PI
	real*8		FKep_Omega, FKep_MeanAnom
	data		PI/ 3.141592653589793 /
	logical		outfile
	character*80	filename

c************************************************************************
c GM is from J. Anderson 1983, private communication to Phil Nicholson.
c
c J2, J4... are from P.D. Nicholson and C. C. Porco 1988.  A New 
c	Constraint on Saturn's Zonal Gravity Harmonics from Voyager
c	Observations of an Eccentric Ringlet.  J. Geophys. Res. 93, 1988.
c************************************************************************

	integer		NJs
	parameter	(NJs = 6)

	real*8		Rs, GM, Js(NJs)
	data		Rs/ 60330.d0 /
	data		GM/ 3.7931200d7 /
	data		Js/ 16297.d-6, -910.d-6, 107.d-6,
     &			      -10.d-6,    2.d-6, -0.5d-6 /

c***********************************************************************
c Specify planetary gravity field
c***********************************************************************

	call FKep_SetPlanet( Rs, GM, Js, NJs )

c***********************************************************************
c Specify satellite orbit
c***********************************************************************

	write(*,10) 'Semimajor axis a (km) = '
	read(*,*) a
	write(*,10) 'Eccentricity e = '
	read(*,*) e
	write(*,10) 'Inclination i (deg) = '
	read(*,*) i
	write(*,10) 'Longitude of pericenter (deg) = '
	read(*,*) peri
	write(*,10) 'Ascending node (deg) = '
	read(*,*) node
	write(*,10) 'True longitude (deg) = '
	read(*,*) trueLon
10	format( 1x, a, $ )

c Convert angles to radians
	i    = i    * PI/180.d0
	peri = peri * PI/180.d0
	node = node * PI/180.d0
	trueLon = trueLon * PI/180.d0
	
c Calculate mean longitude and orbital period
	meanLon = FKep_MeanAnom( e, trueLon-peri ) + peri 
	period = 2.*PI / FKep_Omega(a)

	call FKep_SetOrbit( a, e, i, peri, node, meanLon, 0.d0 )

c***********************************************************************
c MAIN LOOP
c***********************************************************************

100	continue

c Get time limits for listing
	write(*,10) ' '
	write(*,10) 'start, stop, delta time (periods) = '
	read(*,*,end=999) t1, t2, dt
	t2 = t2 + dt/100.

c Open output file
	write(*,10) 'Output filename (or <cr> for terminal): '
	read(*,11,end=999) filename
11	format(a)
	outfile = (filename .ne. ' ')
	if (outfile)
     &		open( 1, file=filename, status='new' )

c Write satellite positions
	do t = t1, t2, dt
		call FKep_Locate( t*period, xyz, vel )
		if (outfile) then
			write(1,20) t, xyz(1), xyz(2), xyz(3)
		else
			write(*,20) t, xyz(1), xyz(2), xyz(3)
		end if
20		format( 1x, f12.4, 3f16.2 )
	end do

	if (outfile) close(1)
	goto 100

c***********************************************************************
c END OF MAIN LOOP
c***********************************************************************

999	continue
	stop
	end

c***********************************************************************