_util.py

import numpy as np
import os
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

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

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
	"""
	nu = VC_A / lambd
	eph = H_PLANK * nu
	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):
	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)
	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])

	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)
	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