;+ ; 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 ; 2007 Sep 29 LAY ; when toprefbin = 0, use hmid[0] for lnpe[toprefbin] ;- 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 if toprefbin GT 1 then begin lnpe[toprefbin] = alog(p0) - $ (r[toprefbin] - r0)/hmid[toprefbin-1] endif else begin lnpe[toprefbin] = alog(p0) - $ (r[toprefbin] - r0)/hmid[toprefbin] endelse 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