function orb_pos, semi, ecc, asc_node, peri_long, inc, mean_anom
;+
; NAME:
;	orb_pos
;
; PURPOSE:
;	This routine returns the 3-D position vector of an orbital body, given the orbital 
;	elements. Uses equations worked out in Murray & Dermott (1999) _Solar System Dynamics_,
;	ch. 2. Available with no warranties whatsoever at
;	http://www.lpl.arizona.edu/~bjackson/idl_code/index.html.
;
; CATEGORY:
;	Astrophysics.
;
; CALLING SEQUENCE:
;
;	Result = orb_pos(semi, ecc, asc_node, peri_long, inc, mean_anom)
;
; INPUTS:
;	semi:		orbital semi-major axis
;	ecc:		eccentricity
;	asc_node:	longitude of ascending node, in degrees
;	peri_long: 	longitude of pericenter, in degrees
;	inc: 		orbital inclination, in degrees
;	mean_anom: 	orbital mean anomaly
;
; OUTPUTS:
;	3 x n_elements(mean_anom) array with x, y, z position of orbital body
;
; RESTRICTIONS:
;	Code doesn't do any checking of parameters, so do them yourself. The code requires a 
;	companion routine keplereq.pro (written by Marc Buie and Joern Wilms, among others) 
;	available here: http://www.lpl.arizona.edu/~bjackson/idl_code/index.html.
;
; EXAMPLE:
;	semi = 1.0
;	ecc = 0.0 ;orbital eccentricity
;	asc_node = 0. ;longitude of planetary ascending node
;	peri_long = 0. ;longitude of planetary pericenter
;	inc = 83.1 ;orbital inclination in degrees
;	mean_anom = 2.*!pi*dindgen(101)
;	r = orb_pos(semi, ecc, asc_node, peri_long, inc, mean_anom)
;
; MODIFICATION HISTORY:
; 	Written by:	Brian Jackson (decaelus@gmail.com), 2011 July 25.
;-

  ;calculate true anomaly as a function of time, i.e. solve Kepler's equation 
  ecc_anom = keplereq(mean_anom, ecc)
  cos_f = (cos(ecc_anom)-ecc)/(1.-ecc*cos(ecc_anom))
  sin_f = (sqrt(1.-ecc^2.)*sin(ecc_anom))/(1.-ecc*cos(ecc_anom))

  ;cos(inclination)
  cos_inc = cos(inc)
  sin_inc = sin(inc)

  ;cos(longitude of ascending node)
  cos_asc = cos(asc_node)
  sin_asc = sin(asc_node)

  ;longitude of pericenter
  cos_peri = cos(peri_long)
  sin_peri = sin(peri_long)

  ;cos(peri_long + f)
  cos_wf = cos_peri*cos_f-sin_peri*sin_f
  sin_wf = sin_peri*cos_f+cos_peri*sin_f

  r = dblarr(3, n_elements(mean_anom))
  r[0,*] = semi*(1.-ecc^2.)/(1.+ecc*cos_f)*$
    (cos_asc*cos_wf-sin_asc*sin_wf*cos_inc)
  r[1,*] = semi*(1.-ecc^2.)/(1.+ecc*cos_f)*$
    (sin_asc*cos_wf+cos_asc*sin_wf*cos_inc)
  r[2,*] = semi*(1.-ecc^2.)/(1.+ecc*cos_f)*sin_wf*sin_inc

  return, r
end;function orb_pos