;+ ; NAME: ; oclc_fwd_dnu.pro ; ; PURPOSE: ; Calculates the derivative of the refractivity, dnudr, from the ; radius, pressure, temperature, temperature derivative with respect ; to radius, refractivity at STP, molecular weight and mass of the planet. ; ; DESCRIPTION: ; ; CALLING SEQUENCE: ; dnu = oclc_fwd_dnu(r, p, t, mu, nuSTP, 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. ; p - An array (corresponding to r) of the pressure, ; in microbar. ; 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). ; nuSTP - the refractivity of the gas at STP ; mp - the planetary mass, in grams. ; ; OOPTIONAL INPUTS: ; dt - the derivative of temperature with radius ; ; OUTPUTS: ; dnu - derivative of 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) ; the call. ; dnu = oclc_fwd_dnu(r, p, t, mu, nuSTP, mp) ; ; ---- 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 ) ; nl = 2.68719e19 ; nu2 = (p2 / (!phys.k * t)) * (nustp/nl) ; dnu2 = -nu2 * (lam) / r ; ; plot, r/1e5, (dnu-dnu2)/dnu2 ; print, minmax( (dnu-dnu2)/dnu2 ) ; -6.0180608e-15 -4.4277798e-15 ; print, mean( (dnu-dnu2)/dnu2 ) ; -3.9870151e-24 ; ; ; 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. ; dnu = oclc_fwd_dnu(r, p, t, mu, nuSTP, mp) ; ; ---- 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 ; nu2 = (p2 / (!phys.k * t)) * (nustp/nl) ; dnu2 = -nu2 * (lam + b) / r ; ; plot, r/1e5, (dnu-dnu2)/dnu2 ; print, minmax( (dnu-dnu2)/dnu2 ) ; -6.0180608e-15 -4.4277798e-15 ; print, mean( (dnu-dnu2)/dnu2 ) ; -3.9870151e-24 ; ; 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_dnu, r_in, p_in, t_in, mu_in, nuSTP, mp, DT = dt_in physconstants k = !phys.k m_u = !phys.m_u nt = n_elements(t_in) nmu = n_elements(mu_in) nz = n_elements(r_in) ; print, nt, nmu, nz 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 rarr = r_in parr = p_in ; Loschmidt's number (molecule/cm^3 at STP) ; number density of an ideal gas at 0 deg C and 1 atmpshere ; nl = 2.68719e19 nl = 1.01325e6/(!phys.k * 273.15) n = parr/(tarr*k) ; molecules/cm^3 nu = n * nuSTP / nl ; refractivity gravacc = !phys.g * mp / (rarr^2) sclhgt= k*tarr/(muarr*m_u*gravacc) if keyword_set(DT) then begin dt = dt_in endif else begin ; calculate dt from t and r nz = n_elements(r_in) dt = dblarr(nz) ; topmost point dt[0] = dydx(rarr[0:2], tarr[0:2], rarr[0]) ; other points for i=1,nz-2 do begin dt[i] = dydx(rarr[i-1:i+1], tarr[i-1:i+1], rarr[i]) end dt[nz-1] = dydx(rarr[nz-3:nz-1], tarr[nz-3:nz-1], rarr[nz-1]) ; since we don't have a point past nz - 1, assume the temperature ; derivative is constant at the end ; dt[ nz-1 ] = dt[ nz-2 ] endelse dnu = -nu * (1/sclhgt + dt/tarr) return, dnu end ;+ ; NAME: ; oclc_fwd_dnuTEST.pro ; ; PURPOSE: ; Test the function oclc_fwd_dnu ; ; DESCRIPTION: ; ; CALLING SEQUENCE: ; oclc_fwd_dnuTEST ; ; INPUTS: ; ; OUTPUTS: ; ; REVISON HISTORY: ; 25-Aug-2006 CBO SwRI ;- pro oclc_fwd_dnuTEST dir = '/Users/colkin/projects/P384.2/TestMaterial/' table = mrdfits(dir + 'testphi.fits',1,h) r_in = table[0:900].r*1.d5 ; cm t_in = table[0:900].t ; K mu_in = 28.01 r0 = table[400].r*1.d5 ; cm p0 = table[400].p ; microbar mp = 1.30497e+25 ; grams km = 1.e-5 p_in = oclc_fwd_hydrostatic(r_in, t_in, mu_in, r0, p0, mp) nuSTP = 3.0000e-4 physconstants k = !phys.k m_u = !phys.m_u nt = n_elements(t_in) nmu = n_elements(mu_in) nz = n_elements(r_in) ; print, nt, nmu, nz 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 rarr = r_in parr = p_in ; Loschmidt's number (molecule/cm^3 at STP) ; number density of an ideal gas at 0 deg C and 1 atmpshere nl = 2.68719e19 n = parr/(tarr*k) ; molecules/cm^3 nu = n * nuSTP / nl ; refractivity gravacc = !phys.g * mp / (rarr^2) sclhgt= k*tarr/(muarr*m_u*gravacc) plot, rarr*km, sclhgt*km oplot, table[0:999].r, table[0:999].h, color='ff0000'xl plot, rarr*km, (sclhgt*km - table[0:999].h)/table[0:999].h ; agrees to 10^-6 in fractional error ; calculate dt from t and r nz = n_elements(r_in) dt = fltarr(nz) ; topmost point dt[0] = dydx(rarr[0:2], tarr[0:2], rarr[0]) ; other points for i=1,nz-2 do dt[i] = dydx(rarr[i-1:i+1], tarr[i-1:i+1], rarr[i]) dt[nz-1] = dydx(rarr[nz-3:nz-1], tarr[nz-3:nz-1], rarr[nz-1]) ; since we don't have a point past nz - 1, assume the temperature derivative is ; constant at the end dt[ nz-1 ] = dt[ nz-2 ] dnudr = -nu * (1/sclhgt + dt/tarr) plot, rarr*km, dt/km oplot, table[0:900].r, table[0:900].dt, color='ff0000'xl plot, rarr*km, (dt/km - table[0:999].dt)/table[0:999].dt ; good to 0.05 in fractional error. Ringing is due to the small step size ; in r (happens in dydx) end