add_sky_background.py

import numpy as np
import galsim
from ObservationSim.Straylight import calculateSkyMap_split_g
from ObservationSim.Instrument import FilterParam

from astropy.time import Time
from datetime import datetime, timezone

def add_sky_background_sci(self, chip, filt, tel, pointing, catalog, obs_param):
    
    # Get exposure time
    if (obs_param) and ("exptime" in obs_param) and (obs_param["exptime"] is not None):
        exptime = obs_param["exptime"]
    else:
        exptime = pointing.exp_time
    
    flat_normal = np.ones_like(chip.img.array)
    if obs_param["flat_fielding"] == True:
        flat_normal = flat_normal * chip.flat_img.array / np.mean(chip.flat_img.array)
    if obs_param["shutter_effect"] == True:
        flat_normal = flat_normal * chip.shutter_img
        flat_normal = np.array(flat_normal, dtype='float32')
        self.updateHeaderInfo(header_flag='ext', keys = ['SHTSTAT'], values = [True])
    else:
        self.updateHeaderInfo(header_flag='ext', keys = ['SHTSTAT','SHTOPEN1','SHTCLOS0'], values = [True,self.h_ext['SHTCLOS1'],self.h_ext['SHTOPEN0']])
    
    if obs_param["enable_straylight_model"]:
        # Filter.sky_background, Filter.zodical_spec will be updated
        filt.setFilterStrayLightPixel(
            jtime = pointing.jdt,
            sat_pos = np.array([pointing.sat_x, pointing.sat_y, pointing.sat_z]),
            pointing_radec = np.array([pointing.ra,pointing.dec]),
            sun_pos = np.array([pointing.sun_x, pointing.sun_y, pointing.sun_z]))
        self.chip_output.Log_info("================================================")
        self.chip_output.Log_info("sky background + stray light pixel flux value: %.5f"%(filt.sky_background))
    
    if chip.survey_type == "photometric":
        sky_map = filt.getSkyNoise(exptime = exptime)
        sky_map = sky_map * np.ones_like(chip.img.array) * flat_normal
        sky_map = galsim.Image(array=sky_map)
    else:
        # chip.loadSLSFLATCUBE(flat_fn='flat_cube.fits')
        sky_map = calculateSkyMap_split_g(
                skyMap=flat_normal,
                blueLimit=filt.blue_limit,
                redLimit=filt.red_limit,
                conf=chip.sls_conf,
                pixelSize=chip.pix_scale,
                isAlongY=0,
                flat_cube=chip.flat_cube, 
                zoldial_spec = filt.zodical_spec)
        sky_map = (sky_map + filt.sky_background)*exptime
    
    # sky_map = sky_map * tel.pupil_area * obs_param["exptime"]
    chip.img += sky_map
    return chip, filt, tel, pointing

def add_sky_flat_calibration(self, chip, filt, tel, pointing, catalog, obs_param):

    if not hasattr(self, 'h_ext'):
        _, _ = self.prepare_headers(chip=chip, pointing=pointing)
    chip_wcs = galsim.FitsWCS(header = self.h_ext)

    # Get exposure time
    if (obs_param) and ("exptime" in obs_param) and (obs_param["exptime"] is not None):
        exptime = obs_param["exptime"]
    else:
        exptime = pointing.exp_time
    
    skyback_level = obs_param["flat_level"]

    filter_param = FilterParam()
    sky_level_filt = obs_param["flat_level_filt"]
    norm_scaler = skyback_level/exptime / filter_param.param[sky_level_filt][5]

    flat_normal = np.ones_like(chip.img.array)
    if obs_param["flat_fielding"] == True:
        flat_normal = flat_normal * chip.flat_img.array / np.mean(chip.flat_img.array)
    if obs_param["shutter_effect"] == True:
        flat_normal = flat_normal * chip.shutter_img
        flat_normal = np.array(flat_normal, dtype='float32')
        if self.overall_config["output_setting"]["shutter_output"] == True:    # output 16-bit shutter effect image with pixel value <=65535
            shutt_gsimg = galsim.ImageUS(chip.shutter_img*6E4)
            shutt_gsimg.write("%s/ShutterEffect_%s_1.fits" % (self.chip_output.subdir, str(chip.chipID).rjust(2, '0')))
        self.updateHeaderInfo(header_flag='ext', keys = ['SHTSTAT'], values = [True])
    else:
        self.updateHeaderInfo(header_flag='ext', keys = ['SHTSTAT','SHTOPEN1','SHTCLOS0'], values = [True,self.h_ext['SHTCLOS1'],self.h_ext['SHTOPEN0']])
    

    if chip.survey_type == "photometric":
        sky_map = flat_normal * np.ones_like(chip.img.array) * norm_scaler * filter_param.param[chip.filter_type][5] / tel.pupil_area * exptime
    elif chip.survey_type == "spectroscopic":
        # flat_normal = np.ones_like(chip.img.array)
        if obs_param["flat_fielding"] == True:
            
            flat_normal = flat_normal * chip.flat_img.array / np.mean(chip.flat_img.array)
        if obs_param["shutter_effect"] == True:
            
            flat_normal = flat_normal * chip.shutter_img
            flat_normal = np.array(flat_normal, dtype='float32')
        sky_map = calculateSkyMap_split_g(
            skyMap=flat_normal,
            blueLimit=filt.blue_limit,
            redLimit=filt.red_limit,
            conf=chip.sls_conf,
            pixelSize=chip.pix_scale,
            isAlongY=0,
            flat_cube=chip.flat_cube)
        sky_map = sky_map * norm_scaler * exptime
    
    chip.img += sky_map
    return chip, filt, tel, pointing

def add_sky_background(self, chip, filt, tel, pointing, catalog, obs_param):
    if not hasattr(self, 'h_ext'):
        _, _ = self.prepare_headers(chip=chip, pointing=pointing)
    chip_wcs = galsim.FitsWCS(header = self.h_ext)

    if "flat_level" not in obs_param or "flat_level_filt" not in obs_param:
        chip, filt, tel, pointing = self.add_sky_background_sci(chip, filt, tel, pointing, catalog, obs_param)
    else:
        if obs_param.get('flat_level') is None or obs_param.get('flat_level_filt')is None:
            chip, filt, tel, pointing = self.add_sky_background_sci(chip, filt, tel, pointing, catalog, obs_param)
        else:
            chip, filt, tel, pointing = self.add_sky_flat_calibration(chip, filt, tel, pointing, catalog, obs_param)

    # renew header info
    datetime_obs = datetime.utcfromtimestamp(pointing.timestamp)
    datetime_obs = datetime_obs.replace(tzinfo=timezone.utc)
    t_obs = Time(datetime_obs)
    
    ##ccd刷新2s,等待0.5s，开始曝光
    t_obs_renew = Time(t_obs.mjd - (2.+0.5) / 86400., format="mjd")

    t_obs_utc = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(t_obs_renew.unix).replace(tzinfo=timezone.utc).timestamp(), 1))
    self.updateHeaderInfo(header_flag='prim', keys = ['DATE-OBS'], values = [t_obs_utc.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5]])

    #dark time : 曝光时间+刷新后等带时间0.5s+关闭快门时间1.5s+管快门后读出前等待0.5s
    self.updateHeaderInfo(header_flag='ext', keys = ['DARKTIME'], values = [0.+0.0+0.0+pointing.exp_time])


    return chip, filt, tel, pointing