SkybackgroundMap.py

from ObservationSim.MockObject.SpecDisperser import SpecDisperser
from ObservationSim.MockObject.SpecDisperser import rotate90

import galsim
import numpy as np
from astropy.table import Table
from scipy import interpolate

import galsim

import os

try:
    import importlib.resources as pkg_resources
except ImportError:
    # Try backported to PY<37 'importlib_resources'
    import importlib_resources as pkg_resources


###calculate sky map by sky SED

def calculateSkyMap_split_g(skyMap=None, blueLimit=4200, redLimit=6500, skyfn='sky_emiss_hubble_50_50_A.dat', conf=[''], pixelSize=0.074, isAlongY=0,
                            split_pos=3685):
    # skyMap = np.ones([yLen, xLen], dtype='float32')
    #
    # if isAlongY == 1:
    #     skyMap = np.ones([xLen, yLen], dtype='float32')

    # for i in range(len(conf)):
    #     conf[i] = os.path.join(SLSSIM_PATH, conf[i])
    conf1 = conf[0]
    conf2 = conf[0]
    if np.size(conf) == 2:
        conf2 = conf[1]

    skyImg = galsim.Image(skyMap, xmin=0, ymin=0)

    tbstart = blueLimit
    tbend = redLimit

    fimg = np.zeros_like(skyMap)

    fImg = galsim.Image(fimg)
    # skyfn = os.path.join(SLSSIM_PATH, skyfn)
    with pkg_resources.path('ObservationSim.MockObject.data', skyfn) as data_path:
        skySpec = np.loadtxt(data_path)
    # skySpec = np.loadtxt(skyfn)
    spec = Table(np.array([skySpec[:, 0], skySpec[:, 1]]).T, names=('WAVELENGTH', 'FLUX'))

    if isAlongY == 0:
        directParm = 0
    if isAlongY ==1:
        directParm = 1

    if split_pos >= skyImg.array.shape[directParm]:
        skyImg1 = galsim.Image(skyImg.array)
        origin1 = [0, 0]
        # sdp = specDisperser.specDisperser(orig_img=skyImg1, xcenter=skyImg1.center.x, ycenter=skyImg1.center.y,
        #                                   full_img=fimg, tar_spec=spec, band_start=tbstart, band_end=tbend,
        #                                   origin=origin1,
        #                                   conf=conf1)
        # sdp.compute_spec_orders()

        
        sdp = SpecDisperser(orig_img=skyImg1, xcenter=skyImg1.center.x, ycenter=skyImg1.center.y, origin=origin1,
                        tar_spec=spec,
                        band_start=tbstart, band_end=tbend,
                        conf=conf2)

        spec_orders = sdp.compute_spec_orders()

        for k, v in spec_orders.items():
            img_s = v[0]
            origin_order_x = v[1]
            origin_order_y = v[2]
            ssImg = galsim.ImageF(img_s)
            ssImg.setOrigin(origin_order_x, origin_order_y)
            bounds = ssImg.bounds & fImg.bounds
            if bounds.area() == 0:
                continue
            fImg[bounds] = fImg[bounds] + ssImg[bounds]
        

        
    else:

        skyImg1 = galsim.Image(skyImg.array[:, 0:split_pos])
        origin1 = [0, 0]
        skyImg2 = galsim.Image(skyImg.array[:, split_pos:-1])
        origin2 = [0, split_pos]

        # sdp = specDisperser.specDisperser(orig_img=skyImg1, xcenter=skyImg1.center.x, ycenter=skyImg1.center.y,
        #                                   full_img=fimg, tar_spec=spec, band_start=tbstart, band_end=tbend,
        #                                   origin=origin1,
        #                                   conf=conf1)

        # sdp.compute_spec_orders()

        sdp = SpecDisperser(orig_img=skyImg1, xcenter=skyImg1.center.x, ycenter=skyImg1.center.y, origin=origin1,
                        tar_spec=spec,
                        band_start=tbstart, band_end=tbend,
                        conf=conf1)

        spec_orders = sdp.compute_spec_orders()

        for k, v in spec_orders.items():
            img_s = v[0]
            origin_order_x = v[1]
            origin_order_y = v[2]
            ssImg = galsim.ImageF(img_s)
            ssImg.setOrigin(origin_order_x, origin_order_y)
            bounds = ssImg.bounds & fImg.bounds
            if bounds.area() == 0:
                continue
            fImg[bounds] = fImg[bounds] + ssImg[bounds]

        
        sdp = SpecDisperser(orig_img=skyImg2, xcenter=skyImg2.center.x, ycenter=skyImg2.center.y, origin=origin2,
                        tar_spec=spec,
                        band_start=tbstart, band_end=tbend,
                        conf=conf2)

        spec_orders = sdp.compute_spec_orders()

        for k, v in spec_orders.items():
            img_s = v[0]
            origin_order_x = v[1]
            origin_order_y = v[2]
            ssImg = galsim.ImageF(img_s)
            ssImg.setOrigin(origin_order_x, origin_order_y)
            bounds = ssImg.bounds & fImg.bounds
            fImg[bounds] = fImg[bounds] + ssImg[bounds]

    if isAlongY == 1:
        fimg, tmx, tmy = rotate90(array_orig=fImg.array, xc=0, yc=0, isClockwise=0)
    else:
        fimg = fImg.array

    fimg = fimg * pixelSize * pixelSize

    return fimg
    
def calculateSkyMap(xLen=9232, yLen=9126, blueLimit=4200, redLimit=6500,
                    skyfn='sky_emiss_hubble_50_50_A.dat', conf='', pixelSize=0.074, isAlongY=0):
    skyMap = np.ones([yLen, xLen], dtype='float32')

    if isAlongY == 1:
        skyMap = np.ones([xLen, yLen], dtype='float32')

    skyImg = galsim.Image(skyMap)

    tbstart = blueLimit
    tbend = redLimit

    fimg = np.zeros_like(skyMap)
    fImg = galsim.Image(fimg)
    with pkg_resources.path('ObservationSim.MockObject.data', skyfn) as data_path:
        skySpec = np.loadtxt(data_path)
    # skySpec = np.loadtxt(skyfn)
    
    spec = Table(np.array([skySpec[:, 0], skySpec[:, 1]]).T, names=('WAVELENGTH', 'FLUX'))

    sdp = SpecDisperser(orig_img=skyImg, xcenter=skyImg.center.x, ycenter=skyImg.center.y, origin=[1, 1],
                        tar_spec=spec,
                        band_start=tbstart, band_end=tbend,
                        conf=conf)

    spec_orders = sdp.compute_spec_orders()

    for k, v in spec_orders.items():
        img_s = v[0]
        origin_order_x = v[1]
        origin_order_y = v[2]
        ssImg = galsim.ImageF(img_s)
        ssImg.setOrigin(origin_order_x, origin_order_y)
        bounds = ssImg.bounds & fImg.bounds
        fImg[bounds] = fImg[bounds] + ssImg[bounds]

    if isAlongY == 1:
        fimg, tmx, tmy = rotate90(array_orig=fImg.array, xc=0, yc=0, isClockwise=0)
    else:
        fimg = fImg.array
        
    fimg = fimg * pixelSize * pixelSize

    return fimg