-
Notifications
You must be signed in to change notification settings - Fork 0
/
SNR.py
152 lines (115 loc) · 4.71 KB
/
SNR.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
152
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
from scipy.interpolate import interp2d
from scipy.integrate import quad
import scipy.constants as const
import pandas as pd
from make_filters import make_filter, init_filters, init_filters_thomas
from scipy.optimize import dual_annealing, fsolve
from reflectance import get_ref
from comet import ref_rock, ref_ice
from camera import Camera
from unibe import *
read_out_noise = 5 # electrons
dark_current_noise = 20 # electrons / second
full_well_capacity = 33000 # electrons
peak_linear_charge = 27000 # electrons
def snr(s):
return np.sqrt(s)
# the above is a simplified version of:
# return s / np.sqrt(np.sqrt(s)**2 + read_out_noise**2 + dark_current_noise**2)
c1 = 595.16
w1 = 86.05
c2 = 690.00
w2 = 103.63
c3 = 788.60
w3 = 93.57
c4 = 899.94
w4 = 129.12
filter_center = 650
filter_width = 100
def get_mirror():
df_mirror = pd.read_csv("data/mirrors_transmission.txt", delimiter="\s")
M = interp1d(df_mirror.wavelength, df_mirror.transmission, fill_value="extrapolate")
# percent
return M
def get_detector():
df_qe = pd.read_csv("data/qe.txt", delimiter=",")
Q = interp1d(df_qe.Wavelength, df_qe.QE / 100, fill_value="extrapolate")
# electrons per photons
return Q
def get_solar():
df_solar = pd.read_csv("data/solar.csv", delimiter=";", skiprows=1)
S = interp1d(df_solar["Wavelength (nm)"], df_solar["Extraterrestrial W*m-2*nm-1"], fill_value="extrapolate")
# W per meter squared per nanometer
return S
if __name__ == "__main__":
M = get_mirror()
Q = get_detector()
S = get_solar()
print(snr(27000))
CoCa = Camera()
def integrand(w, alpha=0):
return w * M(w) * Q(w) * ref_rock(w, alpha).T * S(w)
phase_angle = np.arange(1, 90, 10)
snr_vals = []
snr_vals_80 = []
signals = []
signals_80 = []
centers = range(450, 1000, 50)
t_exp = 0.005
t = []
t_sat = []
df = pd.read_csv("data/texp.csv")
t10 = interp1d(df.alpha, df["texp10"], fill_value="extrapolate")
t80 = interp1d(df.alpha, df["texp80"], fill_value="extrapolate")
for alpha in phase_angle:
i = quad(integrand, filter_center - filter_width / 2, filter_center + filter_width / 2, args=(alpha))[0]
t_exp = t10(alpha)/1000
signal = CoCa.A_Omega / CoCa.G * t_exp * i / (const.h * const.c * CoCa.r_h ** 2) * 1e-9
snr_vals.append(snr(signal * CoCa.G))
signals.append(signal)
t_exp = t80(alpha)/1000
signal = CoCa.A_Omega / CoCa.G * t_exp * i / (const.h * const.c * CoCa.r_h ** 2) * 1e-9
snr_vals_80.append(snr(signal * CoCa.G))
signals_80.append(signal)
def func(t_exp):
i = quad(integrand, filter_center - filter_width / 2, filter_center + filter_width / 2, args=(alpha))[0]
signal = CoCa.A_Omega / CoCa.G * t_exp * i / (const.h * const.c * CoCa.r_h ** 2) * 1e-9
return snr(signal * CoCa.G) - 100
sol = fsolve(func, 0.001)
print(alpha, sol)
t.append(sol[0])
def func(t_exp):
i = quad(integrand, filter_center - filter_width / 2, filter_center + filter_width / 2, args=(alpha))[0]
signal = CoCa.A_Omega / CoCa.G * t_exp * i / (const.h * const.c * CoCa.r_h ** 2) * 1e-9
return snr(signal * CoCa.G) - 164.3
sol = fsolve(func, 0.001)
print(alpha, sol)
t_sat.append(sol[0])
wavelengths = np.linspace(300, 1100, 100)
fig, axes = plt.subplots(nrows=4)
axes[0].plot(phase_angle, np.array(t) * 1000, color=BLACK, ls="-", label="SNR 100")
axes[0].plot(phase_angle, np.array(t_sat) * 1000, color=RED, ls="-", label="saturation")
axes[0].set_xlabel("phase angle [°]")
axes[0].set_ylabel(r"$t_{exp}$ [ms]")
axes[0].legend()
axes[1].plot(phase_angle, snr_vals, color=RED, ls="-")
axes[1].plot(phase_angle, snr_vals_80, color=BLACK, ls="-")
axes[1].axhline(snr(CoCa.G * 2 ** 14), ls="--", color=RED, label="saturation")
axes[1].axhline(np.sqrt(full_well_capacity), ls="--", color=ORANGE, label="full well")
axes[1].axhline(np.sqrt(peak_linear_charge), ls="--", color=GREEN, label="peak linear")
axes[1].legend()
axes[1].set_xlabel("phase angle [°]")
axes[1].set_ylabel("SNR")
axes[2].plot(phase_angle, signals, color=RED, ls="-")
axes[2].plot(phase_angle, signals_80, color=BLACK, ls="-")
axes[2].axhline(2 ** 14, ls="--", color=RED)
axes[2].set_xlabel("phase angle [°]")
axes[2].set_ylabel("signal [DN]")
axes[3].plot(wavelengths, ref_rock(wavelengths, phase_angle).T, color=BLACK, ls="-")
axes[3].set_xlabel("wavelength [nm]")
axes[3].set_ylabel("I/F")
plt.savefig("plots/snr_new.png")
plt.show()