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