;+
; NAME:
;   oclc_fwd_dtheta.pro
;
; PURPOSE:
;   Calculates the derivative of the bending angle from the bending angle.
;
; DESCRIPTION:
;
; CALLING SEQUENCE:
;   dtheta = oclc_fwd_dtheta(r, theta) 
;
; INPUTS: 
; **** All inputs in cgs units ********
;
; r - An array of radius values from the center of the planet, in cm.
; theta - An bending angle in radians.
;
; OUTPUTS:
;   dthetadr - the derivative of the bending angle with respect to radius
;
; COMMENTS:
;    calls dydx
;
; EXAMPLE 1 - scalar temperature, molecular weight
;   r = reverse(dindgen(400)*1d5) + 1200d5  ; array of r, in cm  -- OR
;   r = dindgen(400)*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)
;   theta = oclc_fwd_theta(r, dnu)
;   dtheta = oclc_fwd_dtheta(r, theta)
;   dtheta2 = oclc_ey92_dtheta(r0,nu0,lam0,a,b,order,r) ; exact
;   
;   plot, r/1e5, (dtheta-dtheta2) * distobs
;   print, minmax( dtheta-dtheta2 )  * distobs ;
;
; REVISON HISTORY:
;    25-Aug-2006 CBO SwRI -- modified from LAY's lightcurve.pro
;- 
function oclc_fwd_dtheta, r, theta

    nz = n_elements(r)

    if max(theta) ge 0 then stop
    
    lntheta = alog(-theta)
    
    dlnthetadr = fltarr(nz)
                                ; topmost point
    dlnthetadr[0] = dydx(r[0:2], lntheta[0:2], r[0])
    for i=1,nz-2 do begin
        dlnthetadr[i] = dydx(r[i-1:i+1], lntheta[i-1:i+1], r[i])
    end
    dlnthetadr[nz-1] = dydx(r[nz-3:nz-1], lntheta[nz-3:nz-1], r[nz-1])

    dthetadr = theta*dlnthetadr
    
    return, dthetadr
end