diff --git a/ObservationSim/Config/Header/ImageHeader.py b/ObservationSim/Config/Header/ImageHeader.py index 0281dcddae182eafd299637289f9cf2c3b0ceead..86a46cbc4f634ef8d7d1000e66665c3965f8d5bd 100644 --- a/ObservationSim/Config/Header/ImageHeader.py +++ b/ObservationSim/Config/Header/ImageHeader.py @@ -7,7 +7,7 @@ from astropy.io import fits import astropy.wcs as pywcs from collections import OrderedDict -from scipy import math +# from scipy import math import random import os diff --git a/ObservationSim/MockObject/MockObject.py b/ObservationSim/MockObject/MockObject.py index 5b6dffb9d3582753fe5e36576aa671de533dd428..a8316ab88d612af2396a185127c681c119914303 100755 --- a/ObservationSim/MockObject/MockObject.py +++ b/ObservationSim/MockObject/MockObject.py @@ -256,9 +256,28 @@ class MockObject(object): spec_orders = sdp.compute_spec_orders() for k, v in spec_orders.items(): img_s = v[0] + ######################################################### + # DEBUG + ######################################################### + # print("before convolveGaussXorders, img_s:", img_s) + nan_ids = np.isnan(img_s) + if img_s[nan_ids].shape[0] > 0: + # img_s[nan_ids] = 0 + print("DEBUG: before convolveGaussXorders specImg nan num is", img_s[nan_ids].shape[0]) + ######################################################### img_s, orig_off = convolveGaussXorders(img_s, xOrderSigPlus[k]) origin_order_x = v[1] - orig_off origin_order_y = v[2] - orig_off + + ######################################################### + # DEBUG + ######################################################### + # print("DEBUG: orig_off is", orig_off) + nan_ids = np.isnan(img_s) + if img_s[nan_ids].shape[0] > 0: + img_s[nan_ids] = 0 + print("DEBUG: specImg nan num is", img_s[nan_ids].shape[0]) + ######################################################### specImg = galsim.ImageF(img_s) photons = galsim.PhotonArray.makeFromImage(specImg) photons.x += origin_order_x diff --git a/ObservationSim/MockObject/_util.py b/ObservationSim/MockObject/_util.py index ffb978b910cc63d1bf281148554af29c54292d51..c6cb573d053f59bb2927f9af7c37f3f34c2b6570 100755 --- a/ObservationSim/MockObject/_util.py +++ b/ObservationSim/MockObject/_util.py @@ -545,10 +545,9 @@ def tag_sed(h5file, model_tag, teff=5000, logg=2, feh=0): flux = np.array(h5file["sed"][path][()]).ravel() return path, wave, flux -from astropy.modeling.models import Gaussian2D -from scipy import signal def convolveGaussXorders(img=None, sigma = 1): - + from astropy.modeling.models import Gaussian2D + from scipy import signal offset = int(np.ceil(sigma*10)) g_size = 2*offset+1 diff --git a/ObservationSim/ObservationSim.py b/ObservationSim/ObservationSim.py index 5e79ec38bb51437f32599d6fce48a53d1129bc7a..4426177ddc9334c4a486d5ee63aff7b98713bd44 100755 --- a/ObservationSim/ObservationSim.py +++ b/ObservationSim/ObservationSim.py @@ -7,6 +7,8 @@ import psutil from astropy.io import fits from datetime import datetime +import traceback + from ObservationSim.Config import config_dir, ChipOutput from ObservationSim.Config.Header import generatePrimaryHeader, generateExtensionHeader from ObservationSim.Instrument import Telescope, Filter, FilterParam, FocalPlane, Chip @@ -77,7 +79,7 @@ class Observation(object): chip_output.logger.info(':::::::::::::::::::::::::::END:::::::::::::::::::::::::::::::::::') if self.config["psf_setting"]["psf_model"] == "Gauss": - psf_model = PSFGauss(chip=chip) + psf_model = PSFGauss(chip=chip, psfRa=self.config["psf_setting"]["psf_rcont"]) elif self.config["psf_setting"]["psf_model"] == "Interp": psf_model = PSFInterp(chip=chip, PSF_data_file=self.path_dict["psf_dir"]) else: @@ -207,7 +209,8 @@ class Observation(object): ) except Exception as e: - print(e) + # print(e) + traceback.print_exc() chip_output.logger.error(e) continue @@ -325,7 +328,8 @@ class Observation(object): # print("object omitted", flush=True) continue except Exception as e: - print(e) + # print(e) + traceback.print_exc() chip_output.logger.error(e) pass # Unload SED: diff --git a/ObservationSim/PSF/PSFGauss.py b/ObservationSim/PSF/PSFGauss.py index 508388055bec2e652418f5d41e69685b2d9c0191..7d09c3c67096e318d3935189a4e28b91c6b9deee 100644 --- a/ObservationSim/PSF/PSFGauss.py +++ b/ObservationSim/PSF/PSFGauss.py @@ -6,12 +6,16 @@ from scipy.interpolate import interp1d from ObservationSim.PSF.PSFModel import PSFModel class PSFGauss(PSFModel): - def __init__(self, chip, fwhm=0.187, sigSpin=0., psfRa=0.15): + def __init__(self, chip, fwhm=0.187, sigSpin=0., psfRa=None): self.pix_size = chip.pix_scale self.chip = chip - self.fwhm = fwhm + if psfRa is None: + self.fwhm = fwhm + self.sigGauss = 0.15 + else: + self.fwhm = self.fwhmGauss(r=psfRa) + self.sigGauss = psfRa # 80% light radius self.sigSpin = sigSpin - self.sigGauss = psfRa # 80% light radius self.psf = galsim.Gaussian(flux=1.0,fwhm=fwhm) def perfGauss(self, r, sig): diff --git a/ObservationSim/PSF/PSFInterp.py b/ObservationSim/PSF/PSFInterp.py index 058acca032f174f0dcaf36d75f7f146937e7e29d..45bf211c8d648f9675cb6ab68c3b2f92ec6cef43 100644 --- a/ObservationSim/PSF/PSFInterp.py +++ b/ObservationSim/PSF/PSFInterp.py @@ -350,12 +350,21 @@ class PSFInterp(PSFModel): if galsimGSObject: imPSFt = np.zeros([257,257]) imPSFt[0:256, 0:256] = imPSF + # imPSFt[120:130, 0:256] = 1. img = galsim.ImageF(imPSFt, scale=pixSize) gsp = galsim.GSParams(folding_threshold=folding_threshold) - self.psf = galsim.InterpolatedImage(img, gsparams=gsp).rotate(pointing_pa*galsim.degrees) - - return self.PSFspin(x=px/0.01, y=py/0.01) + ############TEST: START + # Use sheared PSF to test the PSF orientation + # self.psf = galsim.InterpolatedImage(img, gsparams=gsp).shear(g1=0.8, g2=0.) + ############TEST: END + self.psf = galsim.InterpolatedImage(img, gsparams=gsp) + wcs = chip.img.wcs.local(pos_img) + scale = galsim.PixelScale(0.074) + self.psf = wcs.toWorld(scale.toImage(self.psf), image_pos=(pos_img)) + + # return self.PSFspin(x=px/0.01, y=py/0.01) + return self.psf, galsim.Shear(e=0., beta=(np.pi/2)*galsim.radians) return imPSF def PSFspin(self, x, y): diff --git a/config/config_NGP_dev.yaml b/config/config_NGP_dev.yaml index 27bcee0a2c8e72d34d5251badcfee2aede2694b0..9f5ffdbaedf48c0ba4eba4e25458b3d7aa8d7643 100644 --- a/config/config_NGP_dev.yaml +++ b/config/config_NGP_dev.yaml @@ -10,9 +10,9 @@ # Can add some of the command-line arguments here as well; # OK to pass either way or both, as long as they are consistent # work_dir: "/public/home/fangyuedong/sim_code_release/CSST/test/" -work_dir: "/share/home/fangyuedong/csst-simulation/workplace/" +work_dir: "/share/home/fangyuedong/test_psf_rot/csst-simulation/workplace/" data_dir: "/share/simudata/CSSOSDataProductsSims/data/" -run_name: "header_TEST" +run_name: "PSF_TEST" # (Optional) a file of point list # if you just want to run default pointing: @@ -34,7 +34,7 @@ run_option: out_cat_only: NO # Only simulate stars? - star_only: NO + star_only: YES # Only simulate galaxies? galaxy_only: NO @@ -48,7 +48,7 @@ obs_setting: # "Photometric": simulate photometric chips only # "Spectroscopic": simulate slitless spectroscopic chips only # "All": simulate full focal plane - survey_type: "Photometric" + survey_type: "All" # Exposure time [seconds] exp_time: 150. diff --git a/run_NGP.pbs b/run_NGP.pbs index 2763d98eaa45dda2f11a481db4561acd226d1c7f..f4d181ea980833987d570ea233891327598c42a3 100755 --- a/run_NGP.pbs +++ b/run_NGP.pbs @@ -1,7 +1,7 @@ #!/bin/bash #PBS -N SIMS -#PBS -l nodes=wcl-1:ppn=32+wcl-6:ppn=32 +#PBS -l nodes=wcl-1:ppn=8+wcl-2:ppn=8+wcl-3:ppn=8+wcl-4:ppn=8+wcl-5:ppn=8+wcl-6:ppn=8 #PBS -u fangyuedong ###PBS -j oe @@ -12,5 +12,5 @@ date echo $NP # mpirun -np $NP python /public/home/fangyuedong/test/CSST/run_sim.py config_NGP.yaml -c /public/home/fangyuedong/test/CSST/config -mpirun -np $NP python3 /share/home/fangyuedong/csst-simulation/run_sim.py config_NGP_dev.yaml -c /share/home/fangyuedong/csst-simulation/config +mpirun -np $NP python3 /share/home/fangyuedong/test_psf_rot/csst-simulation/run_sim.py config_NGP_dev.yaml -c /share/home/fangyuedong/test_psf_rot/csst-simulation/config