-
Notifications
You must be signed in to change notification settings - Fork 0
/
exptime.py
executable file
·151 lines (135 loc) · 8.59 KB
/
exptime.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
#! /usr/bin/python
import numpy as np
from astropy import units as u
from astropy import constants as c
import rvrms
'''
Exposure time calculator for a list of target stars, desired measurement quality, and optical system specifications. Currently a proof of concept derived from HARPS, but capable to some degree of considering stellar properties. Note: output is in seconds!
Makes simplifications, and but should be valid in optical (400-650 nm), red (650-1000 nm) and near IR (1000-2500 nm).
Guesses at exposure time based on saturating the detector, then scales with t^-0.5 to get the precision. Also considers readout time.
'''
# Beatty spectra (simulated): 100 A chunks; Solar metallicity;
# R = 100k (4000-6500 A), 165000 (6500-10000 A), 88000 (10000A - 25000 A)
# In general need to specify wavelength range, as corrections vary with the one chosen (optical, red, near IR).
# Table from http://www.personal.psu.edu/tgb15/beattygaudi/table1.dat
beatty = np.genfromtxt("table1.dat", dtype=None)
beatty.dtype.names = ('Angstrom', 'Teff', 'Uncertaintykms')
def λ_peak(Teff, λ_min=0 * u.angstrom, λ_max=1e13 * u.angstrom):
λ = (c.b_wien / (Teff * u.K)).si # Star is assumed to be a perfect blackbody
# Correcting for if peak wavelength is outside of detector, angstroms assumed
if (λ > λ_max):
λ = λ_max
if (λ < λ_min):
λ = λ_min
return λ
def extinction(λ):
'''Returns atmospheric extinction coefficient (in magnitudes/airmass), given a wavelength (n Angstroms
best guess, based on measurements at Texas A&M Univeristy, and CFHT in Mauka Kea
http://www.gemini.edu/sciops/telescopes-and-sites/observing-condition-constraints/extinction
http://adsabs.harvard.edu/abs/1994IAPPP..57...12S
Values should not be considered trustworthy outside of ~3000-10000 Angstroms, and really 3500-9500 at that.
Rayleigh scattering, and a general scattering/absorption (handwaving together aerosols and ozone) are considered, but the Water vapor, etc lines in the IR are not!'''
return (0.09 + (3080.0/λ)**4)
def airmass(zenith_angle, site_elevation=0, scale_height=8400):
'''Returns the 'airmass' using the sec(z) approximation, and considering reduced atmospheric pressure from being at altitude. Inputs: zenith_angle (radians), site_elevation (meters, default 0), scale_height (default 8400 m)'''
return np.exp(-site_elevation/scale_height)/np.cos(zenith_angle)
def time_guess(Teff, FeH, logg, vsini, theta_rot, rstar, dstar, v_mac, atmo, efficiency, area, R, gain, read_noise, dark_current, n_pix, λ_peak, well_depth):
'''Generating an initial guess on exposure time. We assume a saturated detector is appropriate, thus putting one in the photon noise limited regime.'''
# We have spectra available every 100 K for 1700-7000 K, but only every 200 for 7000-9800.
temp = int(round(Teff,-2))
if temp > 7000:
temp = int(2*round(Teff/2,-2))
BTSettl = np.genfromtxt(str(temp), dtype=float)
# BT Settl spectra are labled by Teff, and available every 100 K.
# These spectra have a wavelength (Angstroms), and a flux (1 erg/s/cm^2/Angstrom) column
# sum of flux(λ)*Δλ(λ)/1e8 == total flux (in power/area) emitted.
# Multiply bt stellar_radius^2/distance^2 for recieved instead of emitted.
# Spectra downloaded from other sources will use different units!
photons = 0
#Counting up photons for SNR -> exposure time calc
#adding up photons in 100 A chunk centered on peak wavelength:
for i in np.arange(0, len(BTSettl)-1):
if((BTSettl[i][0]>= λ_peak/u.angstrom-50) and (BTSettl[i][0] <= λ_peak/u.angstrom+15)):
power = (BTSettl[i][1]*u.erg/u.cm**2/u.s/u.angstrom) * ((BTSettl[i+1][0]-BTSettl[i][0]) * u.angstrom)
power *= rstar**2/dstar**2
opacity = np.exp(-extinction(BTSettl[i][0]) * atmo)
power *= opacity #rescaling because of atmospheric scattering/absorption.
power = power.si
photons += (power * BTSettl[i][0] * u.angstrom / (c.h * c.c)) * area * efficiency
photons /= (100*u.angstrom / (λ_peak/R)) #100 Angstrom range -> per resolution element
photons /= n_pix #per resolution element -> per pixel
#print("photons/s/pixel", photons)
time_guess = well_depth / photons
SNR_sat = well_depth * n_pix * gain / np.sqrt(well_depth * n_pix * gain + n_pix**2 * ((gain*read_noise*2.2*2)**2 + (gain * dark_current * time_guess)**2))
# compare with I_0[λ][3] = I_0[λ][2]/pix*gain / np.sqrt(I_0[λ][2]/pix*gain + (gain*read_noise*2.2*2*n_expose)**2 + (gain*dark_current*exptime)**2)
#simpler alternate SNR
#SNR_sat = np.sqrt(well_depth * n_pix)
#print("Estimated exposure time", time_guess, time_guess.si)
return time_guess.si, SNR_sat
def time_actual(sigma_v, Teff, FeH, logg, vsini, theta_rot, rstar, dstar, v_mac, atmo, exptime, efficiency, area, R, gain, read_noise, dark_current, n_pix, λ_min, λ_max, Δλ, read_time, t_min, SNR, SNR_sat):
RV_actual, SNR_actual = rvrms.rvcalc(Teff, FeH, logg, vsini, theta_rot, rstar, dstar, v_mac, atmo, exptime, efficiency, area, R, gain, read_noise, dark_current, n_pix, λ_min, λ_max, Δλ, 2)
print("full well rv, snr", RV_actual, SNR_actual)
# from time_guess SNR_actual should equal SNR_sat, so we don't use it
time_actual = exptime * (RV_actual/sigma_v)**2
time_actual_snr = exptime * (SNR/SNR_sat)**2
if (time_actual_snr > time_actual):
print("SNR requirement exceeds RV info requirement")
time_actual = time_actual_snr
else:
print("RV info requirement exceeds SNR requirement")
number_exposures = np.ceil(time_actual/exptime)
# Effectively redoing the exposure time to always between 0.5 and 1 saturation times.
time_exposure = time_actual/number_exposures
reads = read_time * number_exposures
if (time_actual+(number_exposures-1)*read_time < t_min):
print("Warning: nominal exposure time is less than minimum time.")
# Increasing the number of exposures so that the time between the
# start of the first and end of the last is long enough to average
# out over the p-mode oscillations.
number_exposures = np.ceil((t_min+read_time)/(time_exposure+read_time))
time_actual = time_exposure*number_exposures
reads = read_time * number_exposures
RV_final, SNR_final = rvrms.rvcalc(Teff, FeH, logg, vsini, theta_rot, rstar, dstar, v_mac, atmo, time_actual, efficiency, area, R, gain, read_noise, dark_current, n_pix, λ_min, λ_max, Δλ, 2)
return time_actual, reads, time_exposure, number_exposures, SNR_actual, RV_actual, SNR_final, RV_final
if __name__ == '__main__':
# HARPS / ESO 3.6 m at La Silla
telescope = 3.566 * u.m
#area = np.pi * (telescope/2)**2 # telescope area, central obstruction ignored
area = 8.8564 * u.m * u.m
efficiency = 0.02 #could be ~1-3%, depending on what is measured
λ_min = 3800 * u.angstrom
λ_max = 6800 * u.angstrom
Δλ = 100 * u.angstrom
R = 110000 #110k or 120k, depending on source
gain = 1/1.42 # ADU/e-, assume 1:1 photon to electron generation
read_noise = 7.07 #RMS of ± spurious electrons
dark_current = 4 / u.hour
n_pix = 4
well_depth = 30000 #electrons (or photons)
read_time = 25 * u.s
# Alpha Cen B
Teff = 5214 # 5214 K
logg = 4.37
FeH = 0.23
rstar = 0.865 * u.solRad
dstar = 1.3 * u.pc
vsini = 2*np.pi*rstar/(38.7*u.day) * u.s/u.km
theta_rot = 1.13*vsini.si
v_mac = np.float('nan') # km/s, but exact value unclear
sigma_v = 5e-5 # target single measurement photon noise precision in km/s
t_min = 0 * u.s # minimum exposure time in seconds. Typically 300 to even out p-mode oscillations
SNR = 300
#671 is supposed nominal, breakeven with RV is between 935 and 930
# Atmo
zenith_angle = 0 # Can be read in from stellar observations!
scale_height = 8400 * u.m #Can make arguments for `8500 m to ~7500 m, but the higher temperatures typical of the lower atmosphere are more relevant here.
site_elevation = 2400 * u.m #ESO 3.6 m telescope, La Silla, Chile
atmo = airmass(zenith_angle, site_elevation, scale_height)
λ_peak = λ_peak(Teff, λ_min, λ_max)
print("Exposure time test ("+str(sigma_v)+" km/s target precision) with HARPS on Alpha Cen B. 76 seconds on-sky expected.")
guess, SNR_sat = time_guess(Teff, FeH, logg, vsini, theta_rot, rstar, dstar, v_mac, atmo, efficiency, area, R, gain, read_noise, dark_current, n_pix, λ_peak, well_depth)
actual, readout, exposure, number, SNR_actual, RV_actual, SNR_final, RV_final = time_actual(sigma_v, Teff, FeH, logg, vsini, theta_rot, rstar, dstar, v_mac, atmo, guess, efficiency, area, R, gain, read_noise, dark_current, n_pix, λ_min, λ_max, Δλ, read_time, t_min, SNR, SNR_sat)
print("guess, actual exposure, readout(s), total, 'real' individual exposure, number of exposures")
print(guess, actual, readout, actual+readout, exposure, number)
print("SNR in, SNR sat, SNR out", SNR, SNR_sat, SNR_actual)
print("RV_out, SNR_final (estimate), RV_final (estimate)", RV_actual, SNR_final, RV_final)