In [1]:
import numpy as np
import matplotlib.pyplot as plt
from SciServer import CasJobs
from SciServer import Authentication
from io import StringIO
from astropy.table import Table
import requests
Based on answer by Dou Liu
The equations in Stone (1996) relate index of refraction to zenith angle, and give the dependence of index of refraction on local conditions and wavelength. They claim this to be accurate under most conditions better than 0.150 arcsec for zenith angles less than 75 deg, and better than 0.010 arcsec for zenith angles less the 65 deg.
These equations allow you to relate the apparent zenith angle (which is what the star would be at if there were no atmosphere) to the observed zenith angle (which is what the star should be at accounting for the atmosphere). The observed zenith angle is always smaller than the apparent zenith angle. We define the refraction as the apparent minus the observed.
Note that astropy, starlink, IRAF, and other astronomical tools have built into them refraction calculators. Those other tools are the ones you should actually use!
$\kappa$ is the ratio of the gravity at the observing site to the sea-level gravity at the earth's equator. $$ \kappa=\frac{g_0}{g}=1+0.005302\sin^2\phi-5.83\times 10^{-6}\sin^2(2\phi)-3.15\times10^{-7}h $$
In [2]:
def kappa(phi, h):
"""Gravitational constant at observatory relative to equator
Parameters
----------
phi : np.float32
latitude of observatory (deg)
h : np.float32
height of observatory from sea level (meters)
"""
k = 1 + 5.302 / 1E3 * np.sin(phi)**2 - 5.83 / 1E6 * np.sin(2. * phi)**2 - 3.15 / 1E7 * h
return k
The atmosphere height effect is account for with $$ \beta=\frac{H}{R_\text{earth}}=\frac{1}{R_\text{earth}}\int_0^\infty e^{-h/h_0}dh=\frac{0.01254}{273.15}T $$
In [3]:
def beta(T):
"""Scaling of temperature
Parameters
----------
T : np.float32
temperature in K
"""
return(4.5908 / 1E6 * T)
The angular difference in arcsec due to chromatic atmospheric difference refraction is $$R(\lambda)\simeq \kappa(n(\lambda)-1)(1-\beta)\tan z-\kappa(1-n(\lambda)\left(\beta-\frac{1}{2}(n(\lambda)-1)\right)\tan^3z$$ where $z$ is the zenith angle and $n$ is the index of refraction.
This formula accounts for the spherical shape of the Earth and atmosphere, which is why the gravity and the height enter the calculation.
In [4]:
def refraction(lam=5500., T=285., Ps=1000., RH=0.2,
phi=0., h=0., z=0., n=None):
"""Refraction calculation
Parameters
----------
lam : np.float32
wavelength (A) (default 5500.)
T : np.float32
temperature in K (default 285.)
Ps : np.float32
pressure (mbar) (default 1000.)
RH : np.float32
relative humidity (default 0.2)
phi : np.float32
latitude (deg) (default 0.)
h : np.float32
height from sea level (default 0.)
z : np.float32
zenith angle (deg) (default 0.)
n : np.float32
index of refraction (by default sets with local_ref())
Returns:
-------
R : np.float32
refraction in arcsec (difference of apparent to observed zenith angle)
"""
if(n is None):
n = local_ref(lam, T, Ps, RH)
zrad = z * np.pi / 180.
k = kappa(phi, h)
b = beta(T)
R = k * (n - 1) * (1 - b) * np.tan(zrad) - k * (1 - n) * (b - (n - 1) / 2) * np.tan(zrad)**3
R = R * 180. / np.pi * 3600.
return R
The index of refraction of air $n$ is computed accurately with the relations $$ n=1+10^{-8}\left[(2371.34+\frac{683939.7}{130-\sigma^2}+\frac{4547.3}{38.9-\sigma^2}\right]D_s +(6487.31+58.058\sigma^2-0.71150\sigma^4+0.08851\sigma^6)D_w) P_w=RH\times10^{4}\times \exp{77.3450+0.0057T-7235.0/T}/T^{8.2}\\ D_s=\left[1+(P_s-P_w)(57.90\times10^{-8}-\frac{9.3250\times10^{-4}}{T}+\frac{0.25844}{T^2})\right]\frac{P_s-P_w}{T}\\ D_w=\left[1+P_w(1+3.7\times10^{-4}P_w)\left(-2.37321\times10^{-3}+\frac{2.23366}{T}-\frac{710.792}{T^2}+\frac{7.75141\times10^4}{T^3}\right)\right]\frac{P_w}{T} $$ where $\sigma$ is the wave number $$ \sigma=\frac{10^4}{\lambda} $$ where $D_s$ and $D_w$ are the density factors for dry air and water vapor.
In [5]:
def local_ref(lam=5500., T=285., Ps=1000., RH=0.2):
"""Local refraction index
Parameters
----------
lam : np.float32
wavelength (A) (default 5500.)
T : np.float32
temperature in K (default 285.)
Ps : np.float32
pressure (mbar) (default 1000.)
RH : np.float32
relative humidity (default 0.2)
Returns:
-------
n : np.float32
index of refraction (by default sets with local_ref())
"""
sigma=10**4/lam
Pw=RH/1E4*np.exp(77.3450+0.0057*T-7235.0/T)/T**8.2
Ds=(1+(Ps-Pw)*(57.90/1E8-9.3250*1E-4/T+0.25844/T**2))*(Ps-Pw)/T
Dw=(1+Pw*(1+3.7/1E4*Pw)*(-2.37321/1E3+2.23366/T-710.792/T**2+7.75141*10**4/T**3))*Pw/T
n=1+((2371.34+683939.7/(130-sigma)+4547.3/(38.9-sigma**2))*Ds
+(6487.31+58.058*sigma**2-0.71150*sigma**4+0.08851*sigma**6)*Dw)*1E-8
return n
It is instructive to compare this refraction to what we would derive just from Snell's law for a plane-parallel atmosphere. This is much simpler to derive. We start with: $$ n_{\rm space} \sin z_{\rm app} = n_{\rm obs} \sin z_{\rm obs} $$ and noting that $n_{\rm space} = 1$: $$ z_{\rm obs} = \sin^{-1} \left(\frac{\sin z_{\rm app}}{n_{\rm obs}}\right) $$ Note that in the plane-parallel case, only the index of refraction local to the observatory matters.
It is easy to show that for $n \approx 1$: $$ R = z_{\rm app} - z_{\rm obs} \approx (n - 1) \tan z_{\rm app} $$ This shows why the leading term of Stone's equation looks the way that it does.
In [6]:
def refraction_pp(lam=5500., T=285., Ps=1000., RH=0.2,
z=0., n=None):
"""Refraction calculation in plane-parallel approximation
Parameters
----------
lam : np.float32
wavelength (A) (default 5500.)
T : np.float32
temperature in K (default 285.)
Ps : np.float32
pressure (mbar) (default 1000.)
RH : np.float32
relative humidity (default 0.2)
z : np.float32
zenith angle (deg) (default 0.)
n : np.float32
index of refraction (by default sets with local_ref())
Returns:
-------
R : np.float32
refraction in arcsec (difference of apparent to observed zenith angle)
"""
if(n is None):
n = local_ref(lam, T, Ps, RH)
zrad = z * np.pi / 180.
zobsrad = np.arcsin(np.sin(zrad) / n)
R = (zrad - zobsrad) * 180. / np.pi * 3600.
return R
Just to check, we can look at the index of refraction as a function of wavelength. $n-1$ is very close to zero, of order $10^{-4}$, and importantly it changes significantly with wavelength, getting smaller at higher wavelength. This is the source of the chromatic effects. Note that the chromatic effects in this range are about $10^{-2}$ of the total index of refraction.
In [7]:
nlam = 1000
lamstart = 3500.
lamend = 20000.
lam = lamstart + (lamend - lamstart) * np.arange(nlam) / np.float32(nlam - 1)
n = local_ref(lam=lam)
plt.plot(lam, n - 1.)
plt.xlabel('Wavelength (Ang)')
plt.ylabel('n - 1')
Out[7]:
Now let us look at the refraction as a function of zenith angle. We can see a couple of things. First, it is a pretty large effect! The observed images are moved many arcsecs from their apparent position. Second, the plane parallel approximation is pretty good, within a fraction of an arcsec until $z>40$ deg.
In [8]:
nz = 100
z = 75. * np.arange(nz) / np.float32(nz)
Rpp = refraction_pp(z=z)
R = refraction(z=z)
fig, ax = plt.subplots(1, 2, figsize=(9, 4))
ax[0].plot(z, R, color='red', label='Spherical')
ax[0].plot(z, Rpp, color='black', label='Plane-parallel')
ax[0].set_xlabel("Zenith angle (deg)")
ax[0].set_ylabel("Refraction (arcsec)")
ax[0].legend()
ax[1].plot(z, R - Rpp)
ax[1].set_xlabel("Zenith angle (deg)")
ax[1].set_ylabel("Spherical minus Plane-parallal (arcsec)")
Out[8]:
Now let us look at airmass 1.2, which is zenith angle of $z \approx 33$ deg, and consider the refraction as a function of wavelength.
The differential refraction across the optical and near-infrared bands is about 0.5 arcsec, not coincidentally about $10^{-2}$ of the total refraction effect (because this is the fractional variation of $n$ in that range). Between 4000 and 8000 Angstroms, the image shifts about 0.25 arcsec, in the sense that the bluer light is higher in the sky than the redder light.
Aside from a very slight overall difference of $\sim 0.02$ arcsec in the absolute refraction, the plane-parallel approximation captures the variation with wavelength very accurately.
In [9]:
z = 33.
Rpp = refraction_pp(z=z, lam=lam)
R = refraction(z=z, lam=lam)
fig, ax = plt.subplots(1, 2, figsize=(16, 4))
ax[0].plot(lam, R, color='red', label='Spherical')
ax[0].plot(lam, Rpp, color='black', label='Plane-parallel')
ax[0].set_xlabel("Wavelength (Ang)")
ax[0].set_ylabel("Refraction (arcsec)")
ax[0].legend()
ax[1].plot(lam, R - Rpp)
ax[0].set_xlabel("Wavelength (Ang)")
ax[1].set_ylabel("Spherical minus Plane-parallal (arcsec)")
Out[9]:
Based on answer by David Mykytyn
We first define a function to retrieve the data we are interested in for a given run.
In [10]:
def retrieve_camcol(run=756, camcol=3):
# We define the columns we want and their
columns = ('ra', 'dec', 'airmass_g', 'airmass_i', 'sky_g', 'sky_i', 'field')
dtypes= ('f8', 'f8', 'f4', 'f4', 'f4', 'f4', 'i4')
# Now define the query
query = """
SELECT {columns}
FROM Field
WHERE run = {run} and camcol = {camcol}
"""
query = query.format(columns=', '.join(list(columns)),
run=run, camcol=camcol)
# Execute the query (requires internet access)
responseStream = CasJobs.executeQuery(query, "DR12", format="dict")
# convert result into astropy table
result = responseStream['Result'][0]
data = list(map(list, zip(*result['Data'])))
fields = Table(data, names=columns, dtype=dtypes)
isort = np.argsort(fields['field'])
fields = fields[isort]
return(fields)
Now we define a function to plot the asked-for results.
In [11]:
def plot_sky(fields=None):
field = fields['field']
sky_g = 22.5 - 2.5 * np.log10(fields['sky_g'])
sky_i = 22.5 - 2.5 * np.log10(fields['sky_i'])
airmass_g = fields['airmass_g']
airmass_i = fields['airmass_i']
fig, axes = plt.subplots(2, 2, figsize=(10, 10))
axes[0, 0].plot(field, sky_g)
axes[0, 0].set_xlabel('Field Number')
axes[0, 0].set_ylabel('g_band sky (mag/arcsec$^2$)')
axes[0, 1].plot(airmass_g, sky_g, '.')
axes[0, 1].set_xlabel('airmass g')
axes[0, 1].set_ylabel('g_band sky (mag/arcsec$^2$)')
axes[1, 0].plot(field, sky_i)
axes[1, 0].set_xlabel('Field Number')
axes[1, 0].set_ylabel('i_band sky (mag/arcsec$^2$)')
axes[1, 1].plot(airmass_i, sky_i, '.')
axes[1, 1].set_xlabel('airmass i')
axes[1, 1].set_ylabel('i_band sky (mag/arcsec$^2$)')
First, let us consider a run that turns out to be taken at almost fixed airmass the entier time. We see substantial variation in sky brightness in both bands that correlates with field (time) in the same manner for both. This indicates sky that is coming from scattered light that is varying with time.
In [12]:
fields = retrieve_camcol(run=4905)
plot_sky(fields)
But we can also consider a run where there is substantial variation of airmass. There is often a clear correlation of brightness with airmass. As you get further from zenith, you expect the brightness to increase. Note that this does not happen with the exact factor you expect (linear with airmass), because the sky is not uniform and it also varies, even for the $g$ band. For the $i$ band, the OH emission is even more variable in space and time, and so you see extra variation there.
In [13]:
fields = retrieve_camcol(run=4822)
plot_sky(fields)
In [ ]: