_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('observation_sim.instruments.data.throughputs').joinpath(throughputFn) as data_file:
            throughput_f = np.loadtxt(data_file)
    except AttributeError:
        with pkg_resources.path('observation_sim.instruments.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('observation_sim.instruments.data.throughputs').joinpath(skyFn) as data_file:
            skydata = np.loadtxt(data_file)
    except AttributeError:
        with pkg_resources.path('observation_sim.instruments.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