;+
; NAME:
;    oclc_fwd_theta.pro
;
; PURPOSE:
;    Calculates the bending angle, theta, from the radius, r, and
;    the derivative of the refractivity, dnudr
;
; DESCRIPTION:
;    See Chamberlain and Elliot (1997) equation 2
;
; CALLING SEQUENCE:
;    theta = oclc_fwd_theta(r, dnudr) 
;
; INPUTS: 
; **** All inputs in cgs units ********
;  r - An array of radius values from the center of the planet, in cm.
;  dnu - The derivative of the refractivity with respect to the radius.
;
; OUTPUTS:
;  theta  - the bending angle in radians.
;
; COMMENTS;
;   calls function atmint
;
; EXAMPLE 1 - scalar temperature, molecular weight
;   r = reverse(dindgen(1000)*1d5) + 1200d5  ; array of r, in cm  -- OR
;   r = dindgen(1000)*1d5 + 1200d5  ; array of r, in cm
;   km = 1e5
;   distobs = 30.*1.496e8*1.d5 ; 30 AU in cm
;   r0 = 1250d5 ; reference r in cm
;   nu0 = 2d-9
;   lam0 = 60.d
;   a = 0.d
;   b = 0.d
;   order = 4
;   dnu = oclc_ey92_dnu(r0,nu0,lam0,a,b, r)
;   print, systime()
;   theta = oclc_fwd_theta(r, dnu)
;   print, systime()
;   theta2 = oclc_ey92_theta(r0,nu0,lam0,a,b,order,r) ; exact
;   
;   plot, r/1e5, (theta-theta2) * distobs/km
;   print, minmax( theta-theta2 )  * distobs/km ;
;
;  REVISON HISTORY:
;   25-Aug-2006 CBO SwRI -- modified from LAY's lightcurve.pro
;- 
function oclc_fwd_theta, r, dnu

    nz = n_elements(r)
    
    s = reverse(sort(r))

    rsort = r[s]

    invrdnu = dnu[s]/rsort             ; (1/r dnu/dr)

    thetasort = dblarr(nz)
    theta = dblarr(nz)

    for i=1,nz-1 do begin
    
    ; Equation 2 from Chamberlain and Elliot (1997)
        thetasort[i] = -2*rsort[i] * atmint(rsort,-invrdnu,rsort[i],rsort[0], 0.,/exp)
   
    end
   
    for i=1,0,-1 do begin
        thetasort[i] = thetasort[i+1] * $
          (thetasort[i+2]/thetasort[i+1])^( (rsort[i]-rsort[i+1])/(rsort[i+1]-rsort[i+2]))
    end

    theta[s] = thetasort

    return, theta
    
end