In [1]:
from __future__ import division,unicode_literals # Half arsed attempt at Python 2 compatibility
In [2]:
%matplotlib inline
In [3]:
import numpy as np
from scipy.interpolate import RectSphereBivariateSpline, SmoothBivariateSpline, InterpolatedUnivariateSpline
from matplotlib import pyplot as plt
from matplotlib import rcParams
from matplotlib import colors
from astropy.io import fits
from astropy.table import Table
from astropy import units as u
from astropy.constants import h,c
import healpy as hp
In [4]:
rcParams['figure.figsize'] = 12, 9
End-to-end throughput of the instrument from the entrance aperture through to photo-electrons in the image sensor
In [5]:
pupil_area = (np.pi * (90 * u.mm)**2 / 4).si
pupil_area
Out[5]:
In [6]:
n_surfaces = 9 # Doesn't include last surface
throughput = 0.995**n_surfaces * u.dimensionless_unscaled
throughput
Out[6]:
In [7]:
# Load DECam filter data
filters = [Table.read('http://svo2.cab.inta-csic.es/theory/fps/fps.php?ID=CTIO/DECam.u'), \
Table.read('http://svo2.cab.inta-csic.es/theory/fps/fps.php?ID=CTIO/DECam.g'), \
Table.read('http://svo2.cab.inta-csic.es/theory/fps/fps.php?ID=CTIO/DECam.r'), \
Table.read('http://svo2.cab.inta-csic.es/theory/fps/fps.php?ID=CTIO/DECam.i'), \
Table.read('http://svo2.cab.inta-csic.es/theory/fps/fps.php?ID=CTIO/DECam.z')]
for f in filters:
# Set units correctly
f['Wavelength'] = f['Wavelength'].to(u.um)
f['Wavelength'].unit = u.um
f['Transmission'].unit = u.dimensionless_unscaled
In [8]:
for f in filters:
plt.plot(f['Wavelength'], f['Transmission'])
plt.xlim(0.25,1.05)
plt.xlabel('Wavelength/$\mu$m')
plt.ylabel('Transmission')
plt.title('DECam filter transmission')
Baseline concept uses an e2V CIS115, a 2000 x 1504 pixel back-illuminated monolithic CMOS image sensor with 7 micron pixel pitch. Numerical QE data is hard to come by but there are some plots available from which we can extract approximate date with http://arohatgi.info/WebPlotDigitizer/, the most useful of which appears to be in Soman et al (2014), http://proceedings.spiedigitallibrary.org/proceeding.aspx?articleid=1891353, which has full wavelength range predicted QE at both +20$^\circ$C and -40$^\circ$C. These data seem quite conservative as Wang et al (2014), http://proceedings.spiedigitallibrary.org/proceeding.aspx?articleid=1891414, get higher peak QE from both e2v's and their own measurements on a prototype device.
In [11]:
qe_p20C = Table.read('resources/cis115_QE_20C.csv')
qe_m40C = Table.read('resources/cis115_QE_-40C.csv')
qe_p20C['Wavelength'].unit = u.um
qe_m40C['Wavelength'].unit = u.um
qe_p20C['QE'].unit = u.electron / u.photon
qe_m40C['QE'].unit = u.electron / u.photon
In [12]:
plt.plot(qe_p20C['Wavelength'], qe_p20C['QE'], 'r-', label='$+20^\circ$C')
plt.plot(qe_m40C['Wavelength'], qe_m40C['QE'], 'b-', label='$-40^\circ$C')
plt.xlim(0.25,1.05)
plt.ylim(0,1)
plt.xlabel('Wavelength/$\mu$m')
plt.ylabel('QE / e photon$^{-1}$')
plt.legend(loc='best')
plt.title('e2v CIS115 predicted QE (Soman et al)')
plt.savefig('QE.png')
In [13]:
efficiencies = []
for f in filters:
# Interpolate QE data at same wavelengths as filter transmission data
i = np.interp(f['Wavelength'], qe_m40C['Wavelength'], qe_m40C['QE']) * qe_m40C['QE'].unit
efficiencies.append(throughput * f['Transmission'].quantity * i)
In [14]:
for i, f in enumerate(filters):
plt.plot(f['Wavelength'], efficiencies[i])
plt.xlim(0.25,1.05)
plt.ylim(0,1)
plt.xlabel('Wavelength/$\mu$m')
plt.ylabel('Efficiency / e photon$^{-1}$')
plt.title('End-to-end efficiency')
plt.savefig('eff.png')
In [15]:
eff_areas = [pupil_area * efficiency for efficiency in efficiencies]
In [16]:
eff_etendues = [eff_area * (3 * u.arcsecond)**2 / u.pixel for eff_area in eff_areas]
In space the visible wavelength sky background is dominated by zodiacal light. To model the zodiacal light we follow the same prescription as the HST ETC, as described in Giavilisco, Saul & Bohlin (2002), http://www.stsci.edu/hst/wfc3/documents/ISRs/2002/WFC3-2002-12.pdf.
In the visible wavelength range the zodiacal light is the result of scattered sunlight. Consequently the starting point is a solar spectrum from Colina, Bohlin & Castelli (1996), http://adslabs.org/adsabs/abs/1996AJ....112..307C/
In [17]:
sun = fits.open('resources/sun_castelli.fits')
In [18]:
sun.info()
In [19]:
sun[0].header
Out[19]:
In [20]:
sun[1].header
Out[20]:
In [21]:
sun_waves = sun[1].data['WAVELENGTH'] * u.Angstrom
In [22]:
sun_flux = sun[1].data['FLUX'] * u.erg * u.second**-1 * u.cm**-2 * u.Angstrom**-1
To convert the solar spectrum to a zodical light spectrum it is normalised and reddened following the prescription of Leinert et al (1997), http://www.edpsciences.org/10.1051/aas:1998105, with the revised parameters from Aldering (2001), http://www-supernova.lbl.gov/~aldering/zodi.pdf. The result is normalised to the correct surface brightness for the north ecliptic pole (NEP).
In [23]:
# Colina, Bohlin & Castelli solar spectrum normalised to V band flux
# of 184.2 ergs/s/cm^2/A,
solar_normalisation = 184.2 * u.erg * u.second**-1 * u.cm**-2 * u.Angstrom**-1
# Leinert at al NEP zodical light 1.81e-18 erg/s/cm^2/A/arcsec^2 at 0.5 um,
# Aldering suggests using 0.01 dex lower.
zl_nep = 1.81e-18 * u.erg * u.second**-1 * u.cm**-2 * u.Angstrom**-1 * \
u.arcsecond**-2 * 10**(-0.01)
zl_normalisation = zl_nep / solar_normalisation
lambda_c = 0.5 * u.micron
# Aldering reddening parameters
f_blue = 0.9
f_red = 0.48
def solar_to_ZL(waves):
"""
Calculate factors to convert a Colina et alsolar spectrum to spectrum
for the zodical light normalised for NEP. Follows the prescription of
Leinert et al (1997) with the revised parameters from Aldering (2001)
as used in the HST ETC (Giavalsico, Sahi & Bohlin (2002)).
"""
rfactor = np.where(waves < lambda_c, \
1.0 + f_blue * np.log(waves/lambda_c), \
1.0 + f_red * np.log(waves/lambda_c))
return zl_normalisation * rfactor
In [24]:
zl_waves = sun_waves.to(u.micron)
zl_flux = (sun_flux * solar_to_ZL(sun_waves))
zl_flux = zl_flux.to(u.Watt * u.m**-2 * u.arcsecond**-2 * u.micron**-1)
In [25]:
plt.semilogy(zl_waves, zl_flux)
plt.xlim(0.2,2.5)
plt.ylim(1e-18,3e-17)
plt.xlabel('Wavelength/$\mu$m')
plt.ylabel('Surface brightness / W m$^{-2}$ arcsec$^{-2}$ $\mu$m$^{-1}$')
plt.title('Zodical light at North Ecliptic Pole - energy flux surface brightness units');
In [26]:
photon_energies = (h * c / sun_waves) / u.photon
zl_photons = zl_flux / photon_energies
zl_photons = zl_photons.to(u.photon * u.second**-1 * u.m**-2 * u.arcsecond**-2 * u.micron**-1)
In [27]:
plt.semilogy(zl_waves, zl_photons, label='Zodical Light - NEP')
plt.xlim(0.2,2.5)
plt.ylim(2e-0,6e1)
plt.xlabel('Wavelength/$\mu$m')
plt.ylabel('Surface brightness / photons s$^{-1}$ m$^{-2}$ arcsec$^{-2}$ $\mu$m$^{-1}$')
plt.title('Zodical light at North Ecliptic Pole - photon flux surface brightness units');
In [28]:
zl_ep = []
for i, f in enumerate(filters):
# Interpolate etendue data at same wavelengths as zodiacal light data
eff_etendue = np.interp(zl_waves, f['Wavelength'], eff_etendues[i]) * \
eff_etendues[i].unit
# Multiply zodical light surface brightness by effective etendue
zl_epw = zl_photons * eff_etendue
# Integrate over wavelength
zl_ep.append(np.trapz(zl_epw, x=zl_waves))
zl_ep = u.Quantity(zl_ep)
zl_ep
Out[28]:
The zodiacal light background calculated so far is for the North Ecliptic Pole. Tables 16 & 17 of Leinert et al (1997) present data showing how the average surface brightness varies with (ecliptic longitude - solar ecliptic longitude) and ecliptic latitude. To allow evaluation of the zodiacal light surface brightness at arbitrary positions on the sky we will interpolate these tabular data and use it to scale the North Ecliptic Pole value.
The zodical light is generally very smooth, effects such as resonant structures, comet tails and asteroid collisions result in variations of less than 1% which are largely confined to low ecliptic latitudes. We will not consider these here.
Leinert et al also present evidence for variations in the degree of reddening with solar elongation however at visible wavelengths this effect will be small unless the telescope is pointed close to the Sun. As positions close to the Sun would not be observed by ASE we will neglect this effect.
In [29]:
# Data from table 17, Leinert et al (1997).
llsun = np.array([0, 5 ,10, 15, 20, 25, 30, 35, 40, 45, 60, 75, 90, 105, 120, 135, 150, 165, 180])
beta = np.array([0, 5, 10, 15, 20, 25, 30, 45, 60, 75])
zl_scale = np.array([[np.NaN, np.NaN, np.NaN, 3140, 1610, 985, 640, 275, 150, 100], \
[np.NaN, np.NaN, np.NaN, 2940, 1540, 945, 625, 271, 150, 100], \
[np.NaN, np.NaN, 4740, 2470, 1370, 865, 590, 264, 148, 100], \
[11500, 6780, 3440, 1860, 1110, 755, 525, 251, 146, 100], \
[6400, 4480, 2410, 1410, 910, 635, 454, 237, 141, 99], \
[3840, 2830, 1730, 1100, 749, 545, 410, 223, 136, 97], \
[2480, 1870, 1220, 845, 615, 467, 365, 207, 131, 95], \
[1650, 1270, 910, 680, 510, 397, 320, 193, 125, 93], \
[1180, 940, 700, 530, 416, 338, 282, 179, 120, 92], \
[910, 730, 555, 442, 356, 292, 250, 166, 116, 90], \
[505, 442, 352, 292, 243, 209, 183, 134, 104, 86], \
[338, 317, 269, 227, 196, 172, 151, 116, 93, 82], \
[259, 251, 225, 193, 166, 147, 132, 104, 86, 79], \
[212, 210, 197, 170, 150, 133, 119, 96, 82, 77], \
[188, 186, 177, 154, 138, 125, 113, 90, 77, 74], \
[179, 178, 166, 147, 134, 122, 110, 90, 77, 73], \
[179, 178, 165, 148, 137, 127, 116, 96, 79, 72], \
[196, 192, 179, 165, 151, 141, 131, 104, 82, 72], \
[230, 212, 195, 178, 163, 148, 134, 105, 83, 72]]).transpose()
# Normalise to a value of 1.0 at the NEP
zl_scale /= 77
# Expand range of angles to cover the full sphere by using symmetry
beta = np.array(np.concatenate((-np.flipud(beta)[:-1], beta))) * u.degree
llsun = np.array(np.concatenate((llsun, 360 - np.flipud(llsun)[1:-1]))) * u.degree
zl_scale = np.concatenate((np.flipud(zl_scale)[:-1],zl_scale))
zl_scale = np.concatenate((zl_scale,np.fliplr(zl_scale)[:,1:-1]),axis=1)
# Convert angles to radians within the required ranges.
beta = beta.to(u.radian).value + np.pi/2
llsun = llsun.to(u.radian).value
The Scipy spherical interpolation functions can't handle the missing values for positions near to the Sun, either the interpolation fails to run or there are ringing artefacts. To overcome this we use an initial cartesian interpolation pass on a subset of the data near to the sun in order to smoothly fill in the mising values, then use spherical interpolation on the full, filled in data set.
In [30]:
# For initial cartesian interpolation want the hole in the middle of the data,
# i.e. want to remap longitudes onto -180 to 180, not 0 to 360 degrees.
# Only want the region of closely spaced data point near to the Sun.
beta_c = beta[3:-3]
nl = len(llsun)
llsun_c = np.concatenate((llsun[nl//2+9:]-2*np.pi,llsun[:nl//2-8]))
zl_scale_c = np.concatenate((zl_scale[3:-3,nl//2+9:],zl_scale[3:-3,:nl//2-8]),axis=1)
# Convert everthing to 1D arrays (lists) of x, y and z coordinates.
llsuns, betas = np.meshgrid(llsun_c, beta_c)
llsuns_c = llsuns.ravel()
betas_c = betas.ravel()
zl_scale_cflat = zl_scale_c.ravel()
# Indices of the non-NaN points
good = np.where(np.isfinite(zl_scale_cflat))
# 2D cartesian interpolation function
zl_scale_c_interp = SmoothBivariateSpline(betas_c[good], llsuns_c[good], zl_scale_cflat[good])
# Calculate interpolated values
zl_scale_fill = zl_scale_c_interp(beta_c, llsun_c)
# Remap the interpolated values back to original ranges of coordinates
zl_patch = np.zeros(zl_scale.shape)
zl_patch[3:16,0:10] = zl_scale_fill[:,9:]
zl_patch[3:16,-9:] = zl_scale_fill[:,:9]
# Fill the hole in the original data with values from the cartesian interpolation
zl_patched = np.where(np.isfinite(zl_scale),zl_scale,zl_patch)
# Spherical interpolation function from the full, filled data set
zl_scale_interp = RectSphereBivariateSpline(beta, llsun, zl_patched, \
pole_continuity=(False,False), pole_values=(1.0, 1.0), \
pole_exact=True, pole_flat=False)
In [31]:
# Check the process with some plots
n_rows = 4
n_cols = 2
plt.subplot(n_rows,n_cols,1)
plt.imshow(zl_scale, norm=colors.LogNorm(vmin=0.7,vmax=300), \
interpolation='none', cmap='cubehelix_r')
plt.colorbar()
plt.title('Leinert et al data, normalised, full sphere')
plt.subplot(n_rows,n_cols,3)
plt.imshow(zl_scale_c, norm=colors.LogNorm(vmin=0.7,vmax=300), interpolation='none', cmap='cubehelix_r')
plt.colorbar()
plt.title('Centred and clipped')
plt.subplot(n_rows,n_cols,5)
plt.imshow(zl_scale_fill, norm=colors.LogNorm(vmin=0.7,vmax=300), interpolation='none', cmap='cubehelix_r')
plt.colorbar()
plt.title('Centred, clipped, cartesian interpolated')
plt.subplot(n_rows,n_cols,7)
plt.imshow(zl_patch, norm=colors.LogNorm(vmin=0.7,vmax=300), interpolation='none', cmap='cubehelix_r')
plt.colorbar()
plt.title('Cartesian interpolation remapped to original coords')
plt.subplot(n_rows,n_cols,2)
plt.imshow(zl_patched, norm=colors.LogNorm(vmin=0.7,vmax=300), interpolation='none', cmap='cubehelix_r')
plt.colorbar()
plt.title('Leinert et al data, filled with cartesian interpolation')
plt.subplot(n_rows,n_cols,4)
plt.imshow(zl_scale_interp(beta, llsun), norm=colors.LogNorm(vmin=0.7,vmax=300), \
interpolation='none', cmap='cubehelix_r')
plt.colorbar()
plt.title('Spherical interpolation evaluated at original points')
plt.subplot(n_rows,n_cols,6)
plt.imshow((zl_scale_interp(beta, llsun) - zl_scale)/zl_scale, interpolation='none',)
plt.colorbar()
plt.title('Relative deviation of interpolation from original')
plt.subplot(n_rows,n_cols,8)
plt.imshow(zl_scale_interp(np.linspace(0,np.pi,180),np.linspace(0,2*np.pi,360)), \
origin='lower', extent=[0,360,-90,90], norm=colors.LogNorm(vmin=0.7,vmax=300), \
interpolation='none', cmap='cubehelix_r')
plt.colorbar()
plt.title('Spherical interpolation')
plt.xlabel('$\lambda - \lambda_{\odot}$')
plt.ylabel('$\\beta$')
plt.tight_layout()
plt.gcf().set_size_inches(12,12)
We will precalculate an all-sky map of the zodical light surface brightness distribution using the HEALPix pixelation (http://healpix.sourceforge.net/). We can use the healpy
tools to extract projections of this map for signal to noise calculations.
Want to calculate these in galactic coordinates for later combination with dust maps
In [31]:
# Calculate all sky map
# *WARNING* This takes a long time, generally will just want to load it from zl_map.fits instead,
# see cell two below.
nside = 2048
zl_map=np.empty(hp.nside2npix(nside))
for i in range(zl_map.size):
theta, phi = hp.pix2ang(nside, i)
zl_map[i] = zl_scale_interp(theta, phi)
In [32]:
# Write the zodiacal light map to file (it takes a long time to calculate, don't want to have to redo it all the time).
hp.write_map('resources/zl_map.fits', zl_map, fits_IDL=False, coord='E', column_names=('ZL_rel_SB',))
In [34]:
# Read zodical light map from file
zl_map = hp.read_map('resources/zl_map.fits')
In [35]:
hp.mollview(zl_map, coord=('E', 'C'), cmap='cubehelix_r', min=0.7, max=10, norm='log', \
title='Zodiacal light surface brightness distribution, March equinox', \
unit='Surface brightness relative to NEP')
hp.graticule(coord='C', local=True)
hp.projscatter(0, 90, lonlat=True, marker='x', coord='E')
hp.projtext(0, 90, 'NEP', lonlat=True, coord='E')
hp.projscatter(0, -90, lonlat=True, marker='x', coord='E')
hp.projtext(0, -90, 'SEP', lonlat=True, coord='E')
plt.gcf().set_size_inches(12,7.5)
plt.savefig('ZL_March.png')
The main effect is simply due to the movement of the Sun around the ecliptic equator over the course of the year. As a result the distribution of zodiacal light on the sky rotates about an axis through the ecliptic poles.
In addition the zodiacal light exhibits seasonal variations due to the orbit of the Earth taking it above and below the symmetry plane of the zodiacal dust. At high ecliptic latitudes this manifests as an approximately 10% sinusoidal variation in surface brightness. For now we will ignore this.
Otherwise the zodiacal light is essentially time invariant. Resonant structures, comet tails and asteroid collisions result in variations of less than 1% and are largely confined to low ecliptic latitude.
In [36]:
# Very wasteful of RAM but for now just implement this using rotated copies of the ZL map
def map_rotate_polar(input_map, angle):
'''
Function to return a copy of a healpy map rotated about its polar axis
'''
npix = input_map.size
nside = hp.npix2nside(npix)
# Get coordinates of each pixel in original map
indices = np.arange(npix)
theta,phi = hp.pix2ang(nside, indices)
# New empty map for rotated copy
output_map = np.empty_like(input_map)
# Rotator object to transform pixel coordinates
rotator = hp.Rotator(deg=True, rot=[angle,0])
# For each pixel in original map calculate new coordinates and copy
# into appropriate pixel in new map.
output_map[hp.ang2pix(nside, *rotator(theta, phi))] = input_map
return output_map
In [37]:
zl_map_90 = map_rotate_polar(zl_map, 90)
zl_map_180 = map_rotate_polar(zl_map, 180)
zl_map_270 = map_rotate_polar(zl_map, 270)
In [38]:
hp.mollview(zl_map, coord=('E', 'C'), cmap='cubehelix_r', min=0.7, max=10, \
norm='log', sub=(2,2,1), cbar=False, \
title='Zodiacal light surface brightness distribution, March equinox', \
unit='Surface brightness relative to NEP')
hp.graticule(coord='C', local=True)
hp.projscatter(0, 90, lonlat=True, marker='x', coord='E')
hp.projtext(0, 90, 'NEP', lonlat=True, coord='E')
hp.projscatter(0, -90, lonlat=True, marker='x', coord='E')
hp.projtext(0, -90, 'SEP', lonlat=True, coord='E')
hp.mollview(zl_map_270, coord=('E', 'C'), cmap='cubehelix_r', min=0.7, max=10, \
norm='log', sub=(2,2,2), cbar=False, \
title='Zodiacal light surface brightness distribution, June solstice', \
unit='Surface brightness relative to NEP')
hp.graticule(coord='C', local=True)
hp.projscatter(0, 90, lonlat=True, marker='x', coord='E')
hp.projtext(0, 90, 'NEP', lonlat=True, coord='E')
hp.projscatter(0, -90, lonlat=True, marker='x', coord='E')
hp.projtext(0, -90, 'SEP', lonlat=True, coord='E')
hp.mollview(zl_map_180, coord=('E', 'C'), cmap='cubehelix_r', min=0.7, max=10, \
norm='log', sub=(2,2,3), cbar=False, \
title='Zodiacal light surface brightness distribution, September equinox', \
unit='Surface brightness relative to NEP')
hp.graticule(coord='C', local=True)
hp.projscatter(0, 90, lonlat=True, marker='x', coord='E')
hp.projtext(0, 90, 'NEP', lonlat=True, coord='E')
hp.projscatter(0, -90, lonlat=True, marker='x', coord='E')
hp.projtext(0, -90, 'SEP', lonlat=True, coord='E')
hp.mollview(zl_map_90, coord=('E', 'C'), cmap='cubehelix_r', min=0.7, max=10, \
norm='log', sub=(2,2,4), cbar=False, \
title='Zodiacal light surface brightness distribution, December solstice', \
unit='Surface brightness relative to NEP')
hp.graticule(coord='C', local=True)
hp.projscatter(0, 90, lonlat=True, marker='x', coord='E')
hp.projtext(0, 90, 'NEP', lonlat=True, coord='E')
hp.projscatter(0, -90, lonlat=True, marker='x', coord='E')
hp.projtext(0, -90, 'SEP', lonlat=True, coord='E')
plt.gcf().set_size_inches(12,7.5)
plt.savefig('ZL_seasonal.png')