_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

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"]

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

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

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