This page was created by the IDL library routine
mk_html_help. For more information on
this routine, refer to the IDL Online Help Navigator
or type:
? mk_html_help
at the IDL command line prompt.
Last modified: Thu Sep 28 22:47:38 2006.
NAME:
atmint
PURPOSE: (one line)
Integrate a function slantwise through an atmosphere
DESCRIPTION:
CATEGORY:
CALLING SEQUENCE:
integral= atmint (r, f, r0, r1, mu, NOEXTRAP=noextrap, EXPON=expon)
INPUTS:
int - an array of radii. Does NOT need to be sorted
inf - a function tabulated at the radii r
r0 - the lower bound
r1 - the upper bound
mu - the zenith angle (at r=r0) along which to perform the integration
OPTIONAL INPUT PARAMETERS:
KEYWORD PARAMETERS:
NOEXTRAP - if set, the routine will not extrapolate past
the ends of the array
EXPONENTIAL - if set, the array f is assumed to be exponential in r
for the purpose of interpolation (currently the
only implemented mode).
OUTPUTS:
integral - the integral
KEYWORD OUTPUTS:
COMMON BLOCKS:
SIDE EFFECTS:
RESTRICTIONS:
PROCEDURE:
MODIFICATION HISTORY:
(See atmint.pro)
NAME: oclc_fwd_AAREADME PURPOSE Print the file 'oclc_fwd_AAREADME.pro' CATEGORY: Occultation lightcurve (oclc) CALLING SEQUENCE: oclc_fwd_AAREADME PROCEDURE oclc_fwd_AAREADME These routines implement forward models The routines are called oclc_fwd_*. oclc stands for OCcultation LightCurve. fwd stands for "forward." Other occultation lightcurves can be implemented, and will be prefixed oclc__.pro. --------------------------------------------------------------------------- VARIABLES --------------------------------------------------------------------------- distobs D in EY92 Distance between target and observer cm dnu d nu / dr in EY92 derivative of refractivity with respect to planet radius 1/cm dt Derivative of the temperature with respect to planet radius K/cm dtheta d theta/d r in EY92 Derivative of the bending angle with respect to planet radius, positive for an isothermal atmosphere. radian/cm h H_p in EY92 Pressure scale height cm hn H_n in EY92 Number density scale height cm kappa linear absorption coefficient kappa = kappa1 exp( - (r-r1)/(htau1 * (r/r1) ) ) EY92 3.23 1/cm (or cm^2/cm^3, cross section per particle * particle/volume) mp mass of planet (or dwanet) gm n number density molecule/cm^3 nu Refractivity. Unitless phi normalized stellar flux. unitless. p pressure microbar p0 pressure at reference radius, r0 microbar r planet radius (distance from planet center). cm r0 reference planet radius (distance from planet center). cm y shadow radius (distance from shadow center), negative for far-limb contribution. cm. t temperature K t0 temperature at reference radius, r0 K theta bending angle, negative for an isothermal atmosphere. radians. tauobs line-of-sight optical depth unitless. --------------------------------------------------------------------------- FUNCTIONS --------------------------------------------------------------------------- p = oclc_fwd_hydrostatic(r, t, mu, r0, p0, mp) h = oclc_fwd_h(r, t, mu, mp) n = oclc_fwd_n(p, t) hn = oclc_fwd_hn(r, t, mu, mp, dt=dt) nu = oclc_fwd_nu(p, t, nustp) dnu = oclc_fwd_dnu(r, p, t, mu, nuSTP, mp) theta = oclc_fwd_theta(r, dnu) dtheta = oclc_fwd_dtheta(r, theta) y = oclc_fwd_y(r, theta, distobs) phi = oclc_fwd_phiref(r, theta, dtheta, distobs)
(See oclc_fwd_AAREADME.pro)
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)
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 )
loschmidt = 1.01325e6/(!phys.k * 273.15)
nu2 = (p2 / (!phys.k * t)) * (nustp/loschmidt)
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
dt = b*t/r
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)
hnarr = oclc_fwd_hn(r, t, mu, mp)
hn = oclc_fwd_hn(r, t, mu, mp, dt=dt)
dnuarr = oclc_fwd_dnu(r, p, t, mu, nuSTP, mp)
dnu = oclc_fwd_dnu(r, p, t, mu, nuSTP, mp, dt=dt)
; ---- 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)
hn2 = r/(lam + b)
dnu2 = -nu2 * (lam + b) / r
plot, r/1e5, (dnuarr-dnu2)/dnu2
print, minmax( (dnuarr-dnu2)/dnu2 ) ; -1.4146718e-07 2.8278174e-07
print, mean( (dnuarr-dnu2)/dnu2 ) ; -1.0877865e-07
plot, r/1e5, (dnu-dnu2)/dnu2
print, minmax( (dnu-dnu2)/dnu2 ) ;-9.1506314e-15 8.6675652e-15
print, mean( (dnu-dnu2)/dnu2 ) ;-2.4523038e-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
(See oclc_fwd_dnu.pro)
NAME:
oclc_fwd_dtheta.pro
PURPOSE:
Calculates the derivative of the bending angle from the bending angle.
DESCRIPTION:
CALLING SEQUENCE:
dtheta = oclc_fwd_dtheta(r, theta)
INPUTS:
**** All inputs in cgs units ********
r - An array of radius values from the center of the planet, in cm.
theta - An bending angle in radians.
OUTPUTS:
dthetadr - the derivative of the bending angle with respect to radius
COMMENTS:
calls dydx
EXAMPLE 1 - scalar temperature, molecular weight
r = reverse(dindgen(400)*1d5) + 1200d5 ; array of r, in cm -- OR
r = dindgen(400)*1d5 + 1200d5 ; array of r, in cm
km = 1e5
distobs = 30.*1.496e8*1.d5 ; 30 AU in cm
r0 = 1250d5 ; reference r in cm
nu0 = 2d-9
lam0 = 60.d
a = 0.d
b = 0.d
order = 4
dnu = oclc_ey92_dnu(r0,nu0,lam0,a,b, r)
theta = oclc_fwd_theta(r, dnu)
dtheta = oclc_fwd_dtheta(r, theta)
dtheta2 = oclc_ey92_dtheta(r0,nu0,lam0,a,b,order,r) ; exact
plot, r/1e5, (dtheta-dtheta2) * distobs
print, minmax( dtheta-dtheta2 ) * distobs ;
REVISON HISTORY:
25-Aug-2006 CBO SwRI -- modified from LAY's lightcurve.pro
(See oclc_fwd_dtheta.pro)
NAME:
oclc_fwd_h.pro
PURPOSE:
Calculates the pressure scale height, h,
DESCRIPTION:
CALLING SEQUENCE:
h = oclc_fwd_h(r, t, mu, 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).
mp - the planetary mass, in grams.
OUTPUTS:
h - pressure scale height in cm
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
mp = 1.3d25 ; reference planet mass in g
h = oclc_fwd_h(r, t, mu, mp)
; ---- exact calculation, see Elliot and Young 1992
lam = (!phys.g*mp*mu*!phys.m_u/(!phys.k*t*r )) ;energy ratio
h2 = r/lam
plot, r/1e5, (h-h2)/h2
print, minmax( (h-h2)/h2 ) ;
print, mean( (h-h2)/h2 ) ;
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
mp = 1.3d25 ; reference planet mass in g
h = oclc_fwd_h(r, t, mu, mp)
; ---- exact calculation, see Elliot and Young 1992
lam = (!phys.g*mp*mu*!phys.m_u/(!phys.k*t *r )) ;energy ratio
h2 = r/lam
plot, r/1e5, (h-h2)/h2
print, minmax( (h-h2)/h2 ) ;
print, mean( (h-h2)/h2 ) ;
REVISON HISTORY:
25-Aug-2006 CBO SwRI -- modified from LAY's lightcurve.pro
23-Sep-2006 LAY SwRI -- changed some variable names, added example
(See oclc_fwd_h.pro)
NAME:
oclc_fwd_hn.pro
PURPOSE:
Calculates the number density scale height, hn,
DESCRIPTION:
CALLING SEQUENCE:
hn = oclc_fwd_hn(r, t, mu, mp, dt=dt)
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).
mp - the planetary mass, in grams.
KEYWORD INPUTS AND OUTPUTS
dt - the temperature derivative at r, in Kelvin/cm
If dt is not passed, oclc_fwd_hn calculates dt from r and t.
If dt is passed but the variable is undefined, then
the calculated dt gets placed into the passed variable
(see dtarr in example 2)
OUTPUTS:
hn - number density scale height in cm
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
mp = 1.3d25 ; reference planet mass in
h = oclc_fwd_h(r, t, mu, mp)
hn = oclc_fwd_hn(r, t, mu, mp
hnarr = oclc_fwd_hn(r, replicate(t,200), mu, mp)
hndt = oclc_fwd_hn(r, replicate(t,200), mu, mp, dt=replicate(0.,200))
; ---- exact calculation, see Elliot and Young 1992
lam = (!phys.g*mp*mu*!phys.m_u/(!phys.k*t*r )) ;energy ratio
hn2 = r/lam
plot, r/1e5, (hn-hn2)/hn2
print, minmax( (hn-hn2)/hn2 ) ;
print, mean( (hn-hn2)/hn2 ) ;
plot, r/1e5, (hnarr-hn2)/hn2
print, minmax( (hnarr-hn2)/hn2 ) ;
print, mean( (hnarr-hn2)/hn2 ) ;
plot, r/1e5, (hndt-hn2)/hn2
print, minmax( (hndt-hn2)/hn2 ) ;
print, mean( (hndt-hn2)/hn2 ) ;
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
dt = b*t/r
mu = 28.01d ; scalar molecular weight
mp = 1.3d25 ; reference planet mass in g
h = oclc_fwd_h(r, t, mu, mp)
delvarx, dtarr
hnarr = oclc_fwd_hn(r, t, mu, mp, dt=dtarr)
hndt = oclc_fwd_hn(r, t, mu, mp, dt=dt)
; ---- exact calculation, see Elliot and Young 1992
lam = (!phys.g*mp*mu*!phys.m_u/(!phys.k*t *r )) ;energy ratio
hn2 = r/(lam + b)
plot, r/1e5, (hnarr-hn2)/hn2
print, minmax( (hnarr-hn2)/hn2 ) ;
print, mean( (hnarr-hn2)/hn2 ) ;
plot, r/1e5, (hndt-hn2)/hn2
print, minmax( (hndt-hn2)/hn2 ) ;
print, mean( (hndt-hn2)/hn2 ) ;
REVISON HISTORY:
25-Aug-2006 CBO SwRI -- modified from LAY's lightcurve.pro
23-Sep-2006 LAY SwRI -- changed some variable names, added example
(See oclc_fwd_hn.pro)
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
(See oclc_fwd_hydrostatic.pro)
NAME:
oclc_fwd_makelc.pro
PURPOSE:
Computes a lightcurve for a refractive atmosphere in hydrostatic
equilibrium.
DESCRIPTION:
CALLING SEQUENCE:
phi = oclc_fwd_makelc(r, t, mu, r0, p0, mp, nuSTP, distobs)
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.
nuSTP - the refractivity of the gas at STP
distobs - the distance to the observer, in cm.
OPTIONAL INPUTS:
dt - the derivative of temperature with radius
OUTPUTS:
phi - the observed light curve for a refractive atmosphere
COMMENTS:
EXAMPLE 1 - scalar temperature, molecular weight
r = dindgen(400)*1d5 + 1200d5 ; array of r, in cm
km = 1e5
distobs = 30.*1.496e8*1.d5 ; 30 AU in cm
r0 = 1250d5 ; reference r in cm
rsurf = 1000.d5
a = 0.d
b = 0.d
t0 = 104.
t = t0 * (r/r0)^b
p0 = 1.
mu = 28.01d ; scalar molecular weight
nustp = 2.98e-4 ; N2
loschmidt = 1.01325e6/(!phys.k * 273.15)
nu0 = (p0/(!phys.k * t0)) * (nustp/loschmidt)
lam0 = mu * !phys.m_u * !phys.g * mp/(!phys.k * t0 * r0)
order = 4
mp = 1.3d25 ; reference planet mass in g
phiref=oclc_fwd_makelc(r, t, mu, r0, p0, mp, nuSTP, distobs, y = y)
phiref2 = oclc_ey92_phiref_of_r(r0, nu0, lam0, a, b, distobs, rsurf, order, r)
y2 = oclc_ey92_rho_of_r(r0,nu0,lam0,a,b, distobs, order, r)
plot, r/1e5, (phiref-phiref2)
print, minmax(phiref-phiref22 )
REVISON HISTORY:
25-Aug-2006 CBO SwRI -- modified from LAY's lightcurve.pro
(See oclc_fwd_makelc.pro)
NAME:
oclc_fwd_makelcTEST.pro
PURPOSE:
Test oclc_fwd_makelc
DESCRIPTION:
CALLING SEQUENCE:
oclc_fwd_makelcTEST
INPUTS:
OUTPUTS:
REVISON HISTORY:
25-Aug-2006 CBO SwRI
(See oclc_fwd_makelc.pro)
NAME:
oclc_fwd_makelcTEST2.pro
PURPOSE:
Test oclc_fwd_makelc
DESCRIPTION:
CALLING SEQUENCE:
oclc_fwd_makelcTEST
INPUTS:
OUTPUTS:
REVISON HISTORY:
25-Aug-2006 CBO SwRI
(See oclc_fwd_makelc.pro)
NAME:
oclc_fwd_makelcTEST2.pro
PURPOSE:
Test oclc_fwd_makelc
DESCRIPTION:
CALLING SEQUENCE:
oclc_fwd_makelcTEST
INPUTS:
OUTPUTS:
REVISON HISTORY:
25-Aug-2006 CBO SwRI
(See oclc_fwd_makelc.pro)
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
(See oclc_fwd_n.pro)
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
(See oclc_fwd_nu.pro)
NAME:
oclc_fwd_phiref.pro
PURPOSE:
Calculates the refractive flux.
DESCRIPTION:
CALLING SEQUENCE:
phi = oclc_fwd_phiref(r, theta, dtheta, distobs)
INPUTS:
**** All inputs in cgs units ********
r - An array of radius values from the center of the planet, in cm.
theta - An bending angle in radians.
dtheta - The derivative of the bending angle with respect
to readius, in radians.
distobs - the distance to the observer, in cm.
OUTPUTS:
phiref - the observed flux due to refractivity only
EXAMPLE 1 - scalar temperature, molecular weight
r = dindgen(400)*1d5 + 1200d5 ; array of r, in cm
km = 1e5
distobs = 30.*1.496e8*1.d5 ; 30 AU in cm
r0 = 1250d5 ; reference r in cm
rsurf = 1000.d5
nu0 = 2d-9
lam0 = 60.d
a = 0.d
b = 0.d
order = 4
dnu = oclc_ey92_dnu(r0,nu0,lam0,a,b, r)
theta = oclc_fwd_theta(r, dnu)
dtheta = oclc_fwd_dtheta(r, theta)
phiref = oclc_fwd_phiref(r, theta, dtheta, distobs)
phiref2 = oclc_ey92_phiref_of_r(r0, nu0, lam0, a, b, distobs, rsurf, order, r)
plot, r/1e5, (phiref-phiref2)
print, minmax(phiref-phiref22 )
COMMENTS:
REVISON HISTORY:
25-Aug-2006 CBO SwRI -- modified from LAY's lightcurve.pro
(See oclc_fwd_phiref.pro)
NAME:
oclc_fwd_theta.pro
PURPOSE:
Calculates the bending angle, theta, from the radius, r, and
the derivative of the refractivity, dnudr
DESCRIPTION:
See Chamberlain and Elliot (1997) equation 2
CALLING SEQUENCE:
theta = oclc_fwd_theta(r, dnudr)
INPUTS:
**** All inputs in cgs units ********
r - An array of radius values from the center of the planet, in cm.
dnu - The derivative of the refractivity with respect to the radius.
OUTPUTS:
theta - the bending angle in radians.
COMMENTS;
calls function atmint
EXAMPLE 1 - scalar temperature, molecular weight
r = reverse(dindgen(1000)*1d5) + 1200d5 ; array of r, in cm -- OR
r = dindgen(1000)*1d5 + 1200d5 ; array of r, in cm
km = 1e5
distobs = 30.*1.496e8*1.d5 ; 30 AU in cm
r0 = 1250d5 ; reference r in cm
nu0 = 2d-9
lam0 = 60.d
a = 0.d
b = 0.d
order = 4
dnu = oclc_ey92_dnu(r0,nu0,lam0,a,b, r)
print, systime()
theta = oclc_fwd_theta(r, dnu)
print, systime()
theta2 = oclc_ey92_theta(r0,nu0,lam0,a,b,order,r) ; exact
plot, r/1e5, (theta-theta2) * distobs/km
print, minmax( theta-theta2 ) * distobs/km ;
REVISON HISTORY:
25-Aug-2006 CBO SwRI -- modified from LAY's lightcurve.pro
(See oclc_fwd_theta.pro)
NAME: oclc_fwd_y PURPOSE: Calculates the flux. DESCRIPTION: CALLING SEQUENCE: y = oclc_fwd_y(r, theta, distobs) INPUTS: **** All inputs in cgs units ******** r - An array of radius values from the center of the planet, in cm. theta - An bending angle in radians. distobs - the distance to the observer, in cm. OUTPUTS: y - shadow radius (cm) COMMENTS: EXAMPLE 1 - scalar temperature, molecular weight r = dindgen(400)*1d5 + 1200d5 ; array of r, in cm km = 1e5 distobs = 30.*1.496e8*1.d5 ; 30 AU in cm r0 = 1250d5 ; reference r in cm nu0 = 2d-9 lam0 = 60.d a = 0.d b = 0.d order = 4 dnu = oclc_ey92_dnu(r0,nu0,lam0,a,b, r) theta = oclc_fwd_theta(r, dnu) y = oclc_fwd_y(r, theta, distobs) y2 = oclc_ey92_rho_of_r(r0,nu0,lam0,a,b, distobs, order, r) plot, r/1e5, (y-y2) / km print, minmax( y-y2 ) / km ; REVISON HISTORY: 25-Aug-2006 CBO SwRI -- modified from LAY's lightcurve.pro 19-Sep-2006 LAY SwRI. Change _for_ to _fwd_, dist to distobs
(See oclc_fwd_y.pro)