Skip to content
GitLab
Switch branch/tag
Find file Normal viewHistoryPermalink
_util.py 3.47 KB
Newer Older
Fang Yuedong's avatar
Baseline  
Fang Yuedong committed
1
2 
import numpy as np
import os
xin's avatar
syn master code    
xin committed
3
4
5
6
7
8
9
10
11 
import math
from pylab import *
from scipy import interpolate

try:
    import importlib.resources as pkg_resources
except ImportError:
    # Try backported to PY<37 'importlib_resources'
    import importlib_resources as pkg_resources
Fang Yuedong's avatar
Baseline  
Fang Yuedong committed
12
Fang Yuedong's avatar
v0 for testing    
Fang Yuedong committed
13
14
15 
VC_A = 2.99792458e+18  # speed of light: A/s
VC_M = 2.99792458e+8   # speed of light: m/s
H_PLANK = 6.626196e-27 # Plank constant: erg s
Fang Yuedong's avatar
Baseline  
Fang Yuedong committed
16
17
18
19
20
21
22
23
24 

def photonEnergy(lambd):
	""" The energy of photon at a given wavelength

	Parameter:
		lambd: the wavelength in unit of Angstrom
	Return:
		eph: energy of photon in unit of erg
	"""
Fang Yuedong's avatar
v0 for testing    
Fang Yuedong committed
25
26 
	nu = VC_A / lambd
	eph = H_PLANK * nu
xin's avatar
syn master code    
xin committed
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 
	return eph

'''
description: 
param {*} aperture: unit m, default 2 m
param {*} psf_fwhm: psf fwhm, default 0.1969"
param {*} pixelSize: pixel size, default 0.074"
param {*} pmRation: the ratio of souce flux in the limit mag calculation
param {*} throughputFn: throuput file name
param {*} readout: unit, e-/pixel
param {*} skyFn: sky sed file name, average of hst, 'sky_emiss_hubble_50_50_A.dat'
param {*} darknoise: unit, e-/pixel/s
param {*} exTime: exposure time one time, default 150s
param {*} exNum: exposure number, defautl 1
param {*} fw, full well value( or saturation value),default 90000e-/pixel
return {*} limit mag and saturation mag
'''

def calculateLimitMag(aperture = 2.0, psf_fwhm = 0.1969,pixelSize = 0.074, pmRation = 0.8, throughputFn = 'i_throughput.txt', readout = 5.0, skyFn= 'sky_emiss_hubble_50_50_A.dat', darknoise = 0.02,exTime = 150, exNum = 1, fw = 90000):
	with pkg_resources.path('ObservationSim.Instrument.data.throughputs', throughputFn) as data_file:
		throughput_f = np.loadtxt(data_file)
	thr_i = interpolate.interp1d(throughput_f[:,0]/10, throughput_f[:,1]); # wavelength in anstrom
	f_s = 200
	f_e = 1100
	delt_f = 0.5

	data_num = int((f_e-f_s)/delt_f+1)

	eff = np.zeros([data_num,2])
	eff[:,0] = np.arange(f_s,f_e+delt_f,delt_f)
	eff[:,1] = thr_i(eff[:,0])

	wave = np.arange(f_s,f_e+delt_f,delt_f)
	wavey = np.ones(wave.shape[0])

	with pkg_resources.path('ObservationSim.Instrument.data.throughputs', skyFn) as data_file:
		skydata = np.loadtxt(data_file)
	skydatai = interpolate.interp1d(skydata[:,0]/10, skydata[:,1]*10)

	sky_data = np.zeros([data_num,2])
	sky_data[:,0] = np.arange(f_s,f_e+delt_f,delt_f)
	sky_data[:,1] = skydatai(sky_data[:,0])

	flux_sky = trapz((sky_data[:,1])*eff[:,1],sky_data[:,0])
	skyPix = flux_sky*pixelSize*pixelSize*pi*(aperture*aperture/4)


	###limit mag

	r_pix = psf_fwhm*0.7618080243778568/pixelSize # radius RE80, pixel
	cnum = math.pi * r_pix * r_pix
	sn = 5

	d = skyPix*exTime*exNum*cnum + darknoise*exTime*exNum*cnum+readout*readout*cnum*exNum
	a=1
	b=-sn*sn
	c=-sn*sn*d

	flux = (-b+sqrt(b*b-4*a*c))/(2*a)/pmRation
	limitMag = -2.5*log10(flux/(54799275581.04437 * trapz(wavey*eff[:,1]/wave,wave, 0.1)*exTime*exNum*pi*(aperture/2)*(aperture/2)))


	### saturation mag

	from astropy.modeling.models import Gaussian2D
	m_size = int(20 * psf_fwhm/pixelSize)
	if m_size%2 == 0:
		m_size + 1
	m_cen = m_size//2
	psf_sigma = psf_fwhm/2.355/pixelSize

	gaussShape = Gaussian2D(1, m_cen, m_cen, psf_sigma, psf_sigma)
	yp, xp = np.mgrid[0:m_size, 0:m_size]
	psfMap = gaussShape(xp, yp)
	maxRatio = np.amax(psfMap)/np.sum(psfMap)
	# print(maxRatio)

	flux_sat = fw/maxRatio*exNum
	satMag = -2.5*log10(flux_sat/(54799275581.04437 * trapz(wavey*eff[:,1]/wave,wave, 0.1)*exTime*exNum*pi*(aperture/2)*(aperture/2)));


	return limitMag , satMag