;+
; NAME:
;   oclc_fwd_hydrostatic.pro
; PURPOSE:
;   Calculates the atmospheric pressure at the tablulated positions
;   given in the r array.
;
; DESCRIPTION:
;        
;
;
; CALLING SEQUENCE:
;   p = oclc_fwd_hydrostatic(r, t, mu, r0, p0, mp) 
;
; 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).
;   r0    - the reference radius level, in cm.
;   p0    - the pressure at the reference radius, in microbar.
;   mp - the planetary mass, in grams.
;
; OUTPUTS:
;   p     - the atmospheric pressure in the atmsophere at the radii in r.
;
; EXAMPLE
;   
;   NUMERICAL RESULTS IN EXAMPLE MAY DEPEND ON ARCHETECTURE
;   Numbers here:
;   ARCH            STRING    'ppc'
;   OS              STRING    'darwin'
;   OS_FAMILY       STRING    'unix'
;   OS_NAME         STRING    'Mac OS X'
;   RELEASE         STRING    '6.2'
;   BUILD_DATE      STRING    'Jun 20 2005'
;
;   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
;   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.
;   ; 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 )
;   plot, r/1e5, (p-p2)/p2
;   print, minmax( (p-p2)/p2 )  ;  -6.8757831e-07    3.1245414e-07
;   print, mean( (p-p2)/p2 )    ; -2.3825349e-07
;
;   EXAMPLE 2 - variable temperature and molecular weight
;   r = dindgen(200)*1d5 + 1200d5  ; array of r, in cm
;   r0 = 1250d5 ; reference r in cm
;   t0 = 104.d   ; temperature at reference altitude
;   b = -2.d    ; exponenet for temperature, chosen so scale height is constant
;   t = t0 * (r/r0)^b ; temperature array
;   mu0 = 28.01 ; exponent for molecular weight
;   a = 0.0
;   mu = mu0 * (r/r0)^a
;   p0 = 1.0d
;   mp = 1.3d25
;   p = oclc_fwd_hydrostatic(r,t,mu,r0,p0,mp)
;   ; exact
;   lam  =  (!phys.g*mp*mu *!phys.m_u/(!phys.k*t*r)) 
;   lam0  = (!phys.g*mp*mu0*!phys.m_u/(!phys.k*t0*r0)) 
;   p2 = p0 * exp( (lam - lam0)/(1+a+b) )
;   plot, r/1e5, (p-p2)/p2
;   print, minmax( (p-p2)/p2 ) ; 
;   print, mean( (p-p2)/p2 )    ; 
;
; REVISON HISTORY:
;   24-Aug-2006 CBO SwRI
;   2006 Sep 20 LAY.  
;     Edited header.
;     Changed from averaging h to averaging 1/h
;- 

function oclc_fwd_hydrostatic, r_in, t_in, mu_in, r0, p0, mp

    physconstants
    k = !phys.k
    
    nz = n_elements(r_in)

 ; check if t_in and mu_in are scalars and make them arrays   

    nt = n_elements(t_in)
    nmu = n_elements(mu_in)
    
    if nt  EQ 1 then tarr  = replicate(t_in,  nz)
    if nmu EQ 1 then muarr = replicate(mu_in, nz)
    
    if nt  EQ nz then tarr  = t_in
    if nmu EQ nz then muarr = mu_in
    
;    if nmu NE 1 OR nmu NE nz then begin
;        print, 'Error: temperature array has different length than radius array'
;    endif
    
;    if nmu NE 1 OR nmu NE nz then begin
;        print, 'Error: molecular weight has different length than radius array'
;    endif
    
 ; sort the input radius, temperature and molecular weight arrays
     s = sort(r_in)
     r = r_in(s)
     t = tarr(s)
     mu = muarr(s) 
 
     gmmkt = !phys.g * mp * (mu*!phys.m_u)/(k*t) ; GMm/kT

     h = dblarr(nz)
; Calculate the average scale height for the tabulated radii
;    for i=0, nz-2 do h[i] = (1/gmmkt[i]+1/gmmkt[i+1])*0.5*r[i+1]^2 
;    h[nz-1] = h[nz-2] + (h[nz-2] - h[nz-3])

     h = r^2/gmmkt
;    There are two ways to handle the scale height:
;    One is to have H = average of the two H's
;     hmid = (h[0:nz-2] + h[1:nz-1])/2
;    Another is to have 1/H = average of the two 1/H's
;     hmid = 2./(1/h[0:nz-2] + 1/h[1:nz-1])
;    Numerical tests (see examples in header) seem to be
;    to average these two approaches
     hmid = (h[0:nz-2] + h[1:nz-1])/4.d + 1./(1/h[0:nz-2] + 1/h[1:nz-1])

; Check that the reference radius is within the radius array

     if r0 LT min(r) OR r0 GT max(r) then begin
         print, 'Error: reference radius is not within the radius array' 
     endif
     
 ; find closest point to reference radius in r array
 ;  the greater than or equal will catch the case where the reference
 ;  radius is in the list of radii.
 ;     belowindex = where( r lt r0, cntbelow )
     aboveindex = where( r ge r0, cntabove )
 
     toprefbin = aboveindex[0]
      
 ; integrate up from the reference radius and pressure: r0, p0
 
     lnpe = replicate(alog(p0), nz)
    
 ; handle the bin with the reference radius in it separately
 ; integrating the lower part of the bin with r0 in it
     
     lnpe[toprefbin] = alog(p0) - $
       (r[toprefbin] - r0)/hmid[toprefbin-1]

     if toprefbin GT 0 then begin
         lnpe[toprefbin - 1] = alog(p0) - $
           (r[toprefbin - 1] - r0)/hmid[toprefbin - 1]
     endif
    
 ; integrating the remaining bins above the reference radius
     for i= toprefbin+1, nz - 1 do begin
         lnpe[i] = lnpe[i-1] - (r[ i ] - r[i-1])/hmid[i-1]
     end
     
 
; integrating the remaining bins below the reference radius
     for i= toprefbin-2, 0, -1 do begin
         
         lnpe[i] = lnpe[i+1] - (r[ i ] - r[i+1])/hmid[i]
         
     end


    
     pe = exp(double(lnpe))
    
 ; Unsort if necessary
 ;    print, N_ELEMENTS(r), N_ELEMENTS(r_in), r[0], r_in[0], N_ELEMENTS(pe)
     if r[0] NE r_in[0] then pe = reverse(pe)
     
     p = pe[0:nz-1]
     return, p
    
end