;+
; NAME:
;   oclc_fwd_n.pro
;
; PURPOSE:
;   Calculates the number density, n, from the
;   radius, pressure, temperature, temperature derivative with respect
;   to radius, refractivity at STP, molecular weight and mass of the planet.
;
; DESCRIPTION:
;
; CALLING SEQUENCE:
;   n = oclc_fwd_n(p, t)
;
; INPUTS: 
;   **** All inputs in cgs units ********
;
;   p - An array (corresponding to r) of the pressure, 
;       in microbar.
;   t - An array (corresponding to r) or scalar of the temperature, 
;       in Kelvin.
;
; OUTPUTS:
;   n - number density
;
; 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
;   nustp = 2.98e-4  ; N2
;   p0 = 1.0d  ; reference pressure in microbar
;   mp = 1.3d25 ; reference planet mass in g
;   p = oclc_fwd_hydrostatic(r,t,mu,r0,p0,mp) ;
;   n = oclc_fwd_n(p, t)

;   ; ---- exact calculation, see Elliot and Young 1992
;   lam  =  (!phys.g*mp*mu*!phys.m_u/(!phys.k*t*r )) ;energy ratio
;   lam0  = (!phys.g*mp*mu*!phys.m_u/(!phys.k*t*r0)) 
;   p2 = p0 * exp( lam - lam0 )
;   n2 = (p2 / (!phys.k * t))
;   
;   plot, r/1e5, (n-n2)/n2
;   print, minmax( (n-n2)/n2 )  ;-5.2533787e-14   9.7999388e-14 
;   print, mean( (n-n2)/n2 )    ;3.4056848e-14
;
;
; 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
;   mu = 28.01d ; scalar molecular weight
;   nustp = 2.98e-4  ; N2
;   p0 = 1.0d  ; reference pressure in microbar
;   mp = 1.3d25 ; reference planet mass in g
;   p = oclc_fwd_hydrostatic(r,t,mu,r0,p0,mp) ; the call.
;   n = oclc_fwd_n(p, t)

;   ; ---- exact calculation, see Elliot and Young 1992
;   lam  =  (!phys.g*mp*mu*!phys.m_u/(!phys.k*t *r )) ;energy ratio
;   lam0  = (!phys.g*mp*mu*!phys.m_u/(!phys.k*t0*r0)) 
;   p2 = p0 * exp( (lam - lam0)/(1 + b) )
;   nl = 2.68719e19    
;   n2 = (p2 / (!phys.k * t))
;   
;   plot, r/1e5, (n-n2)/n2
;   print, minmax( (n-n2)/n2 )  ;  -9.0324651e-15   8.5583455e-15
;   print, mean( (n-n2)/n2 )    ;  -1.9766112e-16
;
; 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_n, p, t

    physconstants


    n = p / (t * !phys.k)
    
    return, n
end