;+
; NAME:
;   oclc_fwd_nu.pro
;
; PURPOSE:
;   Calculates the number density, n, from the
;   pressure, temperature, refractivity at STP
;
; DESCRIPTION:
;
; CALLING SEQUENCE:
;   nu = oclc_fwd_nu(p, t, nustp)
;
; 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.
;   nustp - refractivity at standard temperature and pressure g(unitless)
;
; 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) ;
;   nu = oclc_fwd_nu(p, t, nuSTP)

;   ; ---- 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 )
;   loschmidt = 1.01325e6/(!phys.k * 273.15)
;   nu2 = (p2 / (!phys.k * t)) * (nustp/loschmidt)
;   
;   plot, r/1e5, (nu-nu2)/nu2
;   print, minmax( (nu-nu2)/nu2 )  ;-5.2533787e-14   9.7999388e-14 
;   print, mean( (nu-nu2)/nu2 )    ;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.
;   nu = oclc_fwd_nu(p, t, nuSTP)

;   ; ---- 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) )
;   loschmidt = 1.01325e6/(!phys.k * 273.15)
;   nu2 = (p2 / (!phys.k * t)) * (nustp/loschmidt)
;   
;   plot, r/1e5, (nu-nu2)/nu2
;   print, minmax( (nu-nu2)/nu2 ) 
;   print, mean( (nu-nu2)/nu2 )   
;

; 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_nu, p, t, nustp

    physconstants

    n = oclc_fwd_n(p, t)
    loschmidt = 1.01325e6/(!phys.k * 273.15)
    nu = n * (nustp/loschmidt)
    
    return, nu
end