In [1]:
import numpy as np
from scipy import stats
import astropy.units as u
from astropy.coordinates import SkyCoord, get_sun, GeocentricTrueEcliptic, Angle
from astropy.wcs import WCS
from astropy.time import Time
from astropy.table import Table
from astropy.visualization import hist
import matplotlib.pyplot as plt
from matplotlib import rcParams
from matplotlib import colors

In [2]:
from imp import reload

In [3]:
%matplotlib inline

In [4]:
rcParams['figure.figsize'] = 14, 9

In [7]:
import astroimsim

In [167]:
reload(astroimsim)


Out[167]:
<module 'astroimsim' from '/mnt/data/ajh/Documents/Spaaaaaaaaaaaaaaaaaaaaaace/HuntsmanInSpace/Simulations/astroimsim/astroimsim.py'>

Zodiacal Light

ZodicalLight class constructs a zodical light model based on that used by the HST ETC.

Initialise zodical light model


In [25]:
zl = astroimsim.ZodiacalLight()

Absolute spectral photon flux density surface brightness


In [26]:
plt.semilogy(zl.waves, zl.photon_sfd, 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');


Relative surface brightness as a function of sky position.


In [27]:
plt.imshow(zl._spatial(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('Spatial variation of Zodiacal Light')
plt.xlabel('$\lambda - \lambda_{\odot}$')
plt.ylabel('$\\beta$')
plt.gcf().set_size_inches(12, 5)


Imager

Imager class represents an imaging instrument

Set up imager parameters

Entrance aperture area (astropy Quantity)

Here we assume an unobstructed circular aperture of 90 mm diameter.


In [28]:
pupil_area = (np.pi * (90 * u.mm)**2 / 4).si
pupil_area


Out[28]:
$0.0063617251 \; \mathrm{m^{2}}$

Optical throughput (astropy Table)

The Table should have 'Wavelength' and 'Throughput' columns with specified units. In real instruments there will be some degree of field dependance, that is not implemented yet.

A very simplistic estimate, neglect bulk absorption and assume a constant 99.5% transmission at each vacuum-glass surface.


In [53]:
n_surfaces = 7
wavelengths = np.array((250, 1050)) * u.nm
t = np.array((0.995**n_surfaces, 0.995**n_surfaces)) * u.dimensionless_unscaled
throughput = Table(names = ('Wavelength', 'Throughput'), data = (wavelengths, t))
throughput


Out[53]:
<Table length=2>
WavelengthThroughput
nm
float64float64
250.00.965520646809
1050.00.965520646809

Filter transmission functions (Python dictionary of astropy Tables)

Keys of the dictionary should be filter name strings, the Tables should have 'Wavelength' and 'Transmission' columns with specified units.

Here we use nominal Space Eye filter profiles derived from Chebyshev Type I filter functions of order 50, ripple factor 0.14, and peak transmission set to 0.95. These have properties similar to a real optical bandpass filter, e.g. they fairly closely resemble the DECam $i'$ and $z'$ filters but with somewhat sharper edges (facilitated by the slower focal ratio).


In [30]:
N = 50
r = 0.14

filters = {"i1":Table((zl.waves, astroimsim.cheby_band(zl.waves, 700*u.nm, 850.0 * u.nm, N, ripple=r)), \
                      names = ['Wavelength', 'Transmission']), \
           "i2":Table((zl.waves, astroimsim.cheby_band(zl.waves, 700*u.nm, 855.5 * u.nm, N, ripple=r)), \
                      names = ['Wavelength', 'Transmission']), \
           "i3":Table((zl.waves, astroimsim.cheby_band(zl.waves, 700*u.nm, 861.0 * u.nm, N, ripple=r)), \
                      names = ['Wavelength', 'Transmission']), \
           "z1":Table((zl.waves, astroimsim.cheby_band(zl.waves, 854.0*u.nm, 1000 * u.nm, N, ripple=r)), \
                      names = ['Wavelength', 'Transmission']), \
           "z2":Table((zl.waves, astroimsim.cheby_band(zl.waves, 859.5*u.nm, 1000 * u.nm, N, ripple=r)), \
                      names = ['Wavelength', 'Transmission']), \
           "z3":Table((zl.waves, astroimsim.cheby_band(zl.waves, 865.0*u.nm, 1000 * u.nm, N, ripple=r)), \
                      names = ['Wavelength', 'Transmission'])}

for f in filters.values():
    # Set units correctly
    f['Wavelength'] = f['Wavelength'].to(u.um)
    f['Wavelength'].unit = u.um
    f['Transmission'].unit = u.dimensionless_unscaled

In [31]:
plt.subplot(2,1,1)

plt.plot(zl.waves, zl.photon_sfd/zl.photon_sfd.max(), 'c-', label='Zodiacal light')
plt.plot(filters['i2']['Wavelength'], filters['i2']['Transmission'], 'b-')
plt.plot(filters['z2']['Wavelength'], filters['z2']['Transmission'], 'r-')

plt.xlim(0.65,1.05)
plt.xlabel('Wavelength/$\mu$m')
plt.ylabel('Transmission')
plt.title('Australian Space Eye filter profiles')

plt.subplot(2,1,2)
plt.plot(filters['i1']['Wavelength'], filters['i1']['Transmission'], 'b--', label='$i1$')
plt.plot(filters['i2']['Wavelength'], filters['i2']['Transmission'], 'b-', label='$i2$')
plt.plot(filters['i3']['Wavelength'], filters['i3']['Transmission'], 'b-.', label='$i3$')
plt.plot(filters['z1']['Wavelength'], filters['z1']['Transmission'], 'r--', label='$z1$')
plt.plot(filters['z2']['Wavelength'], filters['z2']['Transmission'], 'r-', label='$z2$')
plt.plot(filters['z3']['Wavelength'], filters['z3']['Transmission'], 'r-.', label='$z3$')

plt.plot(zl.waves, zl.photon_sfd/zl.photon_sfd.max(), 'c-', label='Zodiacal light')
plt.xlim(0.835,0.88)
plt.xlabel('Wavelength/$\mu$m')
plt.ylabel('Transmission')  
plt.legend(loc='best')
plt.gcf().set_size_inches(14, 18)

#plt.gcf().set_size_inches(8, 12)
#plt.savefig('filters.pdf')


Image sensor quantum efficiency (astropy Table)

The Table should have 'Wavelength' and 'QE' columns with specified units.

Here we use QE data for the e2v CIS115 image sensor at -40$^\circ$C extracted from a plot in Soman et el (2014).


In [32]:
qe_m40C = Table.read('../resources/cis115_QE_-40C.csv')
qe_m40C['Wavelength'].unit = u.um
qe_m40C['QE'].unit = u.electron / u.photon

In [33]:
plt.plot(qe_m40C['Wavelength'], qe_m40C['QE'], 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)')


Out[33]:
<matplotlib.text.Text at 0x7f7cdb84f160>

Image sensor gain


In [34]:
gain = 1.0 * u.electron / u.adu # Well below read noise, easily accommodates dynamic range in 16 bits

Read noise


In [35]:
read_noise = 4.0 * u.electron

Dark current

At visible wavelengths the dominant source of instrumental background will be the dark current of the image sensor. Wang et al (2014), http://proceedings.spiedigitallibrary.org/proceeding.aspx?articleid=1891414, measured the dark current of an e2v CMOS image sensor prototype (CIS107) between 250 and 300 K and found that the dark current was well fit by $D = 10^{0.0488T - 12.772}$ e$^-$ s$^{-1}$ pixel$^{-1}$, where $T$ is the temperature in Kelvin.


In [134]:
Ts = np.arange(250,301,10) * u.Kelvin
Ts2 = np.arange(220,251,10) * u.Kelvin
DSs = 10**(0.0488 * Ts.value - 12.772) * u.electron * u.second**-1 * u.pixel**-1
DSs2 = 10**(0.0488 * Ts2.value - 12.772) * u.electron * u.second**-1 * u.pixel**-1
plt.semilogy(Ts.to(u.deg_C, equivalencies=u.temperature()), DSs)
plt.semilogy(Ts2.to(u.deg_C, equivalencies=u.temperature()), DSs2, 'b--')
plt.xlim(-55,20)
plt.ylim(0.01,40)
plt.xlabel('Temperature / $^\circ$C')
plt.ylabel('Dark current / e$^-$ s$^{-1}$ pixel$^{-1}$')
plt.title('e2v CIS-115 Dark Current Model')
#plt.savefig('CIS115DC.png')
plt.gcf().set_size_inches(8, 6)
plt.savefig('CIS115DC.pdf')


For our nominal operating temperature of -40$^\circ$C we get


In [37]:
temperature = (-40 * u.deg_C).to(u.K, equivalencies=u.temperature())
dark_current = 10**(0.0488 * temperature.value - 12.772) * u.electron * u.second**-1 * u.pixel**-1
dark_current


Out[37]:
$0.040338524 \; \mathrm{\frac{e^{-}}{pix\,s}}$

Create an Imager instance to represent the instrument

This will result in a template WCS being created and effective aperture areas (aperture area x optical throughout x image sensor QE) as a function of wavelength being precalculated for each filter.


In [165]:
ase = astroimsim.Imager(2000, 1504, 3 * u.arcsecond, pupil_area, throughput, filters, \
                        qe_m40C, gain, read_noise, temperature, zl)

In [142]:
ase.wcs


Out[142]:
WCS Keywords

Number of WCS axes: 2
CTYPE : 'RA---TAN'  'DEC--TAN'  
CRVAL : 0.0  0.0  
CRPIX : 1000.5  752.5  
PC1_1 PC1_2  : 1.0  0.0  
PC2_1 PC2_2  : 0.0  1.0  
CDELT : 0.00083333333333333328  0.00083333333333333328  
NAXIS    : 2000 1504

Check effective areas


In [143]:
ase._eff_areas["i1"]['Effective Area'].unit


Out[143]:
$\mathrm{\frac{m^{2}\,e^{-}}{ph}}$

In [144]:
plt.plot(ase._eff_areas["i2"]['Wavelength'], ase._eff_areas["i2"]['Effective Area'], label="$i2$")
plt.plot(ase._eff_areas["z2"]['Wavelength'], ase._eff_areas["z2"]['Effective Area'], label="$z2$")
plt.xlim(0.65, 1.05)
plt.ylim(0,0.005)
plt.xlabel('Wavelength/$\mu$m')
plt.ylabel('Effective aperture area/$m^2 e^- \gamma^{-1}$')
plt.legend()
plt.title('Australian Space Eye Effective Aperture Area')


Out[144]:
<matplotlib.text.Text at 0x7f7cd987a080>

Check observed Zodiacal light normalisation


In [145]:
ase._zl_ep


Out[145]:
{'i1': <Quantity 0.03081828402512441 electron / (arcsec2 s)>,
 'i2': <Quantity 0.03143722550020661 electron / (arcsec2 s)>,
 'i3': <Quantity 0.03205334706260966 electron / (arcsec2 s)>,
 'z1': <Quantity 0.008787292484020487 electron / (arcsec2 s)>,
 'z2': <Quantity 0.008170885335179342 electron / (arcsec2 s)>,
 'z3': <Quantity 0.007547432329835031 electron / (arcsec2 s)>}

Ecliptic pole observed Zodical light photoelectrons per second based on nominal pixel scale


In [146]:
for f in sorted(ase.filters):
    print(f, ase._zl_ep[f] * ase.pixel_scale**2)


i1 0.2773645562261197 electron / s
i2 0.2829350295018595 electron / s
i3 0.28848012356348696 electron / s
z1 0.07908563235618439 electron / s
z2 0.07353796801661408 electron / s
z3 0.06792689096851527 electron / s

In [147]:
plt.bar([1,2,3,4,5, 6], [(ase._zl_ep[f] * ase.pixel_scale**2).value for f in sorted(ase.filters)], \
        align='center', color=('b','b','b','r','r', 'r'))
plt.hlines(ase.dark_current.value, 0.5, 6.5, linestyle='dashed', label='CIS115 dark current at $-40^\circ$C')
plt.xlim(0.5,6.5)
plt.legend(loc='best')
plt.xticks((1,2,3,4,5,6), sorted(ase.filters.keys()))
plt.ylabel('Rate / e$^-$ s$^{-1}$ pixel$^{-1}$')
plt.title('Photoelectron rate from Zodiacal light, North Ecliptic Pole')


Out[147]:
<matplotlib.text.Text at 0x7f7cd9667f98>

SNR in faint source limit $\propto \frac{1}{\sqrt{St + Dt + r^2}}$

Relative SNR (relative to zero instrumental noise) $= (1 + \frac{D + r^2/t}{S})^{-1/2}$


In [148]:
def rel_snr(t, T, f):
    S = (ase._zl_ep[f] * ase.pixel_scale**2).value
    D = 10**(0.0488 * T.to(u.Kelvin, u.temperature()).value - 12.772)
    return (1 + (D + ase.read_noise.value**2/t)/S)**-0.5

In [149]:
ts = np.arange(0,901,10)
plt.plot(ts, rel_snr(ts, -45 * u.Celsius, 'i2'), 'b', label='$i2$, $-45^\circ$C')
plt.plot(ts, rel_snr(ts, -40 * u.Celsius, 'i2'), 'm', label='$i2$, $-40^\circ$C')
plt.plot(ts, rel_snr(ts, -35 * u.Celsius, 'i2'), 'r', label='$i2$, $-35^\circ$C')
plt.plot(ts, rel_snr(ts, -30 * u.Celsius, 'i2'), 'y', label='$i2$, $-30^\circ$C')
plt.plot(ts, rel_snr(ts, -45 * u.Celsius, 'z2'), 'b--', label='$z2$, $-45^\circ$C')
plt.plot(ts, rel_snr(ts, -40 * u.Celsius, 'z2'), 'm--', label='$z2$, $-40^\circ$C')
plt.plot(ts, rel_snr(ts, -35 * u.Celsius, 'z2'), 'r--', label='$z2$, $-35^\circ$C')
plt.plot(ts, rel_snr(ts, -30 * u.Celsius, 'z2'), 'y--', label='$z2$, $-30^\circ$C')
plt.ylim(0.5,1)
plt.legend(loc='best', framealpha=0.5)
plt.ylabel('Relative signal to noise ratio')
plt.xlabel('Sub-exposure time / second')
plt.title('SNR relative to SNR without instrumental noise')
plt.gcf().set_size_inches(8, 6)
plt.savefig('relsnr.pdf')


/mnt/data/ajh/Documents/virtualenvs/python3.4.3/lib/python3.4/site-packages/ipykernel/__main__.py:4: RuntimeWarning: divide by zero encountered in true_divide

Pivot wavelengths


In [153]:
ase._pivot_waves


Out[153]:
{'i1': <Quantity 0.7654012921160522 um>,
 'i2': <Quantity 0.7672586251427996 um>,
 'i3': <Quantity 0.7690646460678219 um>,
 'z1': <Quantity 0.9040549967735372 um>,
 'z2': <Quantity 0.9079801742794213 um>,
 'z3': <Quantity 0.9119022171348669 um>}

Sensitivity integrals


In [168]:
ase._sensitivities


Out[168]:
{'i1': <Quantity 2.190917978179763e+21 electron m um2 / J>,
 'i2': <Quantity 2.2457058612753334e+21 electron m um2 / J>,
 'i3': <Quantity 2.2978265474639208e+21 electron m um2 / J>,
 'z1': <Quantity 8.102300778919353e+20 electron m um2 / J>,
 'z2': <Quantity 7.555676571915522e+20 electron m um2 / J>,
 'z3': <Quantity 7.028721133436624e+20 electron m um2 / J>}

With sensitivity integrals and ecliptic pole Zodical light count rates can do a quick sensitivity limit calculation

SNR, $\sigma = \frac{F_\lambda S_\lambda t}{\sqrt{{F_\lambda S_\lambda t + Zt + Dt + r^2 t/t_{sub}}}}$,

rearranging and solving for $F_\lambda$,

$F_\lambda = \frac{\sigma}{S_\lambda \sqrt{t}} \left(\frac{\sigma}{2 \sqrt{t}} + \sqrt{\frac{\sigma^2}{4 t} + \left( Z + D + r^2/t_{sub} \right)} \right)$


In [234]:
def limit(t, t_sub, sigma, f):
    n = ase._zl_ep[f] * ase.pixel_scale**2 + ase.dark_current + read_noise**2 / (t_sub* u.electron)
    st = sigma / t**0.5 * (u.electron)**0.5
    return ((st / ase._sensitivities[f]) * (st + (st**2 + n)**0.5)).to(u.Watt / (u.m**2 * u.micron))

def limit_sb(t, t_sub, sigma, f):
    return (limit(t, t_sub, sigma, f) / ase.pixel_scale**2).to(u.Watt / (u.m**2 * u.arcsecond**2 * u.micron))

def limit_sum_sb(t, t_sub, sigma, fs):
    z = u.Quantity([ase._zl_ep[f] for f in fs]).sum() * ase.pixel_scale**2
    n = z + ase.dark_current + read_noise**2 / (t_sub* u.electron)
    st = sigma / t**0.5 * (u.electron)**0.5
    sl = u.Quantity([ase._sensitivities[f] for f in fs]).sum()
    return ((st / sl) * (st + (st**2 + n)**0.5) / ase.pixel_scale**2).to(u.Watt / (u.m**2 * u.arcsecond**2 * u.micron))

Total exposure time per filter per hemisphere for a 2 year mission, assuming 3-4 exposures per orbit.


In [247]:
(np.array((3, 4)) * 600 * u.second / (96.7 * u.minute) * u.year / 6).to(u.second)


Out[247]:
$[1631727,~2175636] \; \mathrm{s}$

Surface brightness limits for individual filters at end of 2 year mission.


In [246]:
print(limit_sb(2e6 * u.second, 600 * u.second, 1.0, 'i2'))


2.0720759776160347e-20 W / (arcsec2 m2 micron)

In [245]:
print(limit_sb(2e6 * u.second, 600 * u.second, 1.0, 'z2'))


3.905648359083874e-20 W / (arcsec2 m2 micron)

Surface brightness limits for sum of $i'$, $z'$ filters at end of 2 year mission.


In [250]:
print(limit_sum_sb(2e6 * u.second, 600 * u.second, 1.0, ['i1', 'i2', 'i3']))


1.1172695845296046e-20 W / (arcsec2 m2 micron)

In [251]:
print(limit_sum_sb(2e6 * u.second, 600 * u.second, 1.0, ['z1', 'z2', 'z3']))


1.8595370899333165e-20 W / (arcsec2 m2 micron)

Sensitivity limits in "mag/arcsec^2" for astronomers.


In [257]:
i = limit_sum_sb(2e6 * u.second, 600 * u.second, 1.0, ['i1', 'i2', 'i3']) * u.arcsec**2
i = i.to(u.Jansky, equivalencies=u.spectral_density(ase._pivot_waves['i2']))
i.to(u.ABmag)


Out[257]:
$30.546948 \; \mathrm{mag}$$\mathrm{\left( \mathrm{AB} \right)}$

In [258]:
z = limit_sum_sb(2e6 * u.second, 600 * u.second, 1.0, ['z1', 'z2', 'z3']) * u.arcsec**2
z = z.to(u.Jansky, equivalencies=u.spectral_density(ase._pivot_waves['z2']))
z.to(u.ABmag)


Out[258]:
$29.628158 \; \mathrm{mag}$$\mathrm{\left( \mathrm{AB} \right)}$

Plot!


In [253]:
ts = 600 * 10**np.arange(0,4.01,0.1) * u.second
t_sub = 600 * u.second
sigma = 1.0
plt.loglog(ts, limit_sb(ts, t_sub, sigma, 'i1'), 'b--', label='$i1$')
plt.loglog(ts, limit_sb(ts, t_sub, sigma, 'i2'), 'b-', label='$i2$')
plt.loglog(ts, limit_sb(ts, t_sub, sigma, 'i3'), 'b-.', label='$i3$')
plt.loglog(ts, limit_sum_sb(ts, t_sub, sigma, ['i1', 'i2', 'i3']), 'c-', label="$i'$ sum")
plt.loglog(ts, limit_sb(ts, 600 * u.second, 1.0, 'z1'), 'r--', label='z1')
plt.loglog(ts, limit_sb(ts, 600 * u.second, 1.0, 'z2'), 'r-', label='z2')
plt.loglog(ts, limit_sb(ts, 600 * u.second, 1.0, 'z3'), 'r-.', label='z3')
plt.loglog(ts, limit_sum_sb(ts, t_sub, sigma, ['z1', 'z2', 'z3']), 'm-', label="$z'$ sum")
plt.xlim(600.0, 2e6)
plt.ylim(1e-20, 3e-18)
plt.xlabel('Total exposure time per filter / second')
plt.ylabel('Surface brightness / W m$^{-2}$ arcsec$^{-2}$ $\mu$m$^{-1}$')
plt.legend(loc='best')
plt.title('Sensitivity limit, ecliptic pole, no extinction, SNR = 1 per pixel, $t_{\mathrm{sub}}=600$ s')
plt.gcf().set_size_inches(8, 6)
plt.savefig('sens.pdf')


Check dark frame


In [23]:
ase.dark_current


Out[23]:
$0.040338524 \; \mathrm{\frac{e^{-}}{s}}$

In [24]:
plt.imshow(ase.dark_frame.value, interpolation='none', aspect='equal', cmap='gray')
plt.colorbar()


Out[24]:
<matplotlib.colorbar.Colorbar at 0x7fdf0a25a470>