add something for unittest

ObservationSim/Config/Header/ImageHeader.py

@@ -164,7 +164,7 @@ def calcaluteSLSRotSkyCoor(pix_xy = None,rot_angle = 1, xlen = 9216, ylen = 9232
 
 
 ##9232 9216  898 534 1309 60 -40  -23.4333
-def WCS_def(xlen = 9216, ylen = 9232, gapy = 898.0, gapx1 = 534, gapx2 = 1309, ra = 60, dec = -40, pa = -23.433,psize = 0.074, row_num = 1, col_num = 1, filter = 'GI'):
+def WCS_def(xlen = 9216, ylen = 9232, gapy = 898.0, gapx1 = 534, gapx2 = 1309, ra = 60, dec = -40, pa = -23.433,psize = 0.074, row_num = 1, col_num = 1, filter = 'GI', center_rot = 1):
 
     """ Creat a wcs frame for CCST with multiple extensions
 
@@ -251,7 +251,7 @@ def WCS_def(xlen = 9216, ylen = 9232, gapy = 898.0, gapx1 = 534, gapx2 = 1309, r
 
             # test_center_o = w.wcs_pix2world(np.array([[xlen / 2, ylen / 2]]), 1)
 
-            sls_rot = 1
+            sls_rot = center_rot
             if i > 2:
                 sls_rot = -sls_rot
 
@@ -401,7 +401,7 @@ def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec
 
     return h_prim
 
-def generateExtensionHeader(xlen = 9216, ylen = 9232,ra = 60, dec = -40, pa = -23.433, gain = 1.0, readout = 5.0, dark = 0.02, saturation=90000, psize = 0.074, row_num = 1, col_num = 1, extName='SCI'):
+def generateExtensionHeader(xlen = 9216, ylen = 9232,ra = 60, dec = -40, pa = -23.433, gain = 1.0, readout = 5.0, dark = 0.02, saturation=90000, psize = 0.074, row_num = 1, col_num = 1, extName='SCI', center_rot = 1):
 
     e_header_fn = os.path.split(os.path.realpath(__file__))[0] + '/extension_header.header'
     f = open(os.path.split(os.path.realpath(__file__))[0] + '/filter.lst')
@@ -464,7 +464,7 @@ def generateExtensionHeader(xlen = 9216, ylen = 9232,ra = 60, dec = -40, pa = -2
 
     # h_ext['POS_ANG'] = pa
     header_wcs = WCS_def(xlen=xlen, ylen=ylen, gapy=898.0, gapx1=534, gapx2=1309, ra=ra, dec=dec, pa=pa, psize=psize,
-                         row_num=row_num, col_num=col_num, filter = h_ext['FILTER'])
+                         row_num=row_num, col_num=col_num, filter = h_ext['FILTER'], center_rot = center_rot)
 
     h_ext['CRPIX1'] = header_wcs['CRPIX1']
     h_ext['CRPIX2'] = header_wcs['CRPIX2']

ObservationSim/MockObject/SpecDisperser/SpecDisperser.py

@@ -248,7 +248,7 @@ class SpecDisperser(object):
         if self.isAlongY == 1:
             model, _, _ = rotate90(array_orig=model, isClockwise=0)
 
-        return model, originOut_x, originOut_y
+        return model, originOut_x, originOut_y, dxpix, dypix, lam_beam,ysens
 
     def writerSensitivityFile(self, conffile='', beam='', w=None, sens=None):
         orders = {'A': '1st', 'B': '0st', 'C': '2st', 'D': '-1st', 'E': '-2st'}

setup.py

@@ -51,6 +51,7 @@ setup(name='CSSTSim',
         'ObservationSim.Instrument.data.sls_conf': ['*.conf', '*.fits'],
         'ObservationSim.Instrument.data.flatCube': ['*.fits'],
         'Catalog.data': ['*.fits'],
+        'ObservationSim.Config.Header':['*.header','*.lst'],
     },
     ext_modules = cythonize(extensions),
 )

tests/TestSpecDisperse.py

 import unittest
 from ObservationSim.MockObject.SpecDisperser import rotate90, SpecDisperser
 
-from ObservationSim.Config import ConfigDir, ReadConfig, ChipOutput
+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
@@ -15,6 +15,36 @@ 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_')
@@ -25,7 +55,7 @@ def fit_SingleGauss(xX, yX, contmX, iHa0):
     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.01)
+    pars['g_sigma'].set(12, min=0.0001)
 
     mod = gauss + background
     init = mod.eval(pars, x=xX)
@@ -54,7 +84,7 @@ def fit_SingleGauss(xX, yX, contmX, iHa0):
 
     return sigmaX, err_sigmaX, fwhmX, err_fwhmX, centerX, err_centerX
 
-def produceObj(x,y,chip):
+def produceObj(x,y,chip, ra, dec, pa):
     pos_img = galsim.PositionD(chip.bound.xmin + x, chip.bound.ymin + y)
 
     param = {}
@@ -69,7 +99,22 @@ def produceObj(x,y,chip):
 
     obj = Star(param)
 
-    pos_img, offset, local_wcs = obj.getPosImg_Offset_WCS(img=chip.img, chip=chip)
+    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
@@ -346,28 +391,37 @@ class TestSpecDisperse(unittest.TestCase):
         work_dir = "/public/home/fangyuedong/CSST_unittest/CSST/test/"
         # data_dir = "/Volumes/Extreme SSD/SimData/"
         data_dir = "/data/simudata/CSSOSDataProductsSims/data/"
-        path_dict = ConfigDir(work_dir=work_dir, data_dir=data_dir)
-        config = ReadConfig(path_dict["config_file"])
-
-        filter_param = FilterParam(filter_dir=path_dict["filter_dir"])
-        focal_plane = FocalPlane(survey_type=config["survey_type"])
-        chip = Chip(1, ccdEffCurve_dir=path_dict["ccd_dir"], CRdata_dir=path_dict["CRdata_dir"],
-                    normalize_dir=path_dict["normalize_dir"], sls_dir=path_dict["sls_dir"],
-                    config=config)
+        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(optEffCurve_path=path_dict["mirror_file"])
+        tel = Telescope()
 
         psf_model = PSFGauss(chip=chip)
 
 
-        wcs_fp = focal_plane.getTanWCS(config["ra_center"], config["dec_center"], config["image_rot"], chip.pix_scale)
+        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)
+        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,
@@ -378,9 +432,9 @@ class TestSpecDisperse(unittest.TestCase):
             g1=0,
             g2=0,
             exptime=150,
-            normFilter=chip.normF_star)
+            normFilter=norm_star)
 
-        obj, pos_img = produceObj(3685, 6500, chip)
+        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,
@@ -391,9 +445,9 @@ class TestSpecDisperse(unittest.TestCase):
             g1=0,
             g2=0,
             exptime=150,
-            normFilter=chip.normF_star)
+            normFilter=norm_star)
 
-        obj, pos_img = produceObj(5000, 2500, chip)
+        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,
@@ -404,7 +458,7 @@ class TestSpecDisperse(unittest.TestCase):
             g1=0,
             g2=0,
             exptime=150,
-            normFilter=chip.normF_star)
+            normFilter=norm_star)
 
         print('Spec double disperse test')
         from astropy.io import fits
@@ -414,11 +468,208 @@ class TestSpecDisperse(unittest.TestCase):
         # 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='/data/simudata/CSSOSDataProductsSims/data/CONF/CSST_GI2.conf'
-    throughputf='/data/simudata/CSSOSDataProductsSims/data/CONF/GI.Throughput.1st.fits'
+    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)
@@ -429,6 +680,8 @@ if __name__ == '__main__':
     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()