C+
C NAME:
C	Time2KeplerOrbit
C PURPOSE:
C	Calculate orbital positions and velocities for elliptic and
C	hyperbolic Kepler orbits.
C CATEGORY:
C	Celestial mechanics: two-body problem
C CALLING SEQUENCE:
	subroutine Time2KeplerOrbit(tt,m1i,m2i,OrbitElement,Position,Velocity)
C INPUTS:
C	tt(2)		integer		time difference with epoch t0
C	m1		double precision primary mass (in solar masses)
C					if m1 <=0 then m1 = 1 (solar mass) is used
C	m2		double precision secondary mass (in solar masses)
C					if m2 < 0 then m2 = 0 is used. 
C	OrbitElement(6)	double precision orbital elements at epoch t0:
C					1: semi-major axis (AU)
C					2: eccentricity
C					3: inclination of orbital plane to ecliptic (deg)
C					4: ecliptic longitude ascending node (deg)
C					5: longitude of perihelion (deg)
C					   = angle(longitude ascending node)+
C					     argument of perihelion
C					6: mean longitude at epoch t0 (deg)
C	!! Argument of perihelion = angle(ascending node-perihelion)
C	!! If the orbital inclination (element 3) is zero, then the longitude of the
C		node (element 4) is not used, and the longitude of the perihelion and
C		the mean longitude at t0 become ecliptic longitudes.
C OUTPUTS:
C	Position(3)	double precision position vector relative to m1:
C					1: ecliptic longitude (deg),
C					2: ecliptic latitude (deg),
C					3: radial distance (AU)
C	Velocity(5)	double precision velocity vector relative to m1:
C					1: ecliptic longitude (deg),
C					2: ecliptic latitude (deg),
C					3: velocity (AU/day),
C					4: tangential velocity (AU/day),
C					5: radial velocity (AU/day)
C EXTERNAL TYPE:
	double precision	EqKeplerd
C EXTERNAL:
	external		EqKeplerd
C INCLUDE:
	include		'math.h'
	include		'sun.h'
C COMMON TYPE:
	double precision	M
	double precision	ee
C COMMON BLOCKS:
	common	/KEPLERD/ 	M,ee
C CALLS:
C	Time2Day8, nrZBrentd, rotated
C PROCEDURE:
C >	Basic two-body problem (see Green, p. 163-164)
C >	The time of perihelion passage can be used as epoch t0.
C	In that case the mean longitude (at epoch t0) is equal to
C	the longitude of the perihelion; i.e.
C		OrbitElement(5) = OrbitElement(6)
C MODIFICATION HISTORY:
C	JUN-1997, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C	Adapted from IDL procedure KeplerOrbit
C-
	    integer		tt(2)
	    double precision	m1i
	    double precision	m2i
	    double precision	OrbitElement(6)
	    double precision	Position(3)
	    double precision	Velocity(5)

	double precision	jd
	double precision	m1
	double precision	m2
	double precision	mu
	double precision	nn
	double precision	aa
	double precision	lr
	double precision	lv
	double precision	ii
	double precision	ll
	double precision	OO
	double precision	wb
	double precision	ww
	double precision	AU
	double precision	SpD
	double precision	pm1
	double precision	E
	double precision	r
	double precision	dr
	double precision	vv
	double precision	v
	double precision	vt
	double precision	vr
	double precision	dv
	double precision	g

	double precision	nrZBrentd

	double precision	dsind
	double precision	dcosd
	double precision	datan2d

	call Time2Day8(0,tt,jd)

	aa = OrbitElement(1)			! Semi-major axis
	ee = OrbitElement(2)			! Eccentricity
	ii = OrbitElement(3)			! Orbit plane inclination
	OO = OrbitElement(4)			! Ecliptic longitude ascending node
	wb = OrbitElement(5)			! Longitude of perihelion (degrees)
	ll = OrbitElement(6)			! Mean longitude at epoch t0
						! Only the difference ll-wb is used
	ww = wb
	if (ii .ne. 0.0d0) ww = wb-OO		! Angle ascending node to perihelion
						! (in orbital plane)
	m1 = m1i
	if (m1 .le. 0.0d0) m1 = 1.0d0		! Default: solar mass

	m2  = max(m2i,0.0d0)			! Default: m2 = 0.

	mu  = 0.1d0*SUN__GM*(m1+m2)		! units 10^27 cm^3/s^2

	AU  = SUN__AU				! units 10^13 cm
	SpD = 0.086400d0			! 10^6 Seconds per day

	pm1 = dsign(1.0d0,1.0d0-ee)		! +1 for ellipse; -1 for hyperbole

	nn = AU*aa
	nn = dsqrt(mu/(nn*nn*nn))*SpD*MATH__DPR	! Mean (angular) motion (degrees/day)

	M = nn*jd+ll-wb				! Mean anomaly (degrees)
	M = dmod(dmod(M,360.0d0)+360.0d0,360.0d0)! In [0,360]
						! Eccentric anomaly (degrees)
	E = nrZBrentd(EqKeplerd,0.0d0,360.0d0,0.00001d0)! Solve Kepler's equation

	if (ee .lt. 1.0d0) then			! Elliptical orbit
	    r = 1.0d0-ee*dcosd(E)		! Distance to m1 in units of aa
	    vv = datan2d(dsqrt(1.0d0-ee*ee)*dsind(E),dcosd(E)-ee)! True anomaly (degrees)
	else					! Hyperbolic orbit
	    E = E*MATH__RPD
	    r = ee*dcosh(E)-1.0d0
	    vv = datan2d(dsqrt(ee*ee-1.0d0)*dsinh(E),ee-dcosh(E))! True anomaly (degrees)
	end if

	lr = ww+vv				! Longitude radius vector
	lr = dmod(dmod(lr,360.0d0)+360.0d0,360.0d0)

	dr = 0.0d0				! Latitude radius vector

	v  = 2.0d0/r-pm1			! Velocity^2
	vt = pm1*(1.0d0-ee*ee)/(r*r)		! (Tangential velocity)^2
	vr = max(v-vt,0.0d0)			! (Radial velocity)^2

	v  = dsqrt(v)
	vt = dsqrt(vt)
	vr = dsqrt(vr)*dsign(1.0d0,vv)		! Sign depends on sign of vv

	lv = lr+datan2d(vt,vr)			! Longitude velocity vector
	lv = dmod(dmod(lv,360.0d0)+360.0d0,360.0d0)

	dv = 0.0d0				! Latitude velocity vector
						! Convert from orbital plane to ecliptic
	if (ii .ne. 0) then			! Non-zero orbit inclination
	    call rotated( -90.0d0, -ii, 90.0d0-OO, lr, dr)
	    call rotated( -90.0d0, -ii, 90.0d0-OO, lv, dv)
	end if

	g = (SpD/AU)*dsqrt(mu/(AU*aa))		! Converts velocities to AU/day

	Position(1) = lr
	Position(2) = dr
	Position(3) = aa*r

	Velocity(1) = lv
	Velocity(2) = dv
	Velocity(3) = g*v
	Velocity(4) = g*vt
	Velocity(5) = g*vr

	return
	end