TestSpecDisperse.py

import unittest
from ObservationSim.MockObject.SpecDisperser import rotate90, SpecDisperser

from ObservationSim.Config import ChipOutput, Config
from ObservationSim.Instrument import Telescope, Chip, FilterParam, Filter, FocalPlane
from ObservationSim.MockObject import MockObject, Star
from ObservationSim.PSF import PSFGauss

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

import matplotlib.pyplot as plt

from lmfit.models import LinearModel, GaussianModel

from ObservationSim.Config.Header import generateExtensionHeader
import math
import yaml


def getAngle132(x1=0, y1=0, z1=0, x2=0, y2=0, z2=0, x3=0, y3=0, z3=0):
    cosValue = 0;
    angle = 0;

    x11 = x1 - x3;
    y11 = y1 - y3;
    z11 = z1 - z3;

    x22 = x2 - x3;
    y22 = y2 - y3;
    z22 = z2 - z3;

    tt = np.sqrt((x11 * x11 + y11 * y11 + z11 * z11) * (x22 * x22 + y22 * y22 + z22 * z22));
    if (tt == 0):
        return 0;

    cosValue = (x11 * x22 + y11 * y22 + z11 * z22) / tt;

    if (cosValue > 1):
        cosValue = 1;
    if (cosValue < -1):
        cosValue = -1;
    angle = math.acos(cosValue);
    return angle * 360 / (2 * math.pi);


def fit_SingleGauss(xX, yX, contmX, iHa0):
    background = LinearModel(prefix='line_')
    pars = background.make_params(intercept=yX.max(), slope=0)
    pars = background.guess(yX, x=xX)

    gauss = GaussianModel(prefix='g_')
    pars.update(gauss.make_params())
    pars['g_center'].set(iHa0, min=iHa0 - 3, max=iHa0 + 3)
    pars['g_amplitude'].set(50, min=0)
    pars['g_sigma'].set(12, min=0.0001)

    mod = gauss + background
    init = mod.eval(pars, x=xX)
    outX = mod.fit(yX, pars, x=xX)
    compsX = outX.eval_components(x=xX)
    # print(outX.fit_report(min_correl=0.25))
    # print outX.params['g_center']
    outX.fit_report(min_correl=0.25)
    # print(outX.fit_report(min_correl=0.25))
    line_slopeX = float(outX.fit_report(min_correl=0.25).split('line_slope:')[1].split('+/-')[0]) * contmX
    err_line_slopeX = float(
        outX.fit_report(min_correl=0.25).split('line_slope:')[1].split('+/-')[1].split('(')[0]) * contmX

    line_interceptX = float(outX.fit_report(min_correl=0.25).split('line_intercept:')[1].split('+/-')[0]) * contmX
    err_line_interceptX = float(
        outX.fit_report(min_correl=0.25).split('line_intercept:')[1].split('+/-')[1].split('(')[0]) * contmX

    sigmaX = float(outX.fit_report(min_correl=0.25).split('g_sigma:')[1].split('+/-')[0])
    err_sigmaX = float(outX.fit_report(min_correl=0.25).split('g_sigma:')[1].split('+/-')[1].split('(')[0])

    fwhmX = float(outX.fit_report(min_correl=0.25).split('g_fwhm:')[1].split('+/-')[0])
    err_fwhmX = float(outX.fit_report(min_correl=0.25).split('g_fwhm:')[1].split('+/-')[1].split('(')[0])

    centerX = float(outX.fit_report(min_correl=0.25).split('g_center:')[1].split('+/-')[0])
    err_centerX = float(outX.fit_report(min_correl=0.25).split('g_center:')[1].split('+/-')[1].split('(')[0])

    return sigmaX, err_sigmaX, fwhmX, err_fwhmX, centerX, err_centerX

def produceObj(x,y,chip, ra, dec, pa):
    pos_img = galsim.PositionD(chip.bound.xmin + x, chip.bound.ymin + y)

    param = {}
    param["star"] = 1
    param["id"] = 1
    param["ra"] = chip.img.wcs.posToWorld(pos_img).ra.deg
    param["dec"] = chip.img.wcs.posToWorld(pos_img).dec.deg
    param["z"] = 0
    param["sed_type"] = 1
    param["model_tag"] = 1
    param["mag_use_normal"] = 10

    obj = Star(param)

    header_wcs = generateExtensionHeader(
        xlen=chip.npix_x,
        ylen=chip.npix_y,
        ra=ra,
        dec=dec,
        pa=pa,
        gain=chip.gain,
        readout=chip.read_noise,
        dark=chip.dark_noise,
        saturation=90000,
        psize=chip.pix_scale,
        row_num=chip.rowID,
        col_num=chip.colID,
        extName='raw')

    pos_img, offset, local_wcs, _ = obj.getPosImg_Offset_WCS(img=chip.img, chip=chip, img_header=header_wcs)
    wave = np.arange(2500, 11000.5, 0.5)
    # sedLen = wave.shape[0]
    flux = pow(wave, -2) * 1e8
    flux[200] = flux[200] * 10
    flux[800] = flux[800] * 30
    flux[2000] = flux[2000] * 5

    obj.sed = Table(np.array([wave, flux]).T,
                    names=('WAVELENGTH', 'FLUX'))
    return obj, pos_img


class TestSpecDisperse(unittest.TestCase):

    def __init__(self, methodName='runTest',conff = '', throughputf = ''):
        super(TestSpecDisperse,self).__init__(methodName)
        self.conff = conff
        self.throughputf = throughputf

    def test_rotate901(self):
        m = np.array([[1,2,3,4,5],[6,7,8,9,10],[11,12,13,14,15],[16,17,18,19,20],[21,22,23,24,25]])
        m1 = np.array([[21,16,11,6,1],[22,17,12,7,2],[23,18,13,8,3],[24,19,14,9,4],[25,20,15,10,5]])
        m2 = np.array([[5,10,15,20,25],[4,9,14,19,24],[3,8,13,18,23],[2,7,12,17,22],[1,6,11,16,21]])
        xc = 2
        yc = 2
        isClockwise = 0
        m1, xc1, yc1 = rotate90(array_orig=m, xc=xc, yc=yc, isClockwise=isClockwise)
        self.assertTrue(xc1-xc == 0)
        self.assertTrue(yc1-yc == 0)
        self.assertTrue(np.sum(m-m1) == 0)

    def test_rotate902(self):
        m = np.array([[1,2,3,4,5],[6,7,8,9,10],[11,12,13,14,15],[16,17,18,19,20],[21,22,23,24,25]])
        m1 = np.array([[21,16,11,6,1],[22,17,12,7,2],[23,18,13,8,3],[24,19,14,9,4],[25,20,15,10,5]])
        m2 = np.array([[5,10,15,20,25],[4,9,14,19,24],[3,8,13,18,23],[2,7,12,17,22],[1,6,11,16,21]])
        xc = 2
        yc = 2
        isClockwise =1
        m1, xc1, yc1 = rotate90(array_orig=m, xc=xc, yc=yc, isClockwise=isClockwise)
        self.assertTrue(xc1-xc == 0)
        self.assertTrue(yc1-yc == 0)
        self.assertTrue(np.sum(m-m2) == 0)


    def test_Specdistperse1(self):
        star = galsim.Gaussian(fwhm=0.39)
        g_img = galsim.Image(100, 100, scale=0.074)
        starImg = star.drawImage(image=g_img)

        wave = np.arange(6200, 10000.5, 0.5)
        # sedLen = wave.shape[0]
        flux = pow(wave, -2) * 1e8
        # flux[200] = flux[200] * 10
        # flux[800] = flux[800] * 30
        # flux[2000] = flux[2000] * 5

        sed = Table(np.array([wave, flux]).T,
                    names=('WAVELENGTH', 'FLUX'))
        conff = self.conff
        throughput_f = self.throughputf
        thp = Table.read(throughput_f)
        thp_i = interpolate.interp1d(thp['WAVELENGTH'], thp['SENSITIVITY'])
        sdp = SpecDisperser(orig_img=starImg, xcenter=0, ycenter=0, origin=[100, 100], tar_spec=sed, band_start=6200,
                            band_end=10000, isAlongY=0, conf=conff, gid=0)
        spec = sdp.compute_spec_orders()
        Aimg = spec['A'][0]
        wave_pix = spec['A'][5]
        wave_pos = spec['A'][3]
        sens = spec['A'][6]
        sh = Aimg.shape
        spec_pix = np.zeros(sh[1])
        for i in range(sh[1]):
            spec_pix[i] = sum(Aimg[:, i])
        # figure()
        # imshow(Aimg)

        wave_flux = np.zeros(wave_pix.shape[0])
        for i in np.arange(1, wave_pix.shape[0] - 1):
            w = wave_pix[i]
            thp_w = thp_i(w)
            deltW = (w - wave_pix[i - 1]) / 2 + (wave_pix[i + 1] - w) / 2
            f = spec_pix[wave_pos[0] - 1 + i]

            if 6200 <= w <= 10000:
                f = f / thp_w
            else:
                f = 0
            wave_flux[i] = f / deltW
        sed_i = interpolate.interp1d(sed['WAVELENGTH'], sed['FLUX'])
        ids = wave_pix < 9700
        ids1 = wave_pix[ids] > 6500
        print('Spec disperse flux test')
        self.assertTrue(np.mean((wave_flux[ids][ids1] - sed_i(wave_pix[ids][ids1]))/sed_i(wave_pix[ids][ids1]))<0.004)
        plt.figure()
        plt.plot(wave_pix, wave_flux)
        plt.plot(sed['WAVELENGTH'], sed['FLUX'])
        plt.xlim(6200, 10000)
        plt.ylim(1, 3)
        plt.yscale('log')
        plt.xlabel('$\lambda$')
        plt.ylabel('$F\lambda$')
        plt.legend(['extracted', 'input'])
        plt.show()

    def test_Specdistperse2(self):

        psf_fwhm = 0.39
        pix_scale = 0.074
        star = galsim.Gaussian(fwhm=psf_fwhm)
        g_img = galsim.Image(100, 100, scale=pix_scale)
        starImg = star.drawImage(image=g_img)

        wave = np.arange(6200, 10000.5, 0.5)
        # sedLen = wave.shape[0]
        flux = pow(wave, -2) * 1e8
        flux[200] = flux[200] * 10
        flux[800] = flux[800] * 30
        flux[2000] = flux[2000] * 5

        sed = Table(np.array([wave, flux]).T,
                    names=('WAVELENGTH', 'FLUX'))
        conff = self.conff
        throughput_f = self.throughputf
        thp = Table.read(throughput_f)
        thp_i = interpolate.interp1d(thp['WAVELENGTH'], thp['SENSITIVITY'])
        sdp = SpecDisperser(orig_img=starImg, xcenter=0, ycenter=0, origin=[100, 100], tar_spec=sed, band_start=6200,
                            band_end=10000, isAlongY=0, conf=conff, gid=0)
        spec = sdp.compute_spec_orders()
        Aimg = spec['A'][0]
        wave_pix = spec['A'][5]
        wave_pos = spec['A'][3]
        sens = spec['A'][6]
        sh = Aimg.shape
        spec_pix = np.zeros(sh[1])
        for i in range(sh[1]):
            spec_pix[i] = sum(Aimg[:, i])
        # figure()
        # imshow(Aimg)

        wave_flux = np.zeros(wave_pix.shape[0])
        for i in np.arange(1, wave_pix.shape[0] - 1):
            w = wave_pix[i]
            thp_w = thp_i(w)
            deltW = (w - wave_pix[i - 1]) / 2 + (wave_pix[i + 1] - w) / 2
            f = spec_pix[wave_pos[0] - 1 + i]

            if 6200 <= w <= 10000:
                f = f / thp_w
            else:
                f = 0
            wave_flux[i] = f / deltW
        sed_i = interpolate.interp1d(sed['WAVELENGTH'], sed['FLUX'])
        input_em_lam = 6600
        ids = wave_pix < input_em_lam+200
        ids1 = wave_pix[ids] > input_em_lam-200
        deltLamda_pix = (max(wave_pix[ids][ids1]) - min(wave_pix[ids][ids1])) / (wave_pix[ids][ids1].shape[0] - 1)
        _, _, fwhmx, fwhmx_err, center, center_err = fit_SingleGauss(wave_pix[ids][ids1], wave_flux[ids][ids1], 1.0, 6600)

        print('Emission line position and shape test')

        self.assertTrue(input_em_lam-center < deltLamda_pix)

        self.assertTrue(fwhmx/deltLamda_pix*pix_scale - psf_fwhm  < np.abs(0.01))
        # print('error is ',np.mean((wave_flux[ids][ids1] - sed_i(wave_pix[ids][ids1]))/sed_i(wave_pix[ids][ids1])))
        # self.assertTrue(np.mean((wave_flux[ids][ids1] - sed_i(wave_pix[ids][ids1]))/sed_i(wave_pix[ids][ids1]))<0.004)
        plt.figure()
        plt.plot(wave_pix, wave_flux)
        plt.plot(sed['WAVELENGTH'], sed['FLUX'])
        plt.xlim(6200, 10000)
        plt.ylim(1, 75)
        plt.yscale('log')
        plt.xlabel('$\lambda$')
        plt.ylabel('$F\lambda$')
        plt.legend(['extracted', 'input'])
        plt.show()

    def test_Specdistperse3(self):

        psf_fwhm = 0.39
        pix_scale = 0.074
        star = galsim.Gaussian(fwhm=psf_fwhm)
        g_img = galsim.Image(100, 100, scale=pix_scale)
        starImg = star.drawImage(image=g_img)

        wave = np.arange(6200, 10000.5, 0.5)
        # sedLen = wave.shape[0]
        flux = pow(wave, -2) * 1e8
        flux[200] = flux[200] * 10
        flux[800] = flux[800] * 30
        flux[2000] = flux[2000] * 5

        sed = Table(np.array([wave, flux]).T,
                    names=('WAVELENGTH', 'FLUX'))
        conff = self.conff
        throughput_f = self.throughputf
        thp = Table.read(throughput_f)
        thp_i = interpolate.interp1d(thp['WAVELENGTH'], thp['SENSITIVITY'])
        sdp = SpecDisperser(orig_img=starImg, xcenter=0, ycenter=0, origin=[100, 100], tar_spec=sed, band_start=6200,
                            band_end=8000, isAlongY=0, conf=conff, gid=0)
        sdp1 = SpecDisperser(orig_img=starImg, xcenter=0, ycenter=0, origin=[100, 100], tar_spec=sed, band_start=8000,
                             band_end=10000, isAlongY=0, conf=conff, gid=0)
        spec = sdp.compute_spec_orders()
        spec1 = sdp1.compute_spec_orders()
        Aimg = spec['A'][0] + spec1['A'][0]
        wave_pix = spec['A'][5]
        wave_pos = spec['A'][3]
        sens = spec['A'][6]
        sh = Aimg.shape
        spec_pix = np.zeros(sh[1])
        for i in range(sh[1]):
            spec_pix[i] = sum(Aimg[:, i])


        wave_flux = np.zeros(wave_pix.shape[0])
        for i in np.arange(1, wave_pix.shape[0] - 1):
            w = wave_pix[i]
            thp_w = thp_i(w)
            deltW = (w - wave_pix[i - 1]) / 2 + (wave_pix[i + 1] - w) / 2
            f = spec_pix[wave_pos[0] - 1 + i]

            if 6200 <= w <= 10000:
                f = f / thp_w
            else:
                f = 0
            wave_flux[i] = f / deltW

        sdp2 = SpecDisperser(orig_img=starImg, xcenter=0, ycenter=0, origin=[100, 100], tar_spec=sed, band_start=6200,
                             band_end=10000, isAlongY=0, conf=conff, gid=0)

        spec2 = sdp2.compute_spec_orders()
        Aimg2 = spec2['A'][0]

        spec_pix2 = np.zeros(sh[1])
        for i in range(sh[1]):
            spec_pix2[i] = sum(Aimg2[:, i])

        wave_flux2 = np.zeros(wave_pix.shape[0])
        for i in np.arange(1, wave_pix.shape[0] - 1):
            w = wave_pix[i]
            thp_w = thp_i(w)
            deltW = (w - wave_pix[i - 1]) / 2 + (wave_pix[i + 1] - w) / 2
            f = spec_pix2[wave_pos[0] - 1 + i]

            if 6200 <= w <= 10000:
                f = f / thp_w
            else:
                f = 0
            wave_flux2[i] = f / deltW

        r1_i = interpolate.interp1d(wave_pix, wave_flux2)
        r2_i = interpolate.interp1d(wave_pix, wave_flux)

        print('Spec Splicing test')
        self.assertTrue(r1_i(8000)-r2_i(8000) < np.abs(0.0001))

        # self.assertTrue(fwhmx/deltLamda_pix*pix_scale - psf_fwhm  < np.abs(0.01))
        # print('error is ',np.mean((wave_flux[ids][ids1] - sed_i(wave_pix[ids][ids1]))/sed_i(wave_pix[ids][ids1])))
        # self.assertTrue(np.mean((wave_flux[ids][ids1] - sed_i(wave_pix[ids][ids1]))/sed_i(wave_pix[ids][ids1]))<0.004)
        plt.figure()
        plt.plot(wave_pix, wave_flux2)
        plt.plot(wave_pix, wave_flux)
        # plt.plot(sed['WAVELENGTH'], sed['FLUX'])
        plt.xlim(6200, 10000)
        plt.ylim(1, 4)
        plt.yscale('log')
        plt.xlabel('$\lambda$')
        plt.ylabel('$F\lambda$')
        plt.legend(['one spec', 'split in 8000 A'])
        plt.show()



    def test_double_disperse(self):
        work_dir = "/public/home/fangyuedong/CSST_unittest/CSST/test/"
        # data_dir = "/Volumes/Extreme SSD/SimData/"
        data_dir = "/data/simudata/CSSOSDataProductsSims/data/"
        configFn = '/Users/zhangxin/Work/SlitlessSim/CSST_SIM/CSST_C6/csst-simulation/config/config_C6.yaml'
        normFilterFn = "/Users/zhangxin/Work/SlitlessSim/CSST_SIM/CSST_C6/csst-simulation/Catalog/data/SLOAN_SDSS.g.fits"
        norm_star = Table.read(normFilterFn)
        with open(configFn, "r") as stream:
            try:
                config = yaml.safe_load(stream)
                for key, value in config.items():
                    print(key + " : " + str(value))
            except yaml.YAMLError as exc:
                print(exc)
        # config = Config.read_config(configFn)
        # path_dict = Config.config_dir(config,work_dir=work_dir, data_dir=data_dir)


        filter_param = FilterParam()
        focal_plane = FocalPlane(survey_type=config["obs_setting"]["survey_type"])
        chip = Chip(1, config=config)
        filter_id, filter_type = chip.getChipFilter()
        filt = Filter(filter_id=filter_id, filter_type=filter_type, filter_param=filter_param,
                      ccd_bandpass=chip.effCurve)
        tel = Telescope()

        psf_model = PSFGauss(chip=chip)


        wcs_fp = focal_plane.getTanWCS(float(config["obs_setting"]["ra_center"]), float(config["obs_setting"]["dec_center"]), float(config["obs_setting"]["image_rot"]) * galsim.degrees, chip.pix_scale)
        chip.img = galsim.ImageF(chip.npix_x, chip.npix_y)
        chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin)
        chip.img.wcs = wcs_fp

        obj, pos_img = produceObj(2000,4500, chip,float(config["obs_setting"]["ra_center"]), float(config["obs_setting"]["dec_center"]), float(config["obs_setting"]["image_rot"]))
        obj.drawObj_slitless(
            tel=tel,
            pos_img=pos_img,
            psf_model=psf_model,
            bandpass_list=filt.bandpass_sub_list,
            filt=filt,
            chip=chip,
            g1=0,
            g2=0,
            exptime=150,
            normFilter=norm_star)

        obj, pos_img = produceObj(3685, 6500, chip,float(config["obs_setting"]["ra_center"]), float(config["obs_setting"]["dec_center"]), float(config["obs_setting"]["image_rot"]))
        obj.drawObj_slitless(
            tel=tel,
            pos_img=pos_img,
            psf_model=psf_model,
            bandpass_list=filt.bandpass_sub_list,
            filt=filt,
            chip=chip,
            g1=0,
            g2=0,
            exptime=150,
            normFilter=norm_star)

        obj, pos_img = produceObj(5000, 2500, chip, float(config["obs_setting"]["ra_center"]), float(config["obs_setting"]["dec_center"]), float(config["obs_setting"]["image_rot"]))
        obj.drawObj_slitless(
            tel=tel,
            pos_img=pos_img,
            psf_model=psf_model,
            bandpass_list=filt.bandpass_sub_list,
            filt=filt,
            chip=chip,
            g1=0,
            g2=0,
            exptime=150,
            normFilter=norm_star)

        print('Spec double disperse test')
        from astropy.io import fits
        fits.writeto('test.fits',chip.img.array, overwrite = True)

        # plt.figure()
        # plt.imshow(chip.img.array)
        # plt.show()

    def test_SLSImage_rotation(self):
        from astropy.wcs import WCS
        configFn = '/Users/zhangxin/Work/SlitlessSim/CSST_SIM/CSST_C6/csst-simulation/config/config_C6.yaml'

        with open(configFn, "r") as stream:
            try:
                config = yaml.safe_load(stream)
                for key, value in config.items():
                    print(key + " : " + str(value))
            except yaml.YAMLError as exc:
                print(exc)
        chip = Chip(1, config=config)

        ra=float(config["obs_setting"]["ra_center"])
        dec=float(config["obs_setting"]["dec_center"])
        pa=float(config["obs_setting"]["image_rot"])

        header_wcs1 = generateExtensionHeader(
            xlen=chip.npix_x,
            ylen=chip.npix_y,
            ra=ra,
            dec=dec,
            pa=pa,
            gain=chip.gain,
            readout=chip.read_noise,
            dark=chip.dark_noise,
            saturation=90000,
            psize=chip.pix_scale,
            row_num=chip.rowID,
            col_num=chip.colID,
            extName='raw', center_rot=0)

        center = np.array([chip.npix_x / 2, chip.npix_y / 2])
        h_wcs1 = WCS(header_wcs1)
        x1, y1 = center + [100,0]
        sky_1 = h_wcs1.pixel_to_world(x1,y1)

        rot_angle = 1
        header_wcs2 = generateExtensionHeader(
            xlen=chip.npix_x,
            ylen=chip.npix_y,
            ra=ra,
            dec=dec,
            pa=pa,
            gain=chip.gain,
            readout=chip.read_noise,
            dark=chip.dark_noise,
            saturation=90000,
            psize=chip.pix_scale,
            row_num=chip.rowID,
            col_num=chip.colID,
            extName='raw', center_rot=rot_angle)

        h_wcs2 = WCS(header_wcs2)
        x2, y2 = h_wcs2.world_to_pixel(sky_1)
        angle = getAngle132(x1,y1,0,x2,y2,0,center[0],center[1],0)
        print(angle)
        self.assertTrue(rot_angle - angle < np.abs(0.001))

        rot_angle = 10
        header_wcs2 = generateExtensionHeader(
            xlen=chip.npix_x,
            ylen=chip.npix_y,
            ra=ra,
            dec=dec,
            pa=pa,
            gain=chip.gain,
            readout=chip.read_noise,
            dark=chip.dark_noise,
            saturation=90000,
            psize=chip.pix_scale,
            row_num=chip.rowID,
            col_num=chip.colID,
            extName='raw', center_rot=rot_angle)

        h_wcs2 = WCS(header_wcs2)
        x2, y2 = h_wcs2.world_to_pixel(sky_1)
        angle = getAngle132(x1, y1, 0, x2, y2, 0, center[0], center[1], 0)
        print(angle)
        self.assertTrue(rot_angle - angle < np.abs(0.001))

        rot_angle = 50
        header_wcs2 = generateExtensionHeader(
            xlen=chip.npix_x,
            ylen=chip.npix_y,
            ra=ra,
            dec=dec,
            pa=pa,
            gain=chip.gain,
            readout=chip.read_noise,
            dark=chip.dark_noise,
            saturation=90000,
            psize=chip.pix_scale,
            row_num=chip.rowID,
            col_num=chip.colID,
            extName='raw', center_rot=rot_angle)

        h_wcs2 = WCS(header_wcs2)
        x2, y2 = h_wcs2.world_to_pixel(sky_1)
        angle = getAngle132(x1, y1, 0, x2, y2, 0, center[0], center[1], 0)
        print(angle)
        self.assertTrue(rot_angle - angle < np.abs(0.001))


        chip = Chip(27, config=config)

        ra = float(config["obs_setting"]["ra_center"])
        dec = float(config["obs_setting"]["dec_center"])
        pa = float(config["obs_setting"]["image_rot"])

        header_wcs1 = generateExtensionHeader(
            xlen=chip.npix_x,
            ylen=chip.npix_y,
            ra=ra,
            dec=dec,
            pa=pa,
            gain=chip.gain,
            readout=chip.read_noise,
            dark=chip.dark_noise,
            saturation=90000,
            psize=chip.pix_scale,
            row_num=chip.rowID,
            col_num=chip.colID,
            extName='raw', center_rot=0)

        center = np.array([chip.npix_x / 2, chip.npix_y / 2])
        h_wcs1 = WCS(header_wcs1)
        x1, y1 = center + [100, 0]
        sky_1 = h_wcs1.pixel_to_world(x1, y1)

        rot_angle = 1
        header_wcs2 = generateExtensionHeader(
            xlen=chip.npix_x,
            ylen=chip.npix_y,
            ra=ra,
            dec=dec,
            pa=pa,
            gain=chip.gain,
            readout=chip.read_noise,
            dark=chip.dark_noise,
            saturation=90000,
            psize=chip.pix_scale,
            row_num=chip.rowID,
            col_num=chip.colID,
            extName='raw', center_rot=rot_angle)

        h_wcs2 = WCS(header_wcs2)
        x2, y2 = h_wcs2.world_to_pixel(sky_1)
        angle = getAngle132(x1, y1, 0, x2, y2, 0, center[0], center[1], 0)
        print(angle)
        self.assertTrue(rot_angle - angle < np.abs(0.001))

        rot_angle = 10
        header_wcs2 = generateExtensionHeader(
            xlen=chip.npix_x,
            ylen=chip.npix_y,
            ra=ra,
            dec=dec,
            pa=pa,
            gain=chip.gain,
            readout=chip.read_noise,
            dark=chip.dark_noise,
            saturation=90000,
            psize=chip.pix_scale,
            row_num=chip.rowID,
            col_num=chip.colID,
            extName='raw', center_rot=rot_angle)

        h_wcs2 = WCS(header_wcs2)
        x2, y2 = h_wcs2.world_to_pixel(sky_1)
        angle = getAngle132(x1, y1, 0, x2, y2, 0, center[0], center[1], 0)
        print(angle)
        self.assertTrue(rot_angle - angle < np.abs(0.001))

        rot_angle = 50
        header_wcs2 = generateExtensionHeader(
            xlen=chip.npix_x,
            ylen=chip.npix_y,
            ra=ra,
            dec=dec,
            pa=pa,
            gain=chip.gain,
            readout=chip.read_noise,
            dark=chip.dark_noise,
            saturation=90000,
            psize=chip.pix_scale,
            row_num=chip.rowID,
            col_num=chip.colID,
            extName='raw', center_rot=rot_angle)

        h_wcs2 = WCS(header_wcs2)
        x2, y2 = h_wcs2.world_to_pixel(sky_1)
        angle = getAngle132(x1, y1, 0, x2, y2, 0, center[0], center[1], 0)
        print(angle)
        self.assertTrue(rot_angle - angle < np.abs(0.001))




if __name__ == '__main__':
    conff='/Users/zhangxin/Work/SlitlessSim/CSST_SIM/CSST_C6/csst-simulation/ObservationSim/Instrument/data/sls_conf/CSST_GI2.conf'
    throughputf='/Users/zhangxin/Work/SlitlessSim/CSST_SIM/CSST_C6/csst-simulation/ObservationSim/Instrument/data/sls_conf/GI.Throughput.1st.fits'

    suit = unittest.TestSuite()
    case1 = TestSpecDisperse('test_Specdistperse1',conff,throughputf)
    suit.addTest(case1)
    case2 = TestSpecDisperse('test_Specdistperse2', conff, throughputf)
    suit.addTest(case2)
    case3 = TestSpecDisperse('test_Specdistperse3', conff, throughputf)
    suit.addTest(case3)
    case4 = TestSpecDisperse('test_double_disperse', conff, throughputf)
    suit.addTest(case4)
    case5 = TestSpecDisperse('test_SLSImage_rotation')
    suit.addTest(case5)

    unittest.TextTestRunner(verbosity=2).run(suit)
    # runner = unittest.TextTestRunner()
    # runner.run(suit)