-
Notifications
You must be signed in to change notification settings - Fork 1
/
cdim_sbsens.py
459 lines (359 loc) · 14.8 KB
/
cdim_sbsens.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
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
#########################################################################
##
## CDIM surface brightness sensitivty calculator
## Mike Zemcov ([email protected])
##
## Changelog:
## v0, Late 2016, baseline for NASA APROBES call.
## v1, June 2017, improved commenting for public distribution.
##
#########################################################################
#!/usr/bin/python
import sys
import numpy as np
import pylab as pl
import matplotlib.pyplot as plt
from scipy import constants as cst
from matplotlib.ticker import FormatStrFormatter
#from scipy.signal import savgol_filter
verbose = 2 # control verbosity
if verbose == 2:
fig=plt.figure(figsize=(6.5,5))
ax = fig.add_subplot(1,1,1)
tmpl = [0.8,1.,1.2,1.5,2.,2.5,3.,4.,5.,6.,7.5]
## Physical cst
c0 = cst.value('speed of light in vacuum')*1.e6 #um/s
h0 = cst.value('Planck constant') # J.s
kb = cst.value('Boltzmann constant') # J/K
hc = c0*h0 # um.J
Apert = 0.83 # Aperture in m
fnum = 4.5 # fnumber
R = 300 # lambda / delta lambda spectral resolution
lmax = 7.5 # maximum wavelength, microns
lmin = 0.75 # minumin wavelength, microns
Area_dp = 31 # deep survey area, sq deg
Area_sh = 311 # shallow survey area, sq deg
Tmission = 2. # mission length, years
obs_eff = 0.8 # observation efficiency (that is, fraction of time
# spent surveying)
t_dp = 104 # minutes per band per resolution element, deep survey
t_sh = 14.2 # minutes per band per resolution element, shallow survey
Pitch = 18. # detector pixel pitch, um
ZL_fac = 1.0 # multiplicative factor describing mean Zodiacal Light
# brightness above the NEP minimum
dQ = [10,10.5,10.5,21./np.sqrt(2)] # read noise, single sample, e-
T_samp = 1.5 # sample time, s
t_int = 250. # integration time, s
T_det = 35. # detector temperature, K
T_scope = 70. # telescope temperature, K
n_optics = 3 # number of optics in optical chain
eta_lvf = 0.80 # optical efficiency of lvf
blocking = 1e-5 #1.995e-5 #1e-5 # out of band blocking
blockingfile = 0 # do we use the blocking information from Omega?
if blocking == 0:
OD=0.
else:
OD = - np.log10(blocking)
Npix = 7 # number of pixels contributing to a source
pixelfac = 5335.67/R # the number of pixels that measure a single resolution
# element of width R
nativeR = pixelfac*R # native resolution per pixel column
pix_format = np.array([2048,2048]) # array format
fp_format = np.array([6,4]) # number and layout of detectors
# angle subtended by a single detector
th_pix = 3600.*180./np.pi*(Pitch*1.e-6)/(Apert*fnum)
if verbose:
print "Pixel angle is: " + str(th_pix) + " arcsec."
# number of pixels in the entire array
fp_pix = pix_format * fp_format
if verbose:
print "Focal plane is: " + str(fp_format[0]) + " detectors wide."
# this hard codes the type of detector (H2RG_NN) in each element of the array
fp_det_type = np.zeros(fp_pix[0])
fp_det_type[0:2048] = 1
fp_det_type[2048:3*2048] = 2
fp_det_type[3*2048:5*2048] = 3
fp_det_type[5*2048:6*2048] = 4
fp_cut_off = np.zeros(fp_pix[0])
fp_cut_off[0:2048] = 1.75
fp_cut_off[2048:3*2048] = 2.5
fp_cut_off[3*2048:5*2048] = 5.2
fp_cut_off[5*2048:6*2048] = 8.0
# set up the beginning wavelength
lam = np.zeros(fp_pix[0])
lam[0] = lmin
# loop through array asigning wavelength according to the
# pre-computed wavelength jump
for ipix in range(1,fp_pix[0]):
lam[ipix] = lam[ipix-1] * (1. + 1./nativeR)
# read in optical efficiency of mirrors
text_file = open('lookup/eta_au.txt')
rows = [[float(x) for x in line.split(',')[:]] for line in text_file]
cols = [list(col) for col in zip(*rows)]
text_file.close
# make optical efficiency
eta_au_l = np.asarray(cols[0])
eta_au_e = np.asarray(cols[1])
eta_au = np.interp(lam,eta_au_l,eta_au_e)
# read in optical efficiency of detector
text_file = open('lookup/eta_fpa.txt')
rows = [[float(x) for x in line.split(',')[:]] for line in text_file]
cols = [list(col) for col in zip(*rows)]
text_file.close
# make detective efficiency
eta_fpa_l = np.asarray(cols[0])
eta_fpa_e = np.asarray(cols[1])
eta_fpa = np.interp(lam,eta_fpa_l,eta_fpa_e)
# efficiency assumes
eta_opt = eta_au**n_optics
# build filter function
#filt_lvf = blocking * np.ones(fp_pix[0])
eta_lvf = eta_lvf * np.ones(fp_pix[0])
eta_tot = eta_opt * eta_fpa * eta_lvf
if verbose == 2:
ax.semilogx(lam,eta_opt,linestyle='-',label='AU Optics, 3 reflections')
ax.semilogx(lam,eta_lvf,linestyle='-',label='Dispersive Element')
ax.semilogx(lam,eta_fpa,linestyle='-',label='FPA QE')
ax.semilogx(lam,eta_tot,linestyle='-',marker='',\
label='Total System Efficiency',linewidth=2)
ax.set_xlabel(r'$\lambda$ ($\mu$m)')
ax.set_ylabel(r'Optical Efficiency')
ax.set_xlim([0.75,7.5])
ax.set_ylim([0,1])
ax.xaxis.set_ticks(tmpl)
ax.xaxis.set_major_formatter(FormatStrFormatter('%1.1f'))
plt.legend(loc=4)
plt.tight_layout()
plt.savefig('cdim_sbsens_eta_R'+str(R)+'.pdf')
# compute entendue of a detector
AOmega = 0.25*np.pi*(Apert**2)*(((th_pix/3600)*np.pi/180.)**2) #m^2/sr
if verbose:
print "Entendue is: " + str(AOmega) + " m^2 sr."
# generate an array of appropriate read noise
dQ_lam = np.zeros(fp_pix[0])
dQ_lam[np.where(fp_det_type == 1)] = dQ[0]
dQ_lam[np.where(fp_det_type == 2)] = dQ[1]
dQ_lam[np.where(fp_det_type == 3)] = dQ[2]
dQ_lam[np.where(fp_det_type == 4)] = dQ[3]
# convert to read noise appropriate to sample up the ramp
dQ_rin = dQ_lam*np.sqrt(6.*T_samp/t_int)
if verbose:
whsmpl = (np.abs((lam - 1.0)) == np.min(np.abs(lam - 1.0)))
print "Read noise at 1.0 microns is: " + str(dQ_rin[whsmpl]) + " e-/read."
# compute the sky background assuming two black bodies consistent
# with reflected solar and thermal ZL
sky_bkg = 0.85*6.7e3/lam**4/(np.exp(hc/(kb*5500.*lam))-1.) \
+ 0.85*5.120*1.e9/lam**4/(np.exp(hc/(kb*250.*lam))-1.) # in nW/m^2/sr
# include the "above minimum" factor
sky_bkg *= ZL_fac
if verbose:
print "Sky background at 1.0 microns is: " + str(sky_bkg[whsmpl]) + \
" nW/m^2/sr."
zlout = np.zeros(lam.size,dtype=[('var1',float),('var2',float)])
zlout['var1'] = lam
zlout['var2'] = sky_bkg
np.savetxt('cdim_sbsens_ZL_R'+str(R)+'.txt',zlout,fmt='%f, %f')
# compute the bakcground from the telescope
tele_bkg = 2. * hc * c0 * 1e-12 / ((1e-6*lam)**4*(np.exp(hc/(kb*T_scope*lam)) - 1)) * (1.-eta_au) * 1e9 # in nW/m^2/sr
if verbose:
print "Telescope background at 1.0 microns is: " + str(tele_bkg[whsmpl]) + \
" nW/m^2/sr."
# Compute the photocurrent from these contributions
i_sky = 1.e-9*sky_bkg*AOmega*eta_tot/(R*hc/lam) # e-/s
i_tele = 1.e-9*tele_bkg*np.pi*(Pitch*1e-6)**2/(R*hc/lam) # e-/s
# bug here: formally, I should be including out of band photocurrent
# in this calculation.
i_block = 0.* i_tele
i_pass = i_sky + i_tele
if blocking > 0 or blockingfile == 1:
print "Computing blocking..."
blocking = blocking * np.ones(np.size(lam))
if blockingfile==1:
# read in filter profiles
text_file = open('lookup/filters/filter_profiles.csv')
rows = [[float(x) for x in line.split(',')[:]] for line in text_file]
cols = [list(col) for col in zip(*rows)]
text_file.close
lambda_file = np.asarray(cols[0])/1000.
lambda_trans = np.asarray(cols[1])/100.
lambda_trans_filt = savgol_filter(lambda_trans,51,1)*0.3
blocking = np.interp(lam,lambda_file,lambda_trans_filt,right=0.0)
plt.clf()
ax = fig.add_subplot(1,1,1)
ax.loglog(lam,blocking)
plt.savefig('test.pdf')
i_pass = np.zeros(fp_pix[0])
i_block = np.zeros(fp_pix[0])
for ipix in range(1,fp_pix[0]):
filt = np.exp(-(lam - lam[ipix])**2 / \
((lam[ipix] / R)**2 / (4. * np.log(2))))
filt = filt / np.sum(filt) * eta_lvf
eta_filt = eta_opt * eta_fpa * filt
tp = np.where(lam <= fp_cut_off[ipix])
sumone = np.sum(sky_bkg[tp]*eta_filt[tp]*(lam[tp]/R))
sumtwo = np.sum(tele_bkg[tp]*eta_filt[tp]*(lam[tp]/R))
sumthr = np.sum(sky_bkg[tp]*blocking[tp]*(lam[tp]/R))
sumfor = np.sum(tele_bkg[tp]*blocking[tp]*(lam[tp]/R))
i_pass[ipix] = 1e-9 * AOmega / hc * (sumone + 1.8e-5**2 * sumtwo)
i_block[ipix] = 1e-9 * AOmega / hc * (sumthr + 1.8e-5**2 * sumfor)
blockingratio = i_pass/i_block
blockingratio[0] = blockingratio[1]
if verbose == 2:
plt.clf()
ax = fig.add_subplot(1,1,1)
ax.loglog(lam,blockingratio)
ax.loglog(lam,0*lam+1,linestyle='--',color='black')
ax.set_xlabel(r'$\lambda$ ($\mu$m)')
ax.set_ylabel(r'Ratio of In-band i to Out of Band i at OD%1.1f' % OD)
ax.set_xlim([0.75,7.5])
ax.set_ylim([0.1 * np.min(blockingratio),10*np.max(blockingratio)])
ax.xaxis.set_ticks(tmpl)
ax.xaxis.set_major_formatter(FormatStrFormatter('%1.1f'))
ax.set_xticklabels(['','',''],minor=True)
plt.tight_layout()
plt.savefig('cdim_sbsens_blocking_R'+str(R)+'.pdf')
# compute the total photocurrent
i_photo = i_pass + i_block
if verbose:
print "Photocurrent at 1.0 microns is: " + str(i_photo[whsmpl]) + \
" e-/s."
if verbose == 2:
plt.clf()
ax = fig.add_subplot(1,1,1)
#ax.loglog(lam,i_photo,linestyle='-',marker='',\
# label=r'Total $i_{\rm photo}$')
ax.loglog(lam,i_sky,linestyle='-',marker='',\
label=r'$i_{\rm photo}$ from ZL in $T_{\rm int}=250$s')
ax.loglog(lam,i_tele,linestyle='-',marker='',\
label=r'$i_{\rm photo}$ from 70K telescope')
#ax.loglog(lam,i_block,linestyle='-',marker='',\
# label=r'$i_{\rm photo}$ from filter leakage OD%1.1f' % OD)
ax.loglog(lam,i_block,linestyle='-',marker='',\
label=r'$i_{\rm photo}$ from out of band light')
# compute dark current - these come from empirical measurements
# taken from the literature
dc_one = 0;
dc_two = 1.8e14 * np.exp(-4280/T_det) + 1.1e9 * np.exp(-np.sqrt(57500/T_det)) + 1.2e-3
dc_three = 1.8e15 * np.exp(-2800/T_det) + 2.5e12 * np.exp(-np.sqrt(62000/T_det)) + 1.5e-3
dc_four = np.exp((T_det-42.5)/1.445) + 1.6e-1
# populate the array
DC = np.zeros(fp_pix[0])
DC[np.where(fp_det_type == 1)] = dc_one
DC[np.where(fp_det_type == 2)] = dc_two
DC[np.where(fp_det_type == 3)] = dc_three
DC[np.where(fp_det_type == 4)] = dc_four
if verbose == 2:
ax.loglog(lam,DC,linestyle='-',marker='',label=r'$i_{\rm dark}$ for 35 K detectors')
# total photocurrent measured by the detector
i_tot = i_photo + DC
if verbose == 2:
ax.loglog(lam,i_tot,linestyle='-',marker='',label=r'$i_{\rm total}$')
ax.set_xlabel(r'$\lambda$ ($\mu$m)',fontsize=15)
ax.set_ylabel(r'Current at Detector (e$^{-}$/s)',fontsize=15)
ax.set_xlim([0.75,7.5])
ax.set_ylim([1e-3,1e1])
ax.xaxis.set_ticks(tmpl)
ax.xaxis.set_major_formatter(FormatStrFormatter('%1.1f'))
ax.set_xticklabels(['','',''],minor=True)
plt.legend(loc=2,fontsize=12)
plt.tight_layout(pad=1.5)
plt.savefig('cdim_sbsens_iphoto_R'+str(R)+'.pdf')
# rms noise on the detector per integration
dnIn_ppix = np.sqrt(i_tot*t_int+dQ_rin**2)/t_int/i_sky*sky_bkg
dnIn_ppix_rn = np.sqrt(DC*t_int+dQ_rin**2)/t_int/i_sky*sky_bkg
dnIn_ppix_sky = np.sqrt(i_photo*t_int)/t_int/i_sky*sky_bkg
if verbose == 2:
plt.clf()
ax = fig.add_subplot(1,1,1)
ax.loglog(lam,dnIn_ppix_rn,linestyle='-',marker='',\
label='Detector Components')
ax.loglog(lam,dnIn_ppix_sky,linestyle='-',marker='',label='Photon Noise')
ax.loglog(lam,dnIn_ppix,linestyle='-',marker='',label='Total')
ax.set_xlabel(r'$\lambda$ ($\mu$m)')
ax.set_ylabel(r'$\delta \nu I_{\nu}$ / pixel / 250s (nW m$^{-2}$ sr$^{-1}$, $1 \sigma$)')
ax.set_xlim([0.75,7.5])
ax.set_ylim([1e0,1e4])
ax.xaxis.set_ticks(tmpl)
ax.xaxis.set_major_formatter(FormatStrFormatter('%1.1f'))
plt.legend(loc=3,fontsize=12)
plt.tight_layout()
plt.savefig('cdim_sbsens_dnIn_R'+str(R)+'.pdf')
if verbose == 2:
plt.clf()
ax = fig.add_subplot(1,1,1)
ax.semilogx(lam,dnIn_ppix_sky/dnIn_ppix_rn,linestyle='-',marker='')
ax.semilogx([0.5,7.5],[1,1],marker='',linestyle=':',color='black')
ax.set_xlabel(r'$\lambda$ ($\mu$m)')
ax.set_ylabel(r'Ratio of photon noise to detector noise')
ax.set_xlim([0.75,7.5])
ax.set_ylim([0,3.1])
ax.xaxis.set_ticks(tmpl)
ax.xaxis.set_major_formatter(FormatStrFormatter('%1.1f'))
plt.tight_layout()
plt.savefig('cdim_sbsens_ratio_R'+str(R)+'.pdf')
if verbose:
print "Surface brightness noise at 1.0 microns is: " + \
str(dnIn_ppix[whsmpl]) + " nW/m^2/sr."
# compute flux uncertainty in uJy
dF = 1.e-9*1.e26*1.e6*((np.pi/180.)*(th_pix/3600.))**2*dnIn_ppix*(lam/c0) * np.sqrt(Npix) # uJy
dF_single = dF
n_reps = 1./(4.*6.5)
dF = dF * np.sqrt(n_reps) #t_int / (obs_eff * 60) / t_dp)
# assume some number of sigma
Nsig = 5
# convert to Mab
Mab = -2.5*np.log10(Nsig*1.e-6*dF/3631.)
Mab_single = -2.5*np.log10(Nsig*1e-6*dF_single/3631.)
# convert to line flux
Sline = 1e18*Nsig * 3e-14 * dF / R / (lam * 1e4) * 1e7 * 1e-4
if verbose == 2:
plt.clf()
ax = fig.add_subplot(1,1,1)
ax.loglog(lam,dF,linestyle='-',marker='',label="Deep Survey")
ax.loglog(lam,dF_single,linestyle='-',marker='',label="Single Integration")
ax.set_xlabel(r'$\lambda$ ($\mu$m)')
ax.set_ylabel(r'$\delta F$ ($\mu$Jy, 1$\sigma$)')
ax.set_xlim([0.75,7.5])
ax.set_ylim([0.2,20])
ax.xaxis.set_ticks(tmpl)
ax.xaxis.set_major_formatter(FormatStrFormatter('%1.1f'))
ax.legend()
plt.tight_layout()
plt.savefig('cdim_sbsens_dF_R'+str(R)+'.pdf')
if verbose == 2:
plt.clf()
ax = fig.add_subplot(1,1,1)
ax.semilogx(lam,Mab,linestyle='-',marker='',label="Deep Survey")
ax.semilogx(lam,Mab_single,linestyle='-',marker='',label="Single Integration")
ax.set_xlabel(r'$\lambda$ ($\mu$m)')
ax.set_ylabel(r'$M_{\rm AB}$ (5$\sigma$)')
ax.set_xlim([0.75,7.5])
ax.set_ylim([19,25])
ax.xaxis.set_ticks(tmpl)
ax.xaxis.set_major_formatter(FormatStrFormatter('%1.1f'))
ax.legend()
plt.tight_layout()
plt.savefig('cdim_sbsens_Mab_R'+str(R)+'.pdf')
if verbose == 2:
plt.clf()
ax = fig.add_subplot(1,1,1)
ax.semilogx(lam,Sline,linestyle='-',marker='')
ax.set_xlabel(r'$\lambda$ ($\mu$m)')
ax.set_ylabel(r'Line Sensitivity (1$\times 10^{-18}$ erg s$^{-1}$ cm$^{-2}$, 1$\sigma$)')
ax.set_xlim([0.75,7.5])
ax.set_ylim([0,20])
ax.xaxis.set_ticks(tmpl)
ax.xaxis.set_major_formatter(FormatStrFormatter('%1.1f'))
plt.tight_layout()
plt.savefig('cdim_sbsens_Sline_R'+str(R)+'.pdf')
if verbose == 2:
plt.close()
# save the result of our labors
dout = np.zeros(lam.size,dtype=[('var1',float),('var2',float),('var3',float),('var4',float),('var5',float)])
dout['var1'] = lam
dout['var2'] = dnIn_ppix
dout['var3'] = dF
dout['var4'] = Mab
dout['var5'] = Sline
np.savetxt('cdim_sbsens_out_R'+str(R)+'.txt',dout,fmt='%f, %f, %f, %f, %f')
### return