Star.py

import galsim
import os, sys
import numpy as np
import astropy.constants as cons
from astropy.table import Table
from scipy import interpolate

from ObservationSim.MockObject._util import integrate_sed_bandpass, getNormFactorForSpecWithABMAG, getObservedSED, getABMAG, tag_sed
from ObservationSim.MockObject.MockObject import MockObject

class Star(MockObject):
    def __init__(self, param):
        super().__init__(param)

    def unload_SED(self):
        """(Test) free up SED memory
        """
        del self.sed

    def load_SED(self, survey_type, normFilter=None, target_filt=None, sed_lib=None, sed_path=None):
        if survey_type == "photometric":
            norm_thr_rang_ids = normFilter['SENSITIVITY'] > 0.001
            # spec_lambda = Table.read(sed_path, path=f"/SED/wave_{self.model_tag}")
            # spec_flux = Table.read(sed_path, path=f"/SED/{self.sed_type}")
            # wave, flux = spec_lambda['col0'].data, spec_flux['col0'].data
            _, wave, flux = tag_sed(
                h5file=sed_lib, 
                model_tag=self.param['model_tag'], 
                teff=self.param['teff'],
                logg=self.param['logg'],
                feh=self.param['feh'])
            flux_photon = flux * (wave / (cons.h.value * cons.c.value)) * 1e-13
            sed_photon = Table(np.array([wave, flux_photon]).T, names=('WAVELENGTH', 'FLUX'))
            # Get scaling factor for SED
            sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=self.param['mag_use_normal'],
                spectrum=sed_photon,
                norm_thr=normFilter,
                sWave=np.floor(normFilter[norm_thr_rang_ids][0][0]),
                eWave=np.ceil(normFilter[norm_thr_rang_ids][-1][0]))
            sed_photon = np.array([sed_photon['WAVELENGTH'], sed_photon['FLUX']*sedNormFactor]).T
            # Convert to galsim.SED object
            spec = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(sed_photon[:, 1]), interpolant='nearest')
            self.sed = galsim.SED(spec, wave_type='A', flux_type='1', fast=False)
            # Get magnitude
            interFlux = integrate_sed_bandpass(sed=self.sed, bandpass=target_filt.bandpass_full)
            self.param['mag_%s'%target_filt.filter_type] = getABMAG(
                interFlux=interFlux, 
                bandpass=target_filt.bandpass_full)
            # print('mag_use_normal = ', self.param['mag_use_normal'])
            # print('mag_%s = '%target_filt.filter_type, self.param['mag_%s'%target_filt.filter_type])

        elif survey_type == "spectroscopic":
            # self.sedPhotons(sed_path=sed_path)
            self.sedPhotons(sed_lib=sed_lib)

    def sedPhotons(self, sed_path=None, sed_lib=None):
        # spec_lambda = Table.read(sed_path, path=f"/SED/wave_{self.model_tag}")
        # spec_flux = Table.read(sed_path, path=f"/SED/{self.sed_type}")
        _, w_orig, f_orig = tag_sed(
                h5file=sed_lib, 
                model_tag=self.param['model_tag'], 
                teff=self.param['teff'],
                logg=self.param['logg'],
                feh=self.param['feh'])
        # spec_data = {}
        # f_orig = spec_flux["col0"].data
        # w_orig = spec_lambda["col0"].data
        speci = interpolate.interp1d(w_orig, f_orig)
        lamb = np.arange(2500, 10001 + 0.5, 0.5)
        y = speci(lamb)
        # erg/s/cm2/A --> photo/s/m2/A
        all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13
        self.sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX'))

    def getGSObj(self, psf, g1=0, g2=0, flux=None, filt=None, tel=None, exptime=150.):
        if flux == None:
            flux = self.getElectronFluxFilt(filt, tel, exptime)
        # star = galsim.Gaussian(sigma=1.e-8, flux=1.)
        star = galsim.DeltaFunction()
        star = star.withFlux(flux)
        final = galsim.Convolve(psf, star)
        return final
        
    def getGSObj_multiband(self, tel, psf_list, bandpass_list, filt, nphotons_tot=None, g1=0, g2=0, exptime=150.):
        if len(psf_list) != len(bandpass_list):
            raise ValueError("!!!The number of PSF profiles and the number of bandpasses must be equal.")
        objs = []
        if nphotons_tot == None:
            nphotons_tot = self.getElectronFluxFilt(filt, tel, exptime)

        try:
            full = integrate_sed_bandpass(sed=self.sed, bandpass=filt.bandpass_full)
        except Exception as e:
            print(e)
            return -1

        for i in range(len(bandpass_list)):
            bandpass = bandpass_list[i]
            try:
                sub = integrate_sed_bandpass(sed=self.sed, bandpass=bandpass)
            except Exception as e:
                print(e)
                return -1
        
            ratio = sub/full

            if not (ratio == -1 or (ratio != ratio)):
                nphotons = ratio * nphotons_tot
            else:
                return -1
            star = galsim.DeltaFunction()
            star = star.withFlux(nphotons)
            star = galsim.Convolve(psf, star)
            objs.append(star)
        final = galsim.Sum(objs)
        return final