;+
; 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.
;
; OPTIONAL INPUTS
;  rout - array of radii at which to calculate theta.
;         If not specified, then rout = r
;
; 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, rout=rout

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

    if not keyword_set(rout) then begin

        
        nz = n_elements(r)
    
        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+2]-rsort[i+1]))
        end

        theta[s] = thetasort

    endif else begin

        nz = n_elements(r)
        nzout = n_elements(rout)

        theta = dblarr(nzout)

        ; there are three classes: 
        ; extrapolate to lower r, interpolate, extrapolate to higher r
        ilowextrap = where(rout lt rsort[nz-1], nlowextrap)
        iinterp = where(rout ge rsort[nz-1] and rout lt rsort[1], ninterp)
        ihiextrap = where(rout ge rsort[1], nhiextrap)

        ss = reverse(sort(r[iinterp]))


        for ii=0, ninterp-1 do begin

            i = iinterp[ii]

            ; Equation 2 from Chamberlain and Elliot (1997)
            theta[i] = -2*rout[i] * atmint(rsort,-invrdnu,rout[i],rsort[0], 0.,/exp)

        end

        for ii=0, nhiextrap-1 do begin

            i = ihiextrap[ii]

            theta[i] = theta[ss[3]] * $
              (theta[ss[2]]/theta[ss[3]])^( (rout[i]-rsort[3])/(rsort[2]-rsort[3]))

        end

        for ii=0, nlowextrap-1 do begin

            i = ilowextrap[ii]

            theta[i] = theta[ss[nz-1]] * $
              (theta[ss[nz-2]]/theta[ss[nz-1]])^ $
              ( (rout[i]-rsort[nz-1]) / (rsort[nz-2]-rsort[nz-1]) )
        end

        
    endelse
    
    return, theta
        
end