diff --git a/ObservationSim/Config/Header/ImageHeader.py b/ObservationSim/Config/Header/ImageHeader.py index d5b5399afc0da15b355e429ed32dcfa776084ccc..4c608a6a3e394ffa8d673ab2949e22ffa134a751 100644 --- a/ObservationSim/Config/Header/ImageHeader.py +++ b/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'] diff --git a/ObservationSim/MockObject/SpecDisperser/SpecDisperser.py b/ObservationSim/MockObject/SpecDisperser/SpecDisperser.py index 926c2d2ae0f1f6b30fce92d1c5bcc7710df31d54..07d099a60bde65c74043ab03c34f0fbbf1bad585 100644 --- a/ObservationSim/MockObject/SpecDisperser/SpecDisperser.py +++ b/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'} diff --git a/setup.py b/setup.py index 70464d0ae22275e81a579b97753a2eada60d4f8f..b6c3f599ef14b58c2f8e8a4323c7ec16b3219bdb 100644 --- a/setup.py +++ b/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), ) diff --git a/tests/TestSpecDisperse.py b/tests/TestSpecDisperse.py index 8e0f9cbe39168d3fae09cdd1d7bba446bb38a463..5581dd3e62cc63daebb4408df82513167fb1e473 100644 --- a/tests/TestSpecDisperse.py +++ b/tests/TestSpecDisperse.py @@ -1,7 +1,7 @@ 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()