C+
C NAME:
C	KeplerOrbit
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 KeplerOrbit(JD,m1i,m2i,OrbitElement,Position,Velocity)
C INPUTS:
C	JD		double precision
C					time difference (days) with epoch t0
C	m1		real		primary mass (in solar masses)
C					if m1 <=0 then m1 = 1 (solar mass) is used
C	m2		real		secondary mass (in solar masses)
C					if m2 < 0 then m2 = 0 is used. 
C	OrbitElement(5)	real		orbital elements:
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)	real		position vector relative to m1:
C					1: ecliptic longitude (deg),
C					2: ecliptic latitude (deg),
C					3: radial distance (AU)
C	Velocity(5)	real		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 CALLS:
C	nrZBrent, rotate
C EXTERNAL:
	external	EqKepler
C INCLUDE:
	include		'math.h'
	include		'sun.h'
C COMMON BLOCKS:
C	common /KEPLER/ M,ee
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-
	    double precision	JD
	    real	m1i
	    real	m2i
	    real	OrbitElement(6)
	    real	Position(3)
	    real	Velocity(5)

	real		m1
	real		m2
	real		mu
	real		nn
	real		M
	real		lr
	real		lv
	real		ii
	real		ll

	common		/KEPLER/ M,ee
	real            nrZBrent

	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					! Longitude of perihelion
	if (ii .ne. 0.0) ww = wb-OO		! Argument of perihelion (angle in orbital plane)

	m1 = m1i
	if (m1 .le. 0.0) m1 = 1.0		! Default: solar mass

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

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

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

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

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

	M = nn*JD+ll-wb				! Mean anomaly (degrees)
	M = mod( mod(M,360.0)+360.0, 360.0 )	! In [0,360]
						! Eccentric anomaly (degrees)
	E = nrZBrent(EqKepler,0.0,360.0,0.00001)! Solve Kepler's equation

	if (ee .lt. 1.0) then			! Elliptical orbit
	    r = 1.0-ee*cosd(E)			! Distance to m1 in units of aa
	    vv = atan2d(sqrt(1.0-ee*ee)*sind(E),cosd(E)-ee)! True anomaly (degrees)
	else					! Hyperbolic orbit
	    E = E*MATH__RPD
	    r = ee*cosh(E)-1.0
	    vv = atan2d(sqrt(ee*ee-1.0)*sinh(E),ee-cosh(E))! True anomaly (degrees)
	end if

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

	dr = 0.0				! Latitude radius vector

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

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

	lv = lr+atan2d(vt,vr)			! Longitude velocity vector
	lv = mod( mod(lv,360.0)+360.0, 360.0 )

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

	g = (SpD/AU)*sqrt(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