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 astropy.constants as cons

import os

import time

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, flat_cube = None, zoldial_spec = None):
    # 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)
    try:
        with pkg_resources.files('ObservationSim.MockObject.data').joinpath(skyfn) as data_path:
            skySpec = np.loadtxt(data_path)
    except AttributeError:
        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 zoldial_spec is not None:
        deltL = 0.5
        lamb = np.arange(2000, 11000, deltL)

        speci = interpolate.interp1d(zoldial_spec['WAVELENGTH'], zoldial_spec['FLUX'])

        y = speci(lamb)
        # erg/s/cm2/A --> photo/s/m2/A
        s_flux = y * lamb / (cons.h.value * cons.c.value) * 1e-13
        spec = Table(np.array([lamb, s_flux]).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()

        y_len = skyMap.shape[0]
        x_len = skyMap.shape[1]
        delt_x = 100
        delt_y = 100

        sub_y_start_arr = np.arange(0, y_len, delt_y)
        sub_y_end_arr = sub_y_start_arr + delt_y
        sub_y_end_arr[-1] = min(sub_y_end_arr[-1], y_len)

        sub_x_start_arr = np.arange(0, x_len, delt_x)
        sub_x_end_arr = sub_x_start_arr + delt_x
        sub_x_end_arr[-1] = min(sub_x_end_arr[-1], x_len)

        for i,k1 in enumerate(sub_y_start_arr):
            sub_y_s = k1
            sub_y_e = sub_y_end_arr[i]

            sub_y_center = (sub_y_s+sub_y_e)/2.

            for j,k2 in enumerate(sub_x_start_arr):
                sub_x_s = k2
                sub_x_e = sub_x_end_arr[j]

                skyImg_sub = galsim.Image(skyImg.array[sub_y_s:sub_y_e, sub_x_s:sub_x_e])
                origin_sub = [sub_y_s, sub_x_s]
                sub_x_center = (sub_x_s + sub_x_e) / 2.

                sdp = SpecDisperser(orig_img=skyImg_sub, xcenter=sub_x_center, ycenter=sub_y_center, origin=origin_sub,
                                tar_spec=spec,
                                band_start=tbstart, band_end=tbend,
                                conf=conf2,
                                flat_cube=flat_cube, ignoreBeam=['D','E'])

                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=skyImg1, xcenter=skyImg1.center.x, ycenter=skyImg1.center.y, origin=origin1,
        #                 tar_spec=spec,
        #                 band_start=tbstart, band_end=tbend,
        #                 conf=conf2,
        #                 flat_cube=flat_cube, ignoreBeam=['D','E'])
        #
        # 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:])
        # 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()
        y_len = skyMap.shape[0]
        x_len = skyMap.shape[1]
        delt_x = 500
        delt_y = y_len

        sub_y_start_arr = np.arange(0, y_len, delt_y)
        sub_y_end_arr = sub_y_start_arr + delt_y
        sub_y_end_arr[-1] = min(sub_y_end_arr[-1], y_len)
        
        delt_x = split_pos-0
        sub_x_start_arr = np.arange(0, split_pos, delt_x)
        sub_x_end_arr = sub_x_start_arr + delt_x
        sub_x_end_arr[-1] = min(sub_x_end_arr[-1], split_pos)

        for i,k1 in enumerate(sub_y_start_arr):
            sub_y_s = k1
            sub_y_e = sub_y_end_arr[i]

            sub_y_center = (sub_y_s+sub_y_e)/2.

            for j,k2 in enumerate(sub_x_start_arr):
                sub_x_s = k2
                sub_x_e = sub_x_end_arr[j]
                # print(i,j,sub_y_s, sub_y_e,sub_x_s,sub_x_e)
                T1 = time.time()
                skyImg_sub = galsim.Image(skyImg.array[sub_y_s:sub_y_e, sub_x_s:sub_x_e])
                origin_sub = [sub_y_s, sub_x_s]
                sub_x_center = (sub_x_s + sub_x_e) / 2.

                sdp = SpecDisperser(orig_img=skyImg_sub, xcenter=sub_x_center, ycenter=sub_y_center, origin=origin_sub,
                                tar_spec=spec,
                                band_start=tbstart, band_end=tbend,
                                conf=conf1,
                                flat_cube=flat_cube)

                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]
        
                T2 = time.time()

                print('time: %s ms'% ((T2 - T1)*1000))

        delt_x = x_len-split_pos
        sub_x_start_arr = np.arange(split_pos, x_len, delt_x)
        sub_x_end_arr = sub_x_start_arr + delt_x
        sub_x_end_arr[-1] = min(sub_x_end_arr[-1], x_len)

        for i, k1 in enumerate(sub_y_start_arr):
            sub_y_s = k1
            sub_y_e = sub_y_end_arr[i]

            sub_y_center = (sub_y_s + sub_y_e) / 2.

            for j, k2 in enumerate(sub_x_start_arr):
                sub_x_s = k2
                sub_x_e = sub_x_end_arr[j]
                # print(i,j,sub_y_s, sub_y_e,sub_x_s,sub_x_e)
                
                T1 = time.time()

                skyImg_sub = galsim.Image(skyImg.array[sub_y_s:sub_y_e, sub_x_s:sub_x_e])
                origin_sub = [sub_y_s, sub_x_s]
                sub_x_center = (sub_x_s + sub_x_e) / 2.

                sdp = SpecDisperser(orig_img=skyImg_sub, xcenter=sub_x_center, ycenter=sub_y_center, origin=origin_sub,
                                tar_spec=spec,
                                band_start=tbstart, band_end=tbend,
                                conf=conf2,
                                flat_cube=flat_cube)

                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]
                T2 = time.time()

                print('time: %s ms'% ((T2 - T1)*1000))

    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)
    try:
        with pkg_resources.files('ObservationSim.MockObject.data').joinpath(skyfn) as data_path:
            skySpec = np.loadtxt(data_path)
    except AttributeError:
        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