[35]:

# from astropy.table import Table, QTable, join, vstack
from astropy.table import Table, QTable

# from astropy.io import fits
from astropy import units as u

import os

# import matplotlib
import matplotlib.pyplot as plt

from scipy.integrate import trapezoid
# from scipy.ndimage import uniform_filter1d
# from scipy import interpolate

import numpy as np
# import pandas as pd

import tqdm

import pickle

from svnds import PATH_DATA

The numbers in this notebook is only for forecast!

Survey Plan & Design

  • Site

    • Chile, El Sauce Observatory

  • Surveys

    • Reference (RIS): 20,000 deg^2

    • Wide-Field (Wide)

    • Intensive Monitoring Survey (IMS)

  • Details

    • Moon phases

    • Weather loss

[36]:

###### Telescope
D = 50.8             # Aperture size [cm]
D_obscuration = 29.8   # Central Obscuration (diameter) [cm]
EFL = 1537.3           # Effective FL [mm]

Deff = np.sqrt(D**2 - D_obscuration**2)

###### Detector
array = 'CMOS'       # detector array type - verified, websete
dQ_RN = 3.           # [e], readout noise - estimated, websete
I_dark = 0.01        # [e/s], dark current - seems to be very low ???
pixel_size = 3.76    # [um], "pitch" - verified, websete
theta_pixel = pixel_size / EFL * 206.265  # [arcsec], pixel scale- estimated, websete
nxpix, nypix = 9576, 6388  # [pixels], detector format, approx. 9k x 6k - verified, websete

[37]:

survey_area_per_night = 113.  # [deg^2]
hrs_per_night = 9.            # hours per night

FOV_per_pointing = nxpix*nypix * theta_pixel**2 / 3600**2
T_exposure = hrs_per_night / (survey_area_per_night/FOV_per_pointing) * 3600 / 2

print(f'FOV           = {FOV_per_pointing:10.3g} deg^2')
print(f'Exposure time = {T_exposure:10.3g} sec')

FOV           =        1.2 deg^2
Exposure time =        172 sec

[38]:

fwhm_seeing = 1.5     # [arcsec]
FWHM0 = fwhm_seeing   # analysis aperture size
Tsamp = 180.          # individual exposure time [sec], 3min

# How many pixels does a point source occupy?
# Effective number of pixels for a Gaussian PSF with FWHM0
Npix_ptsrc = np.pi*(FWHM0/theta_pixel)**2
print(f'number of pixels for a point source = {Npix_ptsrc:.1f}')

number of pixels for a point source = 27.8

[39]:

# Constants
h = 6.626e-27 # erg/Hz
c = 3e10      # cm/s

path_skytable = os.path.join(PATH_DATA, "systematics/", "skytable.fits")
sky_tbl = Table.read(path_skytable)

wl_nm = sky_tbl['lam']          # nm
wl_um = wl_nm / 1e3       # micron
wl_cm = wl_um / 1e4       # cm
wl_am = wl_angstrom = wl_nm * 10  # angstrom
nu = 3e18 / wl_angstrom   # Hz

I_lambda = sky_tbl['flux']      # [ph/s/m2/micron/arcsec2] photon reate
f_lambda = I_lambda * (h*c/wl_cm) / (1e2**2) / (1e4)  # erg/s/cm2/A

f_nu = f_lambda * wl_angstrom * (wl_cm/c) / (1e-23 * 1e-6)  # micro Jansky

#####
# with open('filters_corrected', 'rb') as fr:
#         filters_corrected = pickle.load(fr)

Calculate 7DS sensitivity

[40]:

def synth_phot(wave, flux, wave_lvf, resp_lvf, tol=1e-3, return_photonrate = False):
    """
    Quick synthetic photometry routine.

    Parameters
    ----------
    wave : `numpy.ndarray`
        wavelength of input spectrum.
    flux : `numpy.ndarray`
        flux density of input spectrum in f_nu unit
        if `return_countrate` = True, erg/s/cm2/Hz is assumed
    wave_lvf : `numpy.ndarray`
        wavelength of the response function
    resp_lvf : `numpy.ndarray`
        response function. assume that this is a QE.
    tol : float, optional
        Consider only wavelength range above this tolerence (peak * tol).
        The default is 1e-3.

    Returns
    -------
    synthethic flux density in the input unit
        if return_photonrate = True, photon rates [ph/s/cm2]

    """
    index_filt, = np.where(resp_lvf > resp_lvf.max()*tol)

    index_flux, = np.where(np.logical_and( wave > wave_lvf[index_filt].min(),
                                           wave < wave_lvf[index_filt].max() ))

    wave_resamp = np.concatenate( (wave[index_flux], wave_lvf[index_filt]) )
    wave_resamp.sort()
    wave_resamp = np.unique(wave_resamp)
    flux_resamp = np.interp(wave_resamp, wave, flux)
    resp_resamp = np.interp(wave_resamp, wave_lvf, resp_lvf)

    if return_photonrate:
        h_planck = 6.626e-27 # erg/Hz
        return trapezoid(resp_resamp / wave_resamp * flux_resamp, wave_resamp) / h_planck

    return trapezoid(resp_resamp / wave_resamp * flux_resamp, wave_resamp) \
         / trapezoid(resp_resamp / wave_resamp, wave_resamp)

[41]:

def pointsrc_signal2noise(mag_src, Tsamp, wave_sys, resp_sys):
    """
    Calculate SN for a point source

    Input
        mag_src: AB mag of the source, scalar
        Tsamp: individual exposure time [sec], can be scalar or array

    WARNING: !!! ALL VARIABLES ARE GLOBALLY DECLARED !!!
    """
    Naper = Npix_ptsrc

    #source & sky background
    f_nu_src = f_nu*0 + 10**(-0.4*(mag_src + 48.6))  # erg/s/cm2/Hz
    f_nu_sky = f_nu*(1e-23*1e-6)                     # erg/s/cm2/Hz/arcsec2

    photon_rate_src = synth_phot(wl_um, f_nu_src, wave_sys, resp_sys, return_photonrate=True)  # ph/s/cm2
    photon_rate_sky = synth_phot(wl_um, f_nu_sky, wave_sys, resp_sys, return_photonrate=True)  # ph/s/cm2/arcsec2

    I_photo_src = photon_rate_src * (np.pi/4*Deff**2)                     # [e/s] per aperture (no aperture loss)
    I_photo_sky = photon_rate_sky * (np.pi/4*Deff**2) * (theta_pixel**2)  # [e/s] per pixel

    Q_photo_src = I_photo_src * Tsamp
    Q_photo_sky = I_photo_sky * Tsamp
    Q_photo_dark = I_dark * Tsamp

    SN = Q_photo_src / np.sqrt(Q_photo_src + Naper*Q_photo_sky + Naper*Q_photo_dark + Naper*dQ_RN**2)

    return SN

[42]:

lambda_7ds = np.arange(4000., 9000., 125)
Tsamp = 180 * 364/14 * 0.5 * 7

with open('filters_corrected', 'rb') as fr:
    filters_corrected = pickle.load(fr)

[43]:

unit_SB  = u.nW/(u.m)**2/u.sr
unit_cntrate = u.electron / u.s

T_sens = (QTable(
             names=('band', 'wavelength', 'I_photo_sky', 'mag_sky', 'mag_pts'),
             dtype=(np.int16,float,float,float,float,) )
         )
for key in T_sens.colnames:
    T_sens[key].info.format = '.4g'

for iw, wl_cen in enumerate(lambda_7ds):

    wl_cen = int(wl_cen)
    wave_sys = filters_corrected['wave_' + str(wl_cen)]
    resp_sys = filters_corrected['resp_' + str(wl_cen)]

    # photon rate
    photon_rate = synth_phot(wl_um, f_nu*(1e-23*1e-6), wave_sys, resp_sys, return_photonrate=True)

    # SB
    SB_sky = synth_phot(wl_um, f_nu*(1e-23*1e-6), wave_sys, resp_sys)

    # sky photo-current or count rate [e/s]
    I_photo = photon_rate * (np.pi/4*Deff**2) * (theta_pixel**2)

    # noise in count per obs [e].
    Q_photo = (I_photo+I_dark)*Tsamp
    dQ_photo = np.sqrt(Q_photo)

    # noise in count rate [e/s]
    # read-noise (indistinguishable from signal) should be added
    dI_photo = np.sqrt(dQ_photo**2 + dQ_RN**2)/Tsamp

    # surface brightness limit [one pixel]
    dSB_sky = (dI_photo/I_photo)*SB_sky
    mag_sky = -2.5*np.log10(5*dSB_sky) - 48.60

    # point source limit
    dFnu = np.sqrt(Npix_ptsrc) * dSB_sky*(theta_pixel)**2
    mag_pts = -2.5*np.log10(5*dFnu) - 48.60

    # Add data to the table
    T_sens.add_row([iw, wl_cen, I_photo, mag_sky, mag_pts])

[44]:

T_sens

[44]:
QTable length=40
bandwavelengthI_photo_skymag_skymag_pts
int16float64float64float64float64
040000.153423.3323.02
141250.205223.3723.05
242500.222123.4623.14
343750.255723.4823.16
445000.299923.4623.14
546250.322123.4723.15
647500.334223.4723.15
748750.339923.4823.16
850000.335823.4823.16
951250.329423.4723.15
...............
3077500.126221.9221.6
3178750.143821.7721.45
3280000.108121.7621.44
3381250.0628621.8521.53
3482500.0750121.6221.3
3583750.105821.3421.02
3685000.0819321.3421.03
3786250.0637921.3321.02
3887500.0932821.0320.71
3988750.106720.8220.5
[34]:

# path_save = os.path.join(PATH_DATA, "sensitivity/" 'SDS_WFS_7yr_eff_5sig.csv')
# T_sens.write(path_save, overwrite = True)

[ ]: