;+
; NAME:
;    oclc_fwd_alpha.pro
;
; PURPOSE:
;    Calculates the integral of the refractivity, alpha, from the radius, r, and
;    the refractivity, nu
;
; DESCRIPTION:
;
; CALLING SEQUENCE:
;    alpha = oclc_fwd_alpha(r, nu) 
;
; INPUTS: 
; **** All inputs in cgs units ********
;  r - An array of radius values from the center of the planet, in cm.
;  nu - The derivative of the refractivity with respect to the radius.
;
; OPTIONAL INPUTS
;  rout - array of radii at which to calculate alpha.
;         If not specified, then rout = r
;
; OUTPUTS:
;  alpha  - 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
;   nu = oclc_ey92_nu(r0,nu0,lam0,a,b, r)
;   print, systime()
;   alpha = oclc_fwd_alpha(r, nu)
;   print, systime()
;   alpha2 = oclc_ey92_alpha(r0,nu0,lam0,a,b,order,r) ; exact
;   
;   plot, r/1e5, (alpha-alpha2) * distobs/km
;   print, minmax( alpha-alpha2 )  * distobs/km ;
;
;  REVISON HISTORY:
;   25-Aug-2006 CBO SwRI -- modified from LAY's lightcurve.pro
;- 
function oclc_fwd_alpha, r, nu, rout=rout

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

    if not keyword_set(rout) then begin

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

        alpha[s] = alphasort

    endif else begin

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

        alpha = 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)
            alpha[i] = 2* atmint(rsort,nusort,rout[i],rsort[0], 0.,/exp)

        end

        for ii=0, nhiextrap-1 do begin

            i = ihiextrap[ii]

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

        end

        for ii=0, nlowextrap-1 do begin

            i = ilowextrap[ii]

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

        
    endelse
    
    return, alpha
        
end