;+
; NAME:
;   EXOFAST_GETB
;
; PURPOSE: 
;   This function returns the impact parameter as a function of BJD
;   given orbital elements of the planet, assuming a keplerian orbit.
;   Optionally, the 3 space coordinates can be returned.
; 
; CALLING SEQUENCE:
;   result = getb(jds, i=i,a=a,tperiastron=tperiastron,period=p,$
;                 e=e,omega=omega, x=x, y=y, z=z)
;
; INPUTS:
;   bjd         - Barycentric Julians dates for desired impact
;                 parameters (ideally in the target's
;                 barycentric frame; see BJD2TARGET)
;   i           - inclination of orbit (radians)
;   a           - semi-major axis (in units of R_*)
;   tperiastron - periastron passage time (BJD)
;   Period      - Period of orbit (days)
;
; OPTIONAL INPUTS:  
;   e           - eccentricity of orbit (0 if not specified)
;   omega       - argument of periastron of the star's orbit, in radians
;                 -- omega_* is typically quoted from RV
;                 -- required if e is specified
;                 -- assumed to be pi/2 if e not specified
;                 -- omega_* = omega_planet + !dpi
;   lonascnode  - The Longitude of the ascending node
;                 (radians). Assumed to be !dpi if not specified.
;
; OUTPUTS:
;    result     - the impact parameter as a function of BJD, in units
;                 of R_*.
;
; OPTIONAL OUTPUTS:
;    x,y,z - Arrays of the cartesian coordinates of the planet at
;            each BJD, in the units of R_*. 
;              +X is right
;              +Y is up
;              +Z is out of the page (primary transit)
;
; MODIFICATION HISTORY 
;  2009/05/01 -- Jason Eastman (Ohio State University)
;-

function exofast_getb, bjd, i=i, a=a, tperiastron=tperiastron, Period=P, $
                       e=e, omega=omega, x=x, y=y, z=z, $
                       lonascnode=lonascnode

;; calculate the mean anomaly corresponding to each observed time
meananom = (2.d0*!dpi*(1.d0 + (bjd - Tperiastron)/P)) mod (2.d0*!dpi) 

;; if eccentricity is given, integrate the orbit
if keyword_set(e) then begin
    if n_elements(omega) eq 0 then message, $
      'ERROR: omega must be specified if e is specified'
    eccanom = exofast_keplereq(meananom, e)
    trueanom = 2d0*atan(sqrt((1d0 + e)/(1d0 - e))*tan(0.5d0*eccanom))
endif else begin
    trueanom = meananom
    e=0.d0
endelse

;; calculate the corresponding (x,y) coordinates of planet
r = a*(1d0-e^2)/(1d0+e*cos(trueanom))

;; as seen from observer
x = -r*cos(trueanom + omega)
tmp = r*sin(trueanom + omega)
y =  -tmp*cos(i)
z =  tmp*sin(i)

;; Rotate by the Longitude of Ascending Node
;; For transits, it is not constrained, so we assume Omega=!dpi)
if keyword_set(lonascnode) then begin
    xold = x & yold = y
    x = -xold*cos(lonascnode) + yold*sin(lonascnode)
    y = -xold*sin(lonascnode) - yold*cos(lonascnode)
endif

b = sqrt(x^2 + y^2)

return, b

end