pro exofast_rossiter, x, y, u1, p, vsini, lambda, deltarv, z=z
;+
; NAME:
;   EXOFAST_ROSSITER
;
; PURPOSE: 
;   This routine implements the Ohta et al 2005 approximation for the
;   Rossiter-Mclaughlin effect, which assumes linear limb darkening
;   and is accurate to ~1 m/s for slowly rotating stars (see Hirano 2011).
;
; CALLING SEQUENCE:
;   exofast_rossiter, x, y, u1, p, vsini, lambda, deltarv, z=z
;
; INPUTS:
;    x      - array of star-planet distances in the x direction
;             (stellar radii).
;    y      - array or scalar of star-planet distance(s) in the y
;             direction. ie, the impact parameter (stellar radii)
;    u1     - linear limb-darkening coefficient -- epsilon in paper
;    p      - occulting star size (stellar radii) -- gamma in paper
;    vsini  - projected stellar surface velocity (same as output units)
;    lambda - projected spin-orbit angle (radians)
;    
; OPTIONAL INPUTS:
;    z      - array of star-planet distances in the z direction. 
;               -- Positive => in front of star (RM effect)
;               -- Negative => behind star (no effect)
;               -- assumed positive if not specified
; OUTPUTS:
;    deltarv - Approximate Change in RV at each z due to the Rossiter
;              Mclaughlin effect (units same as vsini) 
;                     *** accurate to ~1 m/s ***
;
; MODIFICATION HISTORY
;  2009/04/30 -- Jason Eastman (Ohio State University)
;     Loose translation from Scott Gaudi's (OSU) fortran code
;  2012/08 -- Jason Easmtan (LCOGT)
;     Make sign convention ahere to standard (now consistent with EXOFAST_GETB).
;-

nx = n_elements(x)
if n_elements(z) eq 0 then z = dblarr(nx)+1d0
deltarv = dblarr(nx)
rho = sqrt(x^2+y^2)

xp = x*cos(lambda)-y*sin(lambda)
;yp = x*sin(lambda)+y*cos(lambda)
;theta = dindgen(1000)/1000d0*2d0*!dpi
;plot, cos(theta),sin(theta),/iso, xrange=[-2,2],yrange=[-2,2]
;oplot, xp, yp

etap = rho-1.d0         ;; equation 28

;; transit phase
transit = where(rho le (1d0-p) and z gt 0d0)
if transit[0] ne -1 then begin
    w1 = 0d0                    ;; equation 43
    w2 = sqrt(1-rho[transit]^2) ;; equation 44
    ;; equation 40
    deltarv[transit]=-vsini*xp[transit]*p^2*(1d0-u1*(1d0-w2))/$
      (1d0-p^2-u1*(1d0/3d0-p^2*(1d0-w1)))
endif

;; ingress/egress phase
inegress = where(abs(etap) lt p and z gt 0)
if inegress[0] ne -1 then begin
    ;; equation 33
    x0 = 1.d0 - (p^2-etap[inegress]^2)/(2d0*(1d0+etap[inegress])) ;; equation 33
    zeta = 1.d0 + etap[inegress] - x0 ;; equation 34
    z0 = sqrt(1-x0^2)       ;; equation 33
    xc=x0+(zeta-p)/2d0      ;; equation 50
    w2 = sqrt(2d0*p-p^2)    ;; equation 44
    w3 = 0.d0               ;; equation 48

    ;; equation 51
    g1 = (1d0-xc^2)*asin(sqrt((p^2-(xc-1-etap[inegress])^2)/(1.d0-xc^2))) + $
      sqrt((p^2-(xc-1.d0-etap[inegress])^2)*(1.d0-xc^2-p^2+$
                                             (xc-1.d0-etap[inegress])^2))
    g2 = sqrt(2.d0*p^3*(1.d0-p)) + p*(2.d0-p)*asin(sqrt(p/(2.d0-p)))

    w4 = (!dpi/2.d0)*p*(p-zeta)*xc*g1/g2*w2 ;; equation 49

    ;; equation 45
    num=(1.d0-u1)*(-z0*zeta+p^2*acos(zeta/p))+$
      u1*w4/(1.d0+etap[inegress])

    den1 = !dpi*(1.d0-1.d0/3.d0*u1)-u1*w3
    den2 = -(1.d0-u1)*(asin(z0)-$
               (1.d0+etap[inegress])*z0+p^2*acos(zeta/p))
    deltarv[inegress]=-vsini*xp[inegress]*num/(den1+den2)
endif

end