Index: /trunk/doc/pslib/psLibADD.tex =================================================================== --- /trunk/doc/pslib/psLibADD.tex (revision 4206) +++ /trunk/doc/pslib/psLibADD.tex (revision 4207) @@ -1,3 +1,3 @@ -%%% $Id: psLibADD.tex,v 1.75 2005-06-09 05:16:23 eugene Exp $ +%%% $Id: psLibADD.tex,v 1.76 2005-06-10 03:25:46 eugene Exp $ \documentclass[panstarrs]{panstarrs} @@ -1087,10 +1087,10 @@ For ease of conversion, UTC should be represented as the number of seconds since the UNIX epoch of ``1970-01-01T00:00:00Z'', -non-inclusive of leap-seconds. \tbd{what does this statement -actually mean? what is the source of time (gettimeofday? let's be -explicit here} - -\tbd{Times will always be expressed in the 'UTC timezone'. Use of -the local timezone is forbidden. -- this makes no sense given that we +non-inclusive of leap-seconds. \tbd{what does this statement actually +mean? what is the source of time (gettimeofday? let's be explicit +here} + +\tbd{Times will always be expressed in the 'UTC timezone'. Use of the +local timezone is forbidden. -- this makes no sense given that we define LST. In any case, this statement, or somethign equivalent, belongs in the SDRS not the ADD} @@ -1151,5 +1151,4 @@ \subsubsection{Universal Time (UT1)} -\label{sec:ut1} UT1 is directly tied to the rotation of the Earth. Historically, time @@ -2092,4 +2091,5 @@ \subsubsection{Ray Tidal Model : {\tt psEOC\_PolarTideCorr}} +\label{Raymodel} The Ray Model tidal corrections to X, Y, and dT are given by the the @@ -2134,5 +2134,502 @@ \paragraph{Atmospheric Refraction} -\tbd{add in summary of Ken's paper} +{\em The following discussion is adapted from an article by Ken Chambers} +\newcommand\citep{\em} +\newcommand\citet{\em} + +The hypsometric structure and index of refraction of the Earth's +atmosphere produces: +\begin{itemize} +\item atmospheric refraction, resulting in an apparent positional +displacement of astronomical objects towards smaller topocentric +zenith distances, +\item chromatic dispersion, along great circles intersecting the +topocentric zenith with shorter wavelengths having smaller zenith +distances; and +\item extinction, from scattering and absorption of light by +atmospheric gases (including water vapor) and aerosols. +\end{itemize} + +Atmospheric refraction $R(\lambda) = z_{vac} - z_1$ is the difference +between the topocentric zenith angle {\it in vacuo}, where the index +of refraction is $n \equiv 1$, and the observed refracted zenith angle +at the observatory $z_1$. (All subscripts ``1" in this discussion +indicate the local values of quantities at a given observatory site, +and all units are SI.) There are various ways to express the equation +for zenith angle refraction; common ones include the ``general +equation", e.g. \citep{s1936,s1962}, the Auer \& Standish +transformation \citep{as1979,as2000}, and the Saastamoinen +approximation \citep{saas73a,saas73b,saas73c}. All are derived from +the ``refractive invariant'', the fact that along the refracted light +path described by the locus of points $r(z,n)$, where $r$ is the +radial distance from the center of the Earth and $z$ and $n$ are the +local values of the refracted zenith angle and the index of refraction +of air respectively, the product $ n r {\rm sin} z = \rm {constant} $ +remains invariant. + +We have revisted the Saastamoinen approximation in search of an +improved analytical expression accurate to $\sim 1$ mas up to zenith +angles of 75 degrees. Although different approximations are made at +different stages, both the general equation and the Saastamoinen +approximation lead to an expression for the refraction $R(\lambda)$ in +the form of a series in odd powers of tan$z_1$: +\begin{equation} + R(\lambda) = R_1 {\rm tan} z_1 + R_3 {\rm tan}^3 z_1 + + R_5 {\rm tan}^5 z_1 + R_7 {\rm tan}^7 z_1 + + R_9 {\rm tan}^9 z_1 + \dots +\end{equation} +where we have chosen to label the coefficients with a subscript reflecting +the exponent of the tan$z_1$ terms. +We find (c.f. Stone 1996) +% +\footnote{ The $R_1$ coefficient and the first two terms in the $R_2$ + coefficient in equation 2 above give results similar to that + of Stone's equation (4), if it is modified such that the + $\kappa$ factor multiplies only his $\beta$, rather than the + whole coefficient as is done for both coefficients. We + suspect this is a typographical error; his equation equation + as written gives much greater residuals when compared to the + Pulkovo Refraction Tables than he reports in his Table 2, + whereas the agreement is comparable, if his equation is + modified as described, and computed for the Pulkovo baseline + model of $\lambda = 590$nm, $15\deg C$, 101325 Pa, zero + water vapor, and mean sea level at latitude 45 degrees. An + exact comparison is not possible because the residuals in + his Table 2 are averages computed from a wide but + unspecified range of meteorological conditions. } +% +\begin{equation} + R(\lambda) = \gamma_1(1 - h_1/r_1){\rm tan} z_1 + + [\gamma_1(\gamma_1 / 2 - h_1/r_1) + \delta_1] + {\rm tan}^3 z_1 +\end{equation} +in radians, where we define for convenience the variable +$$ \gamma_1 = n_1 -1 $$ where $n_1(\lambda)$ is the index of +refraction of air at the observatory, which is dependent on the +observatory's altitude and the meteorological conditions at the time +of the observation (see Section 3 below). Each of the other variables +Eq.(2) take detailed discussions which, for the sake of clarity, are +divided into separate subsections. + +\subsection{Observatory height} + +The height of the observatory from the geometric center +of the Earth is +\begin{equation} +\begin{array}{ll} + r_1 = r_e + r_h, \qquad r_h \approx r_n + r_o, +\end{array} +\end{equation} +where +$r_e$ is the local radius of the reference ellipsoid, +$r_h$ is the local height above the reference ellipsoid, +$r_n$ is the local geoid height, normal to the reference ellipsoid, +$r_o$ is the orthometric height, which is the height above the geoid + (or Mean Sea Level) in the direction of normal gravity + (i.e. a local plumb line). + +\subsection{The magnitude of normal gravity at the observatory } + +The local magnitude of normal gravity +% +\footnote{ Stone uses a common + \citep[e.g.][]{allen1973,seid1992,allen2001} expression for the + normal gravity as a function of lattitude and altitude, which + apparently first appeared in \citet{lamb1949}. However its + derivation is unclear -- cited as {\it ``Some notes on the + calculation of the geopotential", unpublished manuscript}, + \citet{lamb1949}. } +% +$g_1$ +is the acceleration due to the combination +of the gravity of the ellipsoid, the local mass distribution, +and the centripetal acceleration from the Earth's rotation, +the vector sum being directed opposite to the local zenith by definition, +being close but not identical to the normal to the ellipsoidal, but +not intersecting the the geometric center of the Earth. +i.e. the atmospheric topocentric zenith differs from the geocentric +astronomical zenith, the impact of this difference on calculating +the atmospheric refraction is minimal, see \citep{seid1992}. +The acceleration at the observatory is +\begin{equation} + g_1(r_h,\phi) = g(\phi)\left[ 1 - 2( 1 + f + m_r - 2 f {\rm sin}^2\phi) + \left(r_h \over a \right) + + 3 \left(r_h \over a \right)^2 \right] +\end{equation} +where +\begin{equation} + g(\phi) = g_e \left( { {1 + k_s {\rm sin}^2 \phi}\over + \sqrt{1 - \epsilon^2 {\rm sin}^2 \phi}} \right) +\end{equation} +and $\phi$ is the latitude of the observatory and the other constants +are given in Table 2. +\begin{table}[!ht] +\caption{WGS-84 World Geodetic System Reference Ellipsoid for GPS} +\smallskip +\begin{center} +{\small +\begin{tabular}{lr} +%\tableline +\noalign{\smallskip} +Semi-major axis $a$ & 6356752.3142 $m$\\ +Flattening $f =(a-b)/a$ & 0.003352811 \\ +Eccentricity $\epsilon = \sqrt{a^2 +b^2}/a$ & 0.081819\\ +Polar gravity $g_p$ & 9.8321849378 $m s^{-2}$ \\ +Equatorial gravity $g_e$ & 9.7803253359 $m s^{-2}$ \\ +Somigliana's Constant $k_s = ((b/a)(g_p/g_e)-1)$ & 0.001931853 \\ +Angular velocity of the Earth $\omega$ & 0.00007292115 $rad/s$\\ +$GM_{earth}$ including atmosphere & 3986004.418$\times 10^8 m^3 +s^{-2}$\\ +Gravity ratio $m_r = \omega^2 a^2 b /(GM)$ & 0.003449787\\ +\noalign{\smallskip} +%\tableline +\end{tabular} +} +\end{center} +\end{table} + + +\subsection{The scale height above the observatory} + +The scale height of the atmosphere above the observatory is +\begin{equation} + h_1 = \mathcal{Z}_1\mathcal{R}T_1/g_1 M_a[1 - x_w(1-M_w/M_a)], +\end{equation} +where +$T_1$ is the local air temperature at the observatory in degrees Kelvin, +$\mathcal{R} = 8.314472 {\rm J \ mol^{-1} K^{-1}}$ +is the gas constant, +\begin{equation} +\begin{array}{ll} \mathcal{Z}_1 = 1 - (P_1/T_1)\left[a_0 + a_1t_1 + a_2t_1^2 + + (b_0 + b_1)t_1 x_w + (c_0+c_1t_1)x_w^2 \right] \\ + \qquad \qquad +(P_1/T_1)^2(d+ex_w^2) + \end{array} +\end{equation} +is the compressibility of moist air at local air temperature +$t_1 = T_1 - 273.15$ +in degrees Celsius. +The values of the other constants in $\mathcal{Z}_1$ are given in +Appendix A. +The quantity +$M_a = 0.0289635 + 1.2001 \times 10^{-8}(x_{c1} - 400)$ +kg/mole +is the molar mass of dry air with a $CO_2$ concentration $x_{c1}$ +in $\mu$mol/mol, +$M_w = 0.018015 $ kg/mole is the molar mass of water vapor, +and +$x_w$ is the molar fraction of water vapor in moist air, which +depends on the local humidity. + +To calculate both the scale height $h_1$ and +the index of refraction of moist air at the observatory $n_1$ +we need to calculate $x_{w1}$, the molar fraction of water vapor +from local measurements of the relative humidity or, much better, the +dew point. Simple formulas such as \citet{davis1992} +are only suitable above $0\deg C$, whereas at Mauna Kea observatories the +temperature is often below $0\deg C$, and occasionally the surrounding +ground is covered in snow or ice, which also alters saturation vapor pressure. +Thus, we adopt the best available equation for the +saturation water vapor pressure $p_{sv}$, +the IAPWS equations \citep{huang1998}. +Haung's equations are: + +\begin{equation} +\begin{array}{ll} + \Omega = T + K_9/(T - K_{10}) \qquad A = \Omega^2 +K_1\Omega + K_2 \\ + B = K_3\Omega^2 + K_4\Omega + K_5 \qquad C = K_6\Omega^2 + K_7\Omega + K_8 \\ + X = - B + \sqrt{B^2 - 4AC} \qquad p_{sv}(t) = 10^6(2C/X)^4 \\ +\end{array} +\end{equation} +For saturation vapor pressure over ice or snow, use +\begin{equation} +\begin{array}{ll} + \Theta = T/273.16 \qquad Y=A_1(1-\Theta^{-1.5})+A_2(1-\Theta^{-1.25}) \\ + p_{sv}(t) = 611.657 e^Y \qquad \\ +\end{array} +\end{equation} +The constants are given in Table 1. +\begin{table}[!ht] +\caption{Constants for Compressibility and Humidity Equations} +\smallskip +\begin{center} +{\small +\begin{tabular}{lcl} +%\tableline +\noalign{\smallskip} + $a_0=1.58123\times 10^{-6}$KPa$^{-1}$ & & $K_1= 1.16705214528E+03$ \\ + $a_1=-2.9331\times 10^{-8}$Pa$^{-1}$ & & $K_2= -7.24213167032E+05$ \\ + $a_2=1.1043\times 10^{-8}$K$^-1$Pa$^{-1}$ & & $K_3= -1.70738469401e+01$ \\ + $b_0=5.707\times 10^{-6}$KPa$^{-1}$ & & $K_4= 1.20208247025E+04$ \\ + $b_1=-2.051\times 10^{-8}$Pa$^{-1}$ & & $K_5= -3.23255503223E+06$ \\ + $c_0=1.9898\times 10^{-4}$KPa$^{-1}$ & & $K_6= 1.49151086135E+01$ \\ + $c_1=-2.376\times 10^{-6}$Pa$^{-1}$ & & $K_7= -4.82326573616E+03$ \\ + $d =1.83\times 10^{-11}$K$^2$Pa$^{-2}$ & & $K_8= 4.05113405421E+05 $ \\ + $e =-0.765\times 10^{-8}$$^2$KPa$^{-2}$ & & $K_9= -2.38555575678E-01$ \\ + & & $K_{10}=6.50175348448E+02 $ \\ + & & $A_1=-13.928169$ \\ + & & $A_2=34.7078238$ \\ +\noalign{\smallskip} +%\tableline +\end{tabular} +} +\end{center} +\end{table} + +\noindent Now, to calculate the mole fraction of water vapor $x_w$, we need the +so called ``enhancement factor" +\begin{equation} + f(p,t)= a' + b' p + c' t^2 +\end{equation} +where $a' = 1.00062$,$b'= 3.14 \times 10^{-8}$, and +$c' = 5.60\times 10^{-7}$, and where $p$ and $t$ are the air pressure in +Pascals and air temperature. + +If you have the more precise measurement of the dew point $t_d$ (or frost point)then +\begin{equation} +x_w = f(p,t) \times p_{sv}(t_d)/p. +\end{equation} +On the other hand if you only have the relative humidity $RH$, a less accurate +expression is +\begin{equation} +x_w = (RH/100) \times f(p,t) \times p_{sv}(t_d)/p. +\end{equation} + +\subsection{The index of refraction of moist air at the observatory} + +The Ciddor equation for the index of refraction of moist air +\citep{ciddor1996} has been adopted by the International Association +of Geodesy (IAG) as the standard as it is believed to provide the most +accurate results under the largest +range of wavelength, temperature, and humidity conditions (300 to 1690 nm; +-40 to 100 C; 0-80\% RH). Note for astronomy, the air temperature at the +tropopause in the 1976 standard atmosphere is 216.6 K, below the stated range ofvalidity. Astronomical observers often observe up to 90\% RH, where +water droplets can form and change the effective index of refraction. +Note that developments in the equation for the index of refraction of moist +air have been poorly tracked in the astronomical literature and it is critical to examine in +every application what equation is actually being used. \footnote{ +Edlen's (1953) original fit to the available data covered the wavelength +range 2752 to 6440 /AA with reasonable residuals. His 1953 constants still +survive in the astronomical literature in the equations of \citet{allen1973}; +%Allen (1973); +\citet{stone1996}; \citet{allen2001}; +%Stone(1996); Allen(2001), +Roe (2002 - who extrapolates the 1953 equation to +K band to correct for dichoric adaptive optics) and in the computer codes +SLALIB and ZEEMAX. None discuss the range of validity. However, Elden himself +revised them \citep{elden1966}, +%(Elden 1966), +and these were further discussed and updated by +\citet{pr1972, bd1993, bd1994, ciddor1996, bp1998}. +\citet{rueg1998} and \citet{sz2004} make convincing arguments +for the Ciddor equation, and the latter's approach is followed here.} + +The index of refraction of air at standard temperature and pressure +is \citep{ciddor1996} +\begin{equation} +\gamma_{as} = 10^{-8} \left( { \left[ { k_1 \over k_0 - \sigma^2} \right] + + \left[ { k_3\over k_2 - \sigma^2} \right] } \right), +\end{equation} +where $\sigma = 1/\lambda$ is the wavenumber of wavelength of light in microns. +Adjusting for the (annually varying and secularly increasing) value of +atmospheric $CO_2$ concentration $x_{CO2}$ in units of $\mu$mole/mole, +the expression for dry air becomes +\begin{equation} +\gamma_{axs} = \gamma_{as} + \left[ 1 + 5.34\times10^{-7} + (x_{CO2} - 450\ \mu{\rm mole/mole} )\right]. +\end{equation} +For water vapor under standard conditions, Ciddor finds +\begin{equation} +\gamma_{ws} = 1.022 \times 10^{-8} + \left[\omega_0 + \omega_1 \sigma^2 + + \omega_2 \sigma^4 + \omega_3 \sigma^6 \right]. +\end{equation} +Following \citep{owens1967}, +the indicies can be combined in proportion to their densities, +thus the index of refraction of moist air at the observatory +$n_1 = \gamma_1 +1 $ is given by +\begin{equation} +\gamma_1 = (\rho_a/ \rho_{axs})\gamma_{axs} + + (\rho_w/ \rho_{ws})\gamma_{ws}, +\end{equation} +where +\begin{equation} +\begin{array}{ll} + \rho_{a} = (1 - x_w) P_1 M_a/({\cal Z}_m {\cal R} T_1) \\ + \rho_{w} = x_w P_1 M_w/( {\cal Z}_m {\cal R} T_1) \\ + \rho_{axs} = P_{STP} M_a/( {\cal Z}_a {\cal R} T_{STP}) \\ +\end{array} +\end{equation} +and $\rho_{ws}$ is given in Table 3. + +\begin{table}[!ht] +\caption{Constants in the Ciddor Eq. for index of refraction of moist air} +\smallskip +\begin{center} +{\small +\begin{tabular}{lll} +%\tableline +\noalign{\smallskip} +$k_0 = 23.0185\ \mu {\rm m}^{-2}$ & +$\omega_0 = 295.235\ \mu {\rm m}^{-2}$ & +$P_{STP} = 101325$ Pa \\ + +$k_1 = 5792105\ \mu{\rm m}^{-2} $ & +$\omega_1 = 2.6422\ \mu {\rm m}^{-2}$ & +$T_{STP} = 288.15 $ K \\ + +$k_2 = 57.362\ \mu{\rm m}^{-2} $ & +$\omega_2 = -0.03238\ \mu {\rm m}^{-4}$ & +${\cal Z}_{a} = 0.9995922115$ \\ + +$k_3 = 167917\ \mu {\rm m}^{-2} $ & +$\omega_3 = 0.004028\ \mu {\rm m}^{-6}$ & +$\rho_{ws} = 0.00985938\ {\rm kg\ m}^3$ \\ +\noalign{\smallskip} +%\tableline +\end{tabular} +} +\end{center} +\end{table} + + +\subsection{The tropopause term in the equation of refraction} + +The final term in the Refraction Equation (2) +is $\delta_1$, which comes from our re-derivation of the +Saastamoinen approximation: +\begin{equation} + \delta_1 = 5 \left( {h_1 \over {r_1 T_1}}\right)^2 + \left[ {{{\gamma_{ft}T_{ft}^2 - \gamma_t T_t^2} \over + {1 - (h_1\beta/T_1)}}} + + \gamma_tT_t^2 \right] +\end{equation} +where $\beta$ is the lapse rate of the troposphere, +and the index of refraction of the free troposphere is given by +\begin{equation} + \gamma_{ft} = \gamma_t(T_{ft}/T_t) ^{-(T_1/h_1 \beta)-1} +\end{equation} +and the index of refraction of the tropopause is given by +\begin{equation} + \gamma_t = \gamma_1 + {\rm exp} \left[ {T_1 (r_t - r_1) \over T_t h_1 }\right] +\end{equation} +where $\beta$ is the lapse rate in K/m, +the temperature of the tropopause $T_t$, and height of the tropopause $r_t$ +are all determined from contemporaneous meterological data +(radiosonde or modern forecast models). Then the temperature of the +free troposphere is given by +\begin{equation} + T_{ft} = T_t - \beta(r_t - r_1). +\end{equation} + +\subsection{Calculating the atmospheric refraction from both + the observed and true zenith angle} + +The monochromatic refraction can now be calculated for any given wavelength +$\lambda_{air}$ +(formally only within the range of validity - 300 to 1670 nm) given +the altitude of the observatory +$h_1$; +contemporaneous meteorological measurements at the observatory of +air temperaure $T_1$ (K); +atmospheric pressure $P_1$ (Pa); +percent relative humidity $RH$, or preferably dew point temperature $t_d$ ($\deg$C); +as well as a small set of additional meteorological data: +the lapse rate of the troposphere $\beta$ (K/m); +the radius of the tropopause $r_t = r_e+h_t$(m), where +$h_t$ is the height above MSL of the tropopause; +the temperature of the tropopause $T_t$(K); +and the atmospheric concentration of CO$_2$ +$x_{CO2}$ ($\mu$mole/mole). +The characterization of the troposphere and tropopause can +be determined either by interpolation between radiosonde measurements +or from hourly updated meteorological models available at many observatory +sites, or, worst case, simply adopting the 1976 Standard Atmosphere +values:\footnote{ +\citet{seid1992} +%Seidelmann 1992 +contains a typographical error quoting this value in units of K/km} +$\beta = -0.0065$ K/m, +$h_t = 11000$m, +$T_t = 216.6$ K and assuming +$x_{C02} = 375 \mu$mole/mole.\footnote{ +Closed rooms have higher CO$_2$ concentration, thus the STP laboratory measurements +have concentrations near 450 $\mu$mole/mole (the prefered unit to parts per million per volume). The secular +increase of atmospheric CO$_2$ in the industrial age is well documented, with +an annual cycle superimposed due to terrestrial biomass and ocean exchange. +} + + +If the number of electrons being created from the illumination from the source in +the interval of wavelength $d\lambda$ is $N_{\lambda}d\lambda$, +then the mean refraction is +\begin{equation} +\bar R = {{\int R(\lambda) N_{\lambda} d \lambda} \over + {\int N_{\lambda} d\lambda}} +\end{equation} + +\section{Atmospheric Dispersion} + +The atmospheric dispersion is then +\begin{equation} +{ \overline{(R- \bar R)^2}} = {{\int (R - \bar R)^2 N_{\lambda} d \lambda} \over + {\int N_{\lambda} d\lambda}} +\end{equation} + +\section{Air Mass and Extinction} + +By Laplace's theorem, the monochromatic airmass (mass per unit area +along the refracted path) is +\begin{equation} +M(z_1) = (P_1/g_1) R(\lambda) / {\rm sin} z_1 +\end{equation} +in kg/m$^2$. Thus +\begin{equation} + M(z_1) = (P_1/g_1) \left( \gamma_1(1 - h_1/r_1){\rm sec} z_1 + + [\gamma_1(\gamma_1 / 2 - h_1/r_1) + \delta_1] + {\rm sec}^3 z_1 \right). +\end{equation} +This is generally normalized to $P_1/g_1$, i.e. in spite of the +daily changes in barometric pressure, and thus daily changes in the +true mass of air over the observatory, the resulting change in extinction is +generally treated as a drift in photometric zeropoint. +For a survey program like Pan-STARRS, one could instead normalize to a +standard barometric pressure (i.e. the altitute pressure), and thus +during a low in atmospheric pressure, the airmass would be less than +1 at zenith, and during periods of high pressure the airmass would +be greater than 1 at zenith. From the variation with zeropoint and +temperature and barometric pressure, this would remove most of the +observed variation in zeropoint in the CFHT legacy program. +(Magnier, private comunication). + +The mean airmass is then +\begin{equation} +\bar M(z_1) = \left( P_1 \over g_1 \right) \int + \left( \gamma_1(1 - h_1/r_1){\rm sec} z_1 + + [\gamma_1(\gamma_1 / 2 - h_1/r_1) + \delta_1] + {\rm sec}^3 z_1 \right) N_{\lambda} d \lambda. +\end{equation} +and depends weakly on the filter bandpass. Use of this more +accurate expression for airmass should lead to improved extinction +corrections at high airmass. + +\section{Limits to ground based relative and absolute astrometry} + +The limits to ground based astrometry may well be our abilitiy to +measure the atmospheric profile along the line of sight of a given +observation, and the systematic limit of the telescope axes encoders +(and sophistication of the telescope mount model.) +The refratction model above requires only the additional data of +the temperature, height, and presure of the tropopause, but much more +detailed atmospheric information will be available for PS1 from our sky +probes which measure atmospheric absorption for each field and even, +phase drifts of GPS clocks from Rubidium or Cesium standard clocks. +These can be converted directly into a nearby line of sight index of refraction +at optical wavelengths. Thus we encourage wide field survey telescopes to +err on the side of over instrumenting the accuracy and repeatibility +of the axes encoders. + + \subsection{Projections}