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]:
In [25]:
zl = astroimsim.ZodiacalLight()
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');
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)
In [28]:
pupil_area = (np.pi * (90 * u.mm)**2 / 4).si
pupil_area
Out[28]:
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]:
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')
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]:
In [34]:
gain = 1.0 * u.electron / u.adu # Well below read noise, easily accommodates dynamic range in 16 bits
In [35]:
read_noise = 4.0 * u.electron
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]:
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]:
In [143]:
ase._eff_areas["i1"]['Effective Area'].unit
Out[143]:
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]:
In [145]:
ase._zl_ep
Out[145]:
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)
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]:
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')
In [153]:
ase._pivot_waves
Out[153]:
In [168]:
ase._sensitivities
Out[168]:
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]:
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'))
In [245]:
print(limit_sb(2e6 * u.second, 600 * u.second, 1.0, 'z2'))
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']))
In [251]:
print(limit_sum_sb(2e6 * u.second, 600 * u.second, 1.0, ['z1', 'z2', 'z3']))
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]:
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]:
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')
In [23]:
ase.dark_current
Out[23]:
In [24]:
plt.imshow(ase.dark_frame.value, interpolation='none', aspect='equal', cmap='gray')
plt.colorbar()
Out[24]:
In [25]:
hist(ase.dark_frame.value.ravel(), bins=np.arange(0,0.2,0.005));
In [26]:
centre = '16h00m00s +57d30m00s'
In [27]:
t = Time.now()
In [28]:
noiseless_i1 = ase.make_noiseless_image(centre, t, "i1", zl=zl)
noiseless_i2 = ase.make_noiseless_image(centre, t, "i2", zl=zl)
noiseless_i3 = ase.make_noiseless_image(centre, t, "i3", zl=zl)
noiseless_z1 = ase.make_noiseless_image(centre, t, "z1", zl=zl)
noiseless_z2 = ase.make_noiseless_image(centre, t, "z2", zl=zl)
noiseless_z3 = ase.make_noiseless_image(centre, t, "z3", zl=zl)
In [29]:
plt.imshow(noiseless_i2, interpolation='none', aspect='equal', cmap='cubehelix')
plt.colorbar()
Out[29]:
In [30]:
real_i2 = ase.make_image_real(noiseless_i2, exp_time = 600 * u.second)
real_i2_sub = ase.make_image_real(noiseless_i2, exp_time = 600 * u.second, subtract_dark=True)
In [31]:
plt.imshow(real_i2, interpolation='none', aspect='equal', cmap='gray', vmax=290, vmin=100)
plt.colorbar()
Out[31]:
In [32]:
plt.imshow(real_i2_sub, interpolation='none', aspect='equal', cmap='gray', vmax=290, vmin=100)
plt.colorbar()
Out[32]:
In [33]:
hist(real_i2.data.ravel(), bins=np.arange(100,290));
hist(real_i2_sub.data.ravel(), bins=np.arange(100,290));
plt.xlim(100,290)
Out[33]:
In [ ]:
In [ ]:
In [54]:
import math
In [55]:
def theta_max(F):
return math.atan(1/(2*F))
In [56]:
def lambda_shift(theta, n_eff):
return math.sqrt(1 - (math.sin(theta)/n_eff)**2)
In [57]:
def NB_filter(wave, wave_c, FWHM, w):
if abs(wave - wave_c) > 0.5*(1 + w)*FWHM:
return 0.0
elif abs(wave - wave_c) < 0.5*(1-w)*FWHM:
return 1.0
else:
return 1 - (abs(wave - wave_c) - 0.5*(1-w)*FWHM)/(w*FWHM)
In [58]:
def fast_NB_filter(wave, wave_c, FWHM, w, F, n_eff, n_thetas):
thetas = np.linspace(0, theta_max(F), n_thetas)
num = math.fsum([NB_filter(wave/lambda_shift(theta, n_eff), wave_c, FWHM, w)*math.sin(theta) for theta in thetas])
denom = math.fsum(math.sin(theta) for theta in thetas)
return num/denom
In [59]:
481/90
Out[59]:
In [60]:
waves = zl.waves.value
rc = (0.71 + 0.859) / 2.0
rf = 0.859 - 0.71
rw = 0.059
In [ ]:
In [ ]: