;+
; NAME:
;   oclc_fwd_hn.pro
;
; PURPOSE:
;   Calculates the number density scale height, hn,
;
; DESCRIPTION:
;
; CALLING SEQUENCE:
;   hn = oclc_fwd_hn(r, t, mu, mp, dt=dt) 
;
; INPUTS: 
;   **** All inputs in cgs units ********
;
;   r - An array of radius values from the center of the planet, in cm.
;       The radius array could be in ascending or descending order.
;   t - An array (corresponding to r) or scalar of the temperature, 
;       in Kelvin.
;   mu - An array (corresponding to r) or scalar of the molecular weight
;         (unitless).
;   mp - the planetary mass, in grams.
;
; KEYWORD INPUTS AND OUTPUTS
;   dt - the temperature derivative at r, in Kelvin/cm
;        If dt is not passed, oclc_fwd_hn calculates dt from r and t.
;        If dt is passed but the variable is undefined, then
;        the calculated dt gets placed into the passed variable
;        (see dtarr in example 2)
; 
; OUTPUTS:
;   hn - number density scale height in cm
;
; EXAMPLE 1 - scalar temperature, molecular weight
;   r = dindgen(200)*1d5 + 1200d5  ; array of r, in cm
;   r0 = 1250d5 ; reference r in cm
;   t = 104.d  ; scalar temperature in K
;   mu = 28.01d ; scalar molecular weight
;   mp = 1.3d25 ; reference planet mass in 
;   h = oclc_fwd_h(r, t, mu, mp)
;   hn = oclc_fwd_hn(r, t, mu, mp
;   hnarr = oclc_fwd_hn(r, replicate(t,200), mu, mp)
;   hndt = oclc_fwd_hn(r, replicate(t,200), mu, mp, dt=replicate(0.,200))

;   ; ---- exact calculation, see Elliot and Young 1992
;   lam  =  (!phys.g*mp*mu*!phys.m_u/(!phys.k*t*r )) ;energy ratio
;   hn2 = r/lam
;   
;   plot, r/1e5, (hn-hn2)/hn2
;   print, minmax( (hn-hn2)/hn2 )  ;
;   print, mean( (hn-hn2)/hn2 )    ;
;
;   plot, r/1e5, (hnarr-hn2)/hn2
;   print, minmax( (hnarr-hn2)/hn2 )  ;
;   print, mean( (hnarr-hn2)/hn2 )    ;
;
;   plot, r/1e5, (hndt-hn2)/hn2
;   print, minmax( (hndt-hn2)/hn2 )  ;
;   print, mean( (hndt-hn2)/hn2 )    ;
;
;
; EXAMPLE 2 - variable temperature
;   r = dindgen(200)*1d5 + 1200d5  ; array of r, in cm
;   r0 = 1250d5 ; reference r in cm
;   b = -2.d
;   t0 = 104.d 
;   t = t0 * (r/r0)^b  ; scalar temperature in K
;   dt = b*t/r
;   mu = 28.01d ; scalar molecular weight
;   mp = 1.3d25 ; reference planet mass in g
;   h = oclc_fwd_h(r, t, mu, mp)
;   delvarx, dtarr
;   hnarr = oclc_fwd_hn(r, t, mu, mp, dt=dtarr)
;   hndt = oclc_fwd_hn(r, t, mu, mp, dt=dt)

;   ; ---- exact calculation, see Elliot and Young 1992
;   lam  =  (!phys.g*mp*mu*!phys.m_u/(!phys.k*t *r )) ;energy ratio
;   hn2 = r/(lam + b)
;   
;   plot, r/1e5, (hnarr-hn2)/hn2
;   print, minmax( (hnarr-hn2)/hn2 )  ;
;   print, mean( (hnarr-hn2)/hn2 )    ;
;
;   plot, r/1e5, (hndt-hn2)/hn2
;   print, minmax( (hndt-hn2)/hn2 )  ;
;   print, mean( (hndt-hn2)/hn2 )    ;
;
; REVISON HISTORY:
; 25-Aug-2006 CBO SwRI -- modified from LAY's lightcurve.pro
; 23-Sep-2006 LAY SwRI -- changed some variable names, added example
; 
;- 
function oclc_fwd_hn, r, t, mu, mp, dt=dt

    h = oclc_fwd_h(r, t, mu, mp) ; pressure scale height

    ; constant temperature: scalar t and no dt passed
    if n_elements(t) eq 1 and n_elements(dt) eq 0 then begin
        hn = h
        return, hn
    endif

    ; explicit derivative: dt passed
    ; t is presumably array (scalar t legal but senseless)
    if n_elements(dt) ne 0 then begin
        hn = 1.d/(1.d/h + dt/t)
        return, hn
    endif

    ; implicit derivative: array t, no dt
    if n_elements(t) gt 1 and n_elements(dt) eq 0 then begin

        nz = n_elements(r)
        dt = dblarr(nz)
        ; topmost point
        dt[0] = dydx(r[0:2], t[0:2], r[0])
        ; other points
        for i=1,nz-2 do begin
            dt[i] = dydx(r[i-1:i+1], t[i-1:i+1], r[i])
        end
        dt[nz-1] = dydx(r[nz-3:nz-1], t[nz-3:nz-1], r[nz-1])

        hn = 1.d/(1.d/h + dt/t)
        return, hn
    endif

end