_util.py 4.35 KB
Newer Older
Fang Yuedong's avatar
Fang Yuedong committed
1
2
import numpy as np
import os
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
Fang Yuedong committed
12

Fang Yuedong's avatar
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
Fang Yuedong committed
16

17
18
19
20
ALL_FILTERS = ["NUV","u", "g", "r", "i","z","y","GU", "GV", "GI", "FGS"]
PHOT_FILTERS = ["NUV", "u", "g", 'r', 'i', 'z', 'y', 'FGS']
SPEC_FILTERS = ["GI", "GV", "GU"]

Fang Yuedong's avatar
Fang Yuedong committed
21
22
23
24
25
26
27
28
29
30
31
def rotate_conterclockwise(x0, y0, x, y, angle):
	"""
    Rotate a point counterclockwise by a given angle around a given origin.

    The angle should be given in radians.
    """
	angle = np.deg2rad(angle)
	qx = x0 + np.cos(angle)*(x - x0) - np.sin(angle) * (y - y0)
	qy = y0 + np.sin(angle)*(x - x0) + np.cos(angle) * (y - y0)
	return qx, qy

Fang Yuedong's avatar
Fang Yuedong committed
32
33
34
35
36
37
38
39
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
Fang Yuedong committed
40
41
	nu = VC_A / lambd
	eph = H_PLANK * nu
42
43
44
	return eph

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):
Fang Yuedong's avatar
Fang Yuedong committed
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
	'''
	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
	'''
60
61
62
63
64
65
	try:
		with pkg_resources.files('ObservationSim.Instrument.data.throughputs').joinpath(throughputFn) as data_file:
			throughput_f = np.loadtxt(data_file)
	except AttributeError:
		with pkg_resources.path('ObservationSim.Instrument.data.throughputs', throughputFn) as data_file:
			throughput_f = np.loadtxt(data_file)
66
67
68
69
70
71
72
73
74
75
76
77
78
79
	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])

80
81
82
83
84
85
	try:
		with pkg_resources.files('ObservationSim.Instrument.data.throughputs').joinpath(skyFn) as data_file:
			skydata = np.loadtxt(data_file)
	except AttributeError:
		with pkg_resources.path('ObservationSim.Instrument.data.throughputs', skyFn) as data_file:
			skydata = np.loadtxt(data_file)
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
	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)
124
	# print(maxRatio)
125
126
127
128
129
130

	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