From a29e7df11bf7eacc6d8d8b38f5b24832bef1848b Mon Sep 17 00:00:00 2001 From: Zhang Xin Date: Sat, 7 Oct 2023 15:44:13 +0800 Subject: [PATCH] modify unittest name spec & straylight --- tests/test_SpecDisperse.py | 688 +++++++++++++++++++++++++++++++++++++ tests/test_Straylight.py | 190 ++++++++++ 2 files changed, 878 insertions(+) create mode 100644 tests/test_SpecDisperse.py create mode 100644 tests/test_Straylight.py diff --git a/tests/test_SpecDisperse.py b/tests/test_SpecDisperse.py new file mode 100644 index 0000000..5581dd3 --- /dev/null +++ b/tests/test_SpecDisperse.py @@ -0,0 +1,688 @@ +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) \ No newline at end of file diff --git a/tests/test_Straylight.py b/tests/test_Straylight.py new file mode 100644 index 0000000..479cbbc --- /dev/null +++ b/tests/test_Straylight.py @@ -0,0 +1,190 @@ +import unittest +from ObservationSim.Straylight import Straylight + +import numpy as np +import math +import astropy.constants as cons +import galsim +from astropy.table import Table +from scipy import interpolate + +import matplotlib.pyplot as plt + +hubbleAverZodiacal = {'nuv':0.0035,'u':0.0163,'g':0.1109,'r':0.1471,'i':0.1568,'z':0.0953,'y':0.0283} +hubbleAverEarthShine = {'nuv':0.00024,'u':0.0051,'g':0.0506,'r':0.0591,'i':0.0568,'z':0.0315,'y':0.0090} + +def transRaDec2D(ra, dec): + x1 = np.cos(dec / 57.2957795) * np.cos(ra / 57.2957795); + y1 = np.cos(dec / 57.2957795) * np.sin(ra / 57.2957795); + z1 = np.sin(dec / 57.2957795); + return np.array([x1, y1, z1]) + + +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 calculateAnglePwithEarth(sat = np.array([0,0,0]), pointing = np.array([0,0,0]), sun = np.array([0,0,0])): + modSat = np.sqrt(sat[0]*sat[0] + sat[1]*sat[1]+sat[2]*sat[2]) + modPoint = np.sqrt(pointing[0]*pointing[0] + pointing[1]*pointing[1] + pointing[2]*pointing[2]) + withLocalZenithAngle = (pointing[0] * sat[0] + pointing[1] * sat[1] + pointing[2] * sat[2]) / (modPoint*modSat) + + innerM_sat_sun = sat[0] * sun[0] + sat[1] * sun[1] + sat[2] * sun[2] + cosAngle = innerM_sat_sun / (modSat * cons.au.value/1000) + isInSunSide = 1 + if (cosAngle < -0.3385737): #cos109.79 + isInSunSide = -1; + elif cosAngle >= -0.3385737 and cosAngle <= 0.3385737: + isInSunSide = 0; + + return math.acos(withLocalZenithAngle)*180/math.pi,isInSunSide + +class TestStraylight(unittest.TestCase): + + def __init__(self, methodName='runTest',datFn = '', filter = 'i', grating = "GI"): + super(TestStraylight,self).__init__(methodName) + self.pointingData = np.loadtxt(datFn, dtype=np.double) + self.filter = filter + self.grating = grating + + + def test_EarthShineFilter(self): + d_sh = self.pointingData.shape + sl_e_pix = np.zeros([d_sh[0],3],dtype=np.double) + + for i in np.arange(d_sh[0]): + # if i > 50: + # continue + ju = self.pointingData[i, 5] + # pointing = transRaDec2D(self.pointingData[i, 0], self.pointingData[i, 1]) + # print(ju, pointing, surveylist[i,3:9]) + sl = Straylight(jtime=ju, sat_pos=self.pointingData[i, 6:9], pointing_radec=np.array([self.pointingData[i, 0], self.pointingData[i, 1]]),sun_pos=self.pointingData[i,9:12]) + e1, py = sl.calculateEarthShineFilter(filter=self.filter) + earthZenithAngle, isInSunSide = calculateAnglePwithEarth(sat=self.pointingData[i, 6:9], pointing= sl.pointing, sun=self.pointingData[i,9:12]) + # e2, _ = sl.calculateZodiacalFilter2(filter='i', sun_pos=sl.sun_pos) + # e3 = sl.calculateStarLightFilter(filter='i', pointYaxis=py) + # e_all = sl.calculateStrayLightFilter(filter='i') + # s_pix, spec = sl.calculateStrayLightGrating(grating='GI') + sl_e_pix[i,0] = e1 + sl_e_pix[i, 1] = earthZenithAngle + sl_e_pix[i, 2] = isInSunSide + median = np.median(sl_e_pix[:,0]) + print(' average Earthshine %s: %e' % (self.filter, median)) + self.assertTrue(median-hubbleAverEarthShine[self.filter] < 0.1) + plt.figure() + ids1 = sl_e_pix[:, 2] == 1 + ids2 = sl_e_pix[:, 2] != 1 + plt.plot(sl_e_pix[ids1, 0], sl_e_pix[ids1, 1], 'r.') + plt.plot(sl_e_pix[ids2, 0], sl_e_pix[ids2, 1], 'b.') + plt.legend(['In Sun Side', 'In Earths shadow']) + plt.xlabel('straylight-earthshine(e-/pixel/s)') + plt.ylabel('Angle with local zenith(degree)') + plt.show() + + def test_ZodiacalFilter(self): + d_sh = self.pointingData.shape + sl_e_pix = np.zeros([d_sh[0],2],dtype=np.double) + + for i in np.arange(d_sh[0]): + ju = self.pointingData[i, 5] + sl = Straylight(jtime=ju, sat_pos=self.pointingData[i, 6:9], pointing_radec=np.array([self.pointingData[i, 0], self.pointingData[i, 1]]),sun_pos=self.pointingData[i,9:12]) + e1, _ = sl.calculateZodiacalFilter2(filter=self.filter, sun_pos=sl.sun_pos) + sl_e_pix[i,0] = e1 + sl_e_pix[i,1] = getAngle132(x1=self.pointingData[i,9], y1=self.pointingData[i,10], z1=self.pointingData[i,11], x2=sl.pointing[0], + y2=sl.pointing[1], z2=sl.pointing[2], x3=0, y3=0, z3=0) + plt.figure() + plt.plot(sl_e_pix[:, 0], sl_e_pix[:, 1], 'r.') + plt.xlabel('straylight-zodiacal(e-/pixel/s)') + plt.ylabel('Angle between pointing and sun(degree)') + plt.show() + median = np.median(sl_e_pix[:,0]) + print(' average Zodiacal %s: %f' % (self.filter, median)) + self.assertTrue(median-hubbleAverZodiacal[self.filter] < 0.1) + + def test_StarFilter(self): + d_sh = self.pointingData.shape + sl_e_pix = np.zeros(d_sh[0],dtype=np.double) + + tnum = 10 + for i in np.arange(tnum): + # if i > 50: + # continue + ju = self.pointingData[i, 5] + # pointing = transRaDec2D(self.pointingData[i, 0], self.pointingData[i, 1]) + # print(ju, pointing, surveylist[i,3:9]) + sl = Straylight(jtime=ju, sat_pos=self.pointingData[i, 6:9], pointing_radec=np.array([self.pointingData[i, 0], self.pointingData[i, 1]]),sun_pos=self.pointingData[i,9:12]) + e1, py = sl.calculateEarthShineFilter(filter=self.filter) + # e2, _ = sl.calculateZodiacalFilter2(filter='i', sun_pos=sl.sun_pos) + e3 = sl.calculateStarLightFilter(filter=self.filter, pointYaxis=py) + # e_all = sl.calculateStrayLightFilter(filter='i') + # s_pix, spec = sl.calculateStrayLightGrating(grating='GI') + sl_e_pix[i] = e3 + median = np.median(sl_e_pix[0:tnum]) + print(' average Earthshine %s: %e' % (self.filter, median)) + self.assertTrue(median-hubbleAverEarthShine[self.filter] < 0.2) + + def test_GratingStraylight(self): + d_sh = self.pointingData.shape + sl_e_pix = np.zeros(d_sh[0],dtype=np.double) + + tnum = 10 + for i in np.arange(tnum): + # if i > 50: + # continue + ju = self.pointingData[i, 5] + # pointing = transRaDec2D(self.pointingData[i, 0], self.pointingData[i, 1]) + # print(ju, pointing, surveylist[i,3:9]) + sl = Straylight(jtime=ju, sat_pos=self.pointingData[i, 6:9], pointing_radec=np.array([self.pointingData[i, 0], self.pointingData[i, 1]]),sun_pos=self.pointingData[i,9:12]) + # e1, py = sl.calculateEarthShineFilter(filter=self.filter) + # e2, _ = sl.calculateZodiacalFilter2(filter='i', sun_pos=sl.sun_pos) + # e3 = sl.calculateStarLightFilter(filter=self.filter, pointYaxis=py) + # e_all = sl.calculateStrayLightFilter(filter='i') + s_pix, spec = sl.calculateStrayLightGrating(grating=self.grating) + sl_e_pix[i] = s_pix + plt.figure() + plt.plot(spec['WAVELENGTH'], spec['FLUX'], 'r') + plt.xlabel('WAVELENGTH') + plt.ylabel('F$\lambda$(erg/s/cm2/A/arcsec2)') + plt.xlim(2000,10000) + plt.show() + median = np.median(sl_e_pix[0:tnum]) + print(' average Earthshine %s: %e' % (self.grating, median)) + self.assertTrue(median < 0.8) + + + + + +if __name__ == '__main__': + + suit = unittest.TestSuite() + # case1 = TestStraylight('test_EarthShineFilter',datFn = 'Straylight_test.dat', filter = 'i') + # suit.addTest(case1) + # case2 = TestStraylight('test_ZodiacalFilter',datFn = 'Straylight_test.dat',filter = 'i') + # suit.addTest(case2) + # case3 = TestStraylight('test_StarFilter', datFn='Straylight_test.dat', filter='i') + # suit.addTest(case3) + case4 = TestStraylight('test_GratingStraylight', datFn='Straylight_test.dat', grating = 'GI') + suit.addTest(case4) + unittest.TextTestRunner(verbosity=2).run(suit) \ No newline at end of file -- GitLab