From 64eb7956671f54d4fa6de10f6047fbea7e1bf6af Mon Sep 17 00:00:00 2001 From: Fang Yuedong Date: Thu, 1 Aug 2024 01:58:30 +0000 Subject: [PATCH] Revert "Merge branch 'milky_way_extinction' into 'release_v3.0'" This reverts merge request !28 --- catalog/C9_Catalog.py | 44 +-- config/config_overall.yaml | 10 +- config/obs_config_SCI.yaml | 2 +- observation_sim/ObservationSim.py | 45 +-- observation_sim/config/ChipOutput.py | 4 +- observation_sim/config/header/ImageHeader.py | 2 +- observation_sim/instruments/chip/effects.py | 2 +- .../instruments/data/filters/fgs_sub.list | 2 +- .../instruments/data/filters/g_sub.list | 4 +- .../instruments/data/filters/i_sub.list | 4 +- .../instruments/data/filters/nuv_sub.list | 4 +- .../instruments/data/filters/r_sub.list | 4 +- .../instruments/data/filters/u_sub.list | 4 +- .../instruments/data/filters/y_sub.list | 4 +- .../instruments/data/filters/z_sub.list | 4 +- observation_sim/mock_objects/ExtinctionMW.py | 96 ----- observation_sim/mock_objects/Galaxy.py | 127 ++----- observation_sim/mock_objects/MockObject.py | 345 ++++++------------ .../SpecDisperser/SpecDisperser.py | 66 +--- .../SpecDisperser/disperse_c/disperse.pyx | 49 +-- observation_sim/mock_objects/__init__.py | 1 - observation_sim/psf/PSFGauss.py | 4 +- observation_sim/psf/PSFInterpSLS.py | 224 +----------- observation_sim/sim_steps/__init__.py | 12 +- observation_sim/sim_steps/add_objects.py | 27 +- .../sim_steps/add_pattern_noise.py | 4 +- observation_sim/sim_steps/readout_output.py | 8 +- submit_jobs.slurm | 6 +- tools/get_pointing.py | 22 +- tools/get_pointing_accuracy.py | 165 --------- tools/target_location_check.py | 9 +- 31 files changed, 236 insertions(+), 1068 deletions(-) delete mode 100644 observation_sim/mock_objects/ExtinctionMW.py delete mode 100644 tools/get_pointing_accuracy.py diff --git a/catalog/C9_Catalog.py b/catalog/C9_Catalog.py index 146205d..948fc96 100644 --- a/catalog/C9_Catalog.py +++ b/catalog/C9_Catalog.py @@ -12,7 +12,7 @@ from astropy.table import Table from scipy import interpolate from datetime import datetime -from observation_sim.mock_objects import CatalogBase, Star, Galaxy, Quasar, ExtinctionMW +from observation_sim.mock_objects import CatalogBase, Star, Galaxy, Quasar from observation_sim.mock_objects._util import tag_sed, getObservedSED, getABMAG, integrate_sed_bandpass, comoving_dist from observation_sim.astrometry.Astrometry_util import on_orbit_obs_position @@ -110,16 +110,6 @@ class Catalog(CatalogBase): self.max_size = 0. - # [TODO] Milky Way extinction - if "enable_mw_ext_gal" in config["catalog_options"] and config["catalog_options"]["enable_mw_ext_gal"]: - if "planck_ebv_map" not in config["catalog_options"]: - raise ValueError( - "Planck dust map must be given to enable Milky Way extinction calculation for galaxies.") - self.mw_ext = ExtinctionMW() - self.mw_ext.init_ext_model(model_name="odonnell") - self.mw_ext.load_Planck_ext( - file_path=config["catalog_options"]["planck_ebv_map"]) - if "star_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["star_cat"] and not config["catalog_options"]["galaxy_only"]: # Get the cloest star catalog file star_file_name = get_star_cat( @@ -258,12 +248,6 @@ class Catalog(CatalogBase): input_time_str=time_str ) - # [TODO] get Milky Way extinction AVs - if "enable_mw_ext_gal" in self.config["catalog_options"] and self.config["catalog_options"]["enable_mw_ext_gal"]: - MW_Av_arr = self.mw_ext.Av_from_Planck(ra=ra_arr, dec=dec_arr) - else: - MW_Av_arr = np.zeros(len(ra_arr)) - for igals in range(ngals): # # (TEST) # if igals > 100: @@ -275,9 +259,6 @@ class Catalog(CatalogBase): param['ra_orig'] = gals['ra'][igals] param['dec_orig'] = gals['dec'][igals] - # [TODO] get Milky Way extinction AVs - param['mw_Av'] = MW_Av_arr[igals] - if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200): continue @@ -549,30 +530,11 @@ class Catalog(CatalogBase): speci = interpolate.interp1d(wave, flux) lamb = np.arange(2000, 11001+0.5, 0.5) y = speci(lamb) - - # [TODO] Apply Milky Way extinction - if obj.type != 'star' and ("enable_mw_ext_gal" in self.config["catalog_options"] and self.config["catalog_options"]["enable_mw_ext_gal"]): - y = self.mw_ext.apply_extinction(y, Av=obj.mw_Av) - # erg/s/cm2/A --> photon/s/m2/A all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13 sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX')) - # if obj.type == 'quasar': - # # integrate to get the magnitudes - # sed_photon = np.array([sed['WAVELENGTH'], sed['FLUX']]).T - # sed_photon = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array( - # sed_photon[:, 1]), interpolant='nearest') - # sed_photon = galsim.SED( - # sed_photon, wave_type='A', flux_type='1', fast=False) - # interFlux = integrate_sed_bandpass( - # sed=sed_photon, bandpass=self.filt.bandpass_full) - # obj.param['mag_use_normal'] = getABMAG( - # interFlux, self.filt.bandpass_full) - # # mag = getABMAG(interFlux, self.filt.bandpass_full) - # # print("mag diff = %.3f"%(mag - obj.param['mag_use_normal'])) - - # integrate to get the magnitudes - if obj.type == 'quasar' or obj.type == 'galaxy': + if obj.type == 'quasar': + # integrate to get the magnitudes sed_photon = np.array([sed['WAVELENGTH'], sed['FLUX']]).T sed_photon = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array( sed_photon[:, 1]), interpolant='nearest') diff --git a/config/config_overall.yaml b/config/config_overall.yaml index d22dc68..8624cc0 100644 --- a/config/config_overall.yaml +++ b/config/config_overall.yaml @@ -11,7 +11,7 @@ # 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/project/workplace/" -run_name: "ext_on" +run_name: "no_nonlinearity_test" # Project cycle and run counter are used to name the outputs project_cycle: 9 @@ -49,10 +49,6 @@ catalog_options: # rotate galaxy ellipticity rotateEll: 0. # [degree] - # Whether to apply milky way extinction to galaxies - enable_mw_ext_gal: YES - planck_ebv_map: "/public/home/fangyuedong/project/ext_maps/planck/HFI_CompMap_ThermalDustModel_2048_R1.20.fits" - ############################################### # Observation setting ############################################### @@ -69,10 +65,10 @@ obs_setting: # - give a list of indexes of pointings: [ip_1, ip_2...] # - run all pointings: null # Note: only valid when a pointing list is specified - run_pointings: [0, 1, 2, 3, 4] + run_pointings: [1] # Whether to enable astrometric modeling - enable_astrometric_model: YES + enable_astrometric_model: True # Cut by saturation magnitude in which band? cut_in_band: "z" diff --git a/config/obs_config_SCI.yaml b/config/obs_config_SCI.yaml index ef82bcd..405950f 100644 --- a/config/obs_config_SCI.yaml +++ b/config/obs_config_SCI.yaml @@ -16,7 +16,7 @@ obs_id: "00000001" # this setting will only be used if pointing list file is not # Define list of chips # run_chips: [6,7,8,9,11,12,13,14,15,16,17,18,19,20,22,23,24,25] # Photometric chips #run_chips: [1,2,3,4,5,10,21,26,27,28,29,30] # Spectroscopic chips -run_chips: [17, 22] +run_chips: [1] # Define observation sequence call_sequence: diff --git a/observation_sim/ObservationSim.py b/observation_sim/ObservationSim.py index cb060ab..051d6d1 100755 --- a/observation_sim/ObservationSim.py +++ b/observation_sim/ObservationSim.py @@ -23,7 +23,7 @@ class Observation(object): self.filter_param = FilterParam() self.Catalog = Catalog - def prepare_chip_for_exposure(self, chip, ra_cen, dec_cen, pointing, wcs_fp=None, slsPSFOptim=False): + def prepare_chip_for_exposure(self, chip, ra_cen, dec_cen, pointing, wcs_fp=None): # Get WCS for the focal plane if wcs_fp == None: wcs_fp = self.focal_plane.getTanWCS( @@ -34,27 +34,6 @@ class Observation(object): chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin) chip.img.wcs = wcs_fp - chip.slsPSFOptim = slsPSFOptim - if chip.chipID in [1, 2, 3, 4, 5, 10, 21, 26, 27, 28, 29, 30] and slsPSFOptim: - chip.img_stack = {} - for id1 in np.arange(2): - gn = chip_utils.getChipSLSGratingID(chip.chipID)[id1] - orders = {} - # for id2 in ['-2','-1','0','1','2']: - for id2 in ['0', '1']: - o_n = "order"+id2 - allbands = {} - for id3 in ['1', '2', '3', '4']: - w_n = "w"+id3 - allbands[w_n] = galsim.ImageF(chip.npix_x, chip.npix_y) - allbands[w_n].setOrigin( - chip.bound.xmin, chip.bound.ymin) - allbands[w_n].wcs = wcs_fp - orders[o_n] = allbands - chip.img_stack[gn] = orders - else: - chip.img_stack = {} - # Get random generators for this chip chip.rng_poisson, chip.poisson_noise = chip_utils.get_poisson( seed=int(self.config["random_seeds"]["seed_poisson"]) + pointing.id*30 + chip.chipID, sky_level=0.) @@ -66,7 +45,7 @@ class Observation(object): chip.shutter_img = np.ones_like(chip.img.array) else: chip.shutter_img = effects.ShutterEffectArr( - chip.img, t_exp=pointing.exp_time, t_shutter=1.3, dist_bearing=735, dt=1E-3) + chip.img, t_shutter=1.3, dist_bearing=735, dt=1E-3) chip.prnu_img = effects.PRNU_Img(xsize=chip.npix_x, ysize=chip.npix_y, sigma=0.01, seed=int(self.config["random_seeds"]["seed_prnu"]+chip.chipID)) @@ -112,24 +91,17 @@ class Observation(object): input_date_str=date_str, input_time_str=time_str ) - ra_offset, dec_offset = pointing.ra - \ - ra_cen[0], pointing.dec - dec_cen[0] + ra_cen, dec_cen = ra_cen[0], dec_cen[0] else: - ra_offset, dec_offset = 0., 0. - ra_cen = pointing.ra - dec_cen = pointing.dec + ra_cen = pointing.ra + dec_cen = pointing.dec - slsPSFOpt = False # Prepare necessary chip properties for simulation - chip = self.prepare_chip_for_exposure( - chip, ra_cen, dec_cen, pointing, slsPSFOptim=slsPSFOpt) + chip = self.prepare_chip_for_exposure(chip, ra_cen, dec_cen, pointing) # Initialize SimSteps sim_steps = SimSteps(overall_config=self.config, - chip_output=chip_output, - all_filters=self.all_filters, - ra_offset=ra_offset, - dec_offset=dec_offset) + chip_output=chip_output, all_filters=self.all_filters) for step in pointing.obs_param["call_sequence"]: if self.config["run_option"]["out_cat_only"]: @@ -158,9 +130,6 @@ class Observation(object): chip_output.Log_info("check running:1: pointing-%d chip-%d pid-%d memory-%6.2fGB" % (pointing.id, chip.chipID, os.getpid(), (psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024 / 1024))) del chip.img - del chip.flat_img - del chip.prnu_img - del chip.shutter_img def runExposure_MPI_PointingList(self, pointing_list, chips=None): comm = MPI.COMM_WORLD diff --git a/observation_sim/config/ChipOutput.py b/observation_sim/config/ChipOutput.py index 7db5e12..5efd061 100755 --- a/observation_sim/config/ChipOutput.py +++ b/observation_sim/config/ChipOutput.py @@ -78,12 +78,12 @@ class ChipOutput(object): self.hdr += "\n" self.cat.write(self.hdr) - def cat_add_obj(self, obj, pos_img, pos_shear, ra_offset=0., dec_offset=0.): + def cat_add_obj(self, obj, pos_img, pos_shear): ximg = obj.real_pos.x + 1.0 yimg = obj.real_pos.y + 1.0 line = self.fmt % ( - obj.id, int(self.chip_label), self.filt.filter_type, ximg, yimg, obj.ra + ra_offset, obj.dec + dec_offset, obj.ra_orig, obj.dec_orig, obj.z, obj.getMagFilter( + obj.id, int(self.chip_label), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.ra_orig, obj.dec_orig, obj.z, obj.getMagFilter( self.filt), obj.type, obj.pmra, obj.pmdec, obj.rv, obj.parallax) line += obj.additional_output_str diff --git a/observation_sim/config/header/ImageHeader.py b/observation_sim/config/header/ImageHeader.py index 25a2867..7fc3444 100644 --- a/observation_sim/config/header/ImageHeader.py +++ b/observation_sim/config/header/ImageHeader.py @@ -434,7 +434,7 @@ def generatePrimaryHeader(xlen=9216, ylen=9232, pointing_id='00000001', pointing co = coord.SkyCoord(ra, dec, unit='deg') ra_hms = format(co.ra.hms.h, '02.0f') + format(co.ra.hms.m, - '02.0f') + format(co.ra.hms.s, '04.1f') + '02.0f') + format(co.ra.hms.s, '02.1f') dec_hms = format(co.dec.dms.d, '02.0f') + format(abs(co.dec.dms.m), '02.0f') + format(abs(co.dec.dms.s), '02.0f') if dec >= 0: diff --git a/observation_sim/instruments/chip/effects.py b/observation_sim/instruments/chip/effects.py index 21735ed..7a6b7a6 100644 --- a/observation_sim/instruments/chip/effects.py +++ b/observation_sim/instruments/chip/effects.py @@ -793,7 +793,7 @@ def ShutterEffectArr(GSImage, t_exp=150, t_shutter=1.3, dist_bearing=735, dt=1E- sizey = ymax-ymin+1 xnewgrid = np.mgrid[xmin:(xmin+sizex)] expeffect = interpolate.splev(xnewgrid, intp, der=0) - expeffect /= t_exp + expeffect /= np.max(expeffect) exparrnormal = np.tile(expeffect, (sizey,1)) # GSImage *= exparrnormal diff --git a/observation_sim/instruments/data/filters/fgs_sub.list b/observation_sim/instruments/data/filters/fgs_sub.list index c02a44f..5a826ee 100644 --- a/observation_sim/instruments/data/filters/fgs_sub.list +++ b/observation_sim/instruments/data/filters/fgs_sub.list @@ -1,4 +1,4 @@ -2000 +3000 4500 4750 5000 diff --git a/observation_sim/instruments/data/filters/g_sub.list b/observation_sim/instruments/data/filters/g_sub.list index 1384ade..f5dda57 100644 --- a/observation_sim/instruments/data/filters/g_sub.list +++ b/observation_sim/instruments/data/filters/g_sub.list @@ -1,4 +1,4 @@ -2000 +3800 4217 4432 4631 @@ -6,4 +6,4 @@ 5002 5179 5354 -11000 \ No newline at end of file +5799 \ No newline at end of file diff --git a/observation_sim/instruments/data/filters/i_sub.list b/observation_sim/instruments/data/filters/i_sub.list index 0449403..715f74a 100644 --- a/observation_sim/instruments/data/filters/i_sub.list +++ b/observation_sim/instruments/data/filters/i_sub.list @@ -1,4 +1,4 @@ -2000 +6600 7061 7255 7448 @@ -6,4 +6,4 @@ 7833 8027 8226 -11000 \ No newline at end of file +8999 \ No newline at end of file diff --git a/observation_sim/instruments/data/filters/nuv_sub.list b/observation_sim/instruments/data/filters/nuv_sub.list index b6c60c5..c89b5dc 100644 --- a/observation_sim/instruments/data/filters/nuv_sub.list +++ b/observation_sim/instruments/data/filters/nuv_sub.list @@ -1,4 +1,4 @@ -2000 +2513 2621 2716 2805 @@ -6,4 +6,4 @@ 2969 3050 3132 -11000 \ No newline at end of file +3499 \ No newline at end of file diff --git a/observation_sim/instruments/data/filters/r_sub.list b/observation_sim/instruments/data/filters/r_sub.list index ec0852e..0abf191 100644 --- a/observation_sim/instruments/data/filters/r_sub.list +++ b/observation_sim/instruments/data/filters/r_sub.list @@ -1,4 +1,4 @@ -2000 +5100 5642 5821 6001 @@ -6,4 +6,4 @@ 6363 6547 6735 -11000 \ No newline at end of file +7199 \ No newline at end of file diff --git a/observation_sim/instruments/data/filters/u_sub.list b/observation_sim/instruments/data/filters/u_sub.list index 3ec5a58..f3a2b72 100644 --- a/observation_sim/instruments/data/filters/u_sub.list +++ b/observation_sim/instruments/data/filters/u_sub.list @@ -1,4 +1,4 @@ -2000 +3000 3277 3380 3485 @@ -6,4 +6,4 @@ 3703 3813 3918 -11000 \ No newline at end of file +4499 \ No newline at end of file diff --git a/observation_sim/instruments/data/filters/y_sub.list b/observation_sim/instruments/data/filters/y_sub.list index 69a05ee..81dc243 100644 --- a/observation_sim/instruments/data/filters/y_sub.list +++ b/observation_sim/instruments/data/filters/y_sub.list @@ -1,4 +1,4 @@ -2000 +9000 9322 9405 9489 @@ -6,4 +6,4 @@ 9695 9832 10024 -11000 \ No newline at end of file +10590 \ No newline at end of file diff --git a/observation_sim/instruments/data/filters/z_sub.list b/observation_sim/instruments/data/filters/z_sub.list index f105711..b0191a5 100644 --- a/observation_sim/instruments/data/filters/z_sub.list +++ b/observation_sim/instruments/data/filters/z_sub.list @@ -1,4 +1,4 @@ -2000 +7800 8494 8638 8790 @@ -6,4 +6,4 @@ 9141 9363 9663 -11000 \ No newline at end of file +10590 \ No newline at end of file diff --git a/observation_sim/mock_objects/ExtinctionMW.py b/observation_sim/mock_objects/ExtinctionMW.py deleted file mode 100644 index 85ae692..0000000 --- a/observation_sim/mock_objects/ExtinctionMW.py +++ /dev/null @@ -1,96 +0,0 @@ -import os -import numpy as np -import healpy as hp -from astropy.io import fits -from astropy import units as u -from astropy.coordinates import SkyCoord - - -class ExtinctionMW(object): - def __init__(self): - self.Rv = 3.1 - self.lamb = np.arange(2000, 11001+0.5, 0.5) - - @staticmethod - def radec2pix(ra, dec, NSIDE=2048): - return hp.pixelfunc.ang2pix(nside=NSIDE, theta=ra, phi=dec, nest=True, lonlat=True) - - def init_ext_model(self, model_name="odonnell", Av=1.0, Rv=3.1, lamb=None): - self.model_name = model_name - self.Av = Av - self.Rv = Rv - if lamb is not None: - self.lamb = lamb - - if self.model_name == "odonnell": - alpha = np.zeros(self.lamb.shape) - beta = np.zeros(self.lamb.shape) - x = 1.e4 / self.lamb - def f_alpha_1(x): return 0.574 * (x**1.61) - def f_beta_1(x): return -0.527 * (x**1.61) - - def f_alpha_2(x): - y = x - 1.82 - return 1 + 0.104*y - 0.609*(y**2) + 0.701*(y**3) + \ - 1.137*(y**4) - 1.718*(y**5) - 0.827 * \ - (y**6) + 1.647*(y**7) - 0.505*(y**8) - - def f_beta_2(x): - y = x - 1.82 - return 1.952*y + 2.908*(y**2) - 3.989*(y**3) - 7.985 * \ - (y**4) + 11.102*(y**5) + 5.491 * \ - (y**6) - 10.805*(y**7) + 3.347*(y**8) - - def f_alpha_3(x): return 1.752 - 0.316*x - \ - 0.104 / ((x - 4.67)**2 + 0.341) - - def f_beta_3(x): return -3.090 + 1.825*x + \ - 1.206 / ((x - 4.62)**2 + 0.262) - - def f_alpha_4(x): return f_alpha_3(x) + -0.04473 * \ - (x - 5.9)**2 - 0.009779 * (x - 5.9)**3 - - def f_beta_4(x): return f_beta_3(x) + 0.2130 * \ - (x - 5.9)**2 + 0.1207 * (x - 5.9)**3 - - def f_alpha_5(x): return -1.073 - 0.628*(x - 8) + \ - 0.137*(x - 8)**2 - 0.070*(x - 8)**3 - - def f_beta_5(x): return 13.670 + 4.257*(x - 8) - \ - 0.420*(x - 8)**2 + 0.374 * (x - 8)**3 - - alpha = np.piecewise( - x, [x <= 1.1, (x > 1.1)*(x <= 3.3), (x > 3.3)*(x <= 5.9), (x > 5.9)*(x <= 8.), (x > 8.)], [f_alpha_1, f_alpha_2, f_alpha_3, f_alpha_4, f_alpha_5]) - beta = np.piecewise( - x, [x <= 1.1, (x > 1.1)*(x <= 3.3), (x > 3.3)*(x <= 5.9), (x > 5.9)*(x <= 8.), (x > 8.)], [f_beta_1, f_beta_2, f_beta_3, f_beta_4, f_beta_5]) - - self.ext = (alpha + beta / self.Rv) * self.Av - - def load_Planck_ext(self, file_path): - hdu = fits.open(file_path) - self.ebv_planck = hdu[1].data["EBV"] - - def Av_from_Planck(self, ra, dec): - if not hasattr(self, 'ebv_planck'): - raise ValueError( - "Need to load planck dust map first") - # Convert to galactic coordinates - c = SkyCoord(ra=ra * u.degree, dec=dec * - u.degree, frame='icrs').galactic - l, b = c.l.radian, c.b.radian - NSIDE = hp.pixelfunc.get_nside(self.ebv_planck) - pix = hp.pixelfunc.ang2pix( - nside=NSIDE, theta=np.pi/2. - b, phi=l, nest=True) - return self.ebv_planck[pix] * self.Rv - - def apply_extinction(self, spec, Av=1.0): - if len(spec) != len(self.lamb): - raise ValueError( - "Dimension of spec do not match that of ExtinctionMW.lamb: try initilize (init_ext_model) by the actual size of spec which you want to apply the extinction to") - if not hasattr(self, 'ext'): - raise ValueError( - "Need to initialize the extinction model (init_ext_model) first") - scale = 10**(-.4 * self.ext * Av) - # print("scale = ", scale) - spec *= scale - return spec diff --git a/observation_sim/mock_objects/Galaxy.py b/observation_sim/mock_objects/Galaxy.py index b0404e7..c220621 100755 --- a/observation_sim/mock_objects/Galaxy.py +++ b/observation_sim/mock_objects/Galaxy.py @@ -323,73 +323,28 @@ class Galaxy(MockObject): # # if fd_shear is not None: # # gal = gal.shear(fd_shear) - galImg_List = [] - try: - pos_img_local = [0,0] - x_start = chip.x_cen/chip.pix_size - chip.npix_x / 2. - y_start = chip.y_cen/chip.pix_size - chip.npix_y / 2. - pos_img_local[0] = pos_img.x - x_start - pos_img_local[1] = pos_img.y - y_start - nnx = 0 - nny = 0 - for order in ["A","B"]: - psf, pos_shear = psf_model.get_PSF( - chip, pos_img_local=pos_img_local, bandNo=i+1, galsimGSObject=True, g_order=order, grating_split_pos=grating_split_pos) - star_p = galsim.Convolve(psf, gal) - if nnx == 0: - galImg = star_p.drawImage(wcs=chip_wcs_local, offset=offset) - nnx = galImg.xmax - galImg.xmin + 1 - nny = galImg.ymax - galImg.ymin + 1 - else: - galImg = star_p.drawImage(nx = nnx, ny = nny, wcs=chip_wcs_local, offset=offset) - galImg.setOrigin(0, 0) - # n1 = np.sum(np.isinf(galImg.array)) - # n2 = np.sum(np.isnan(galImg.array)) - # if n1>0 or n2 > 0: - # print("DEBUG: Galaxy, inf:%d, nan:%d"%(n1, n2)) - if np.sum(np.isnan(galImg.array)) > 0: - # ERROR happens - return 2, pos_shear - galImg_List.append(galImg) - for order in ["C","D","E"]: - galImg_List.append(galImg) - except: - psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img) - star_p = galsim.Convolve(psf, gal) - galImg = star_p.drawImage(wcs=chip_wcs_local, offset=offset) - galImg.setOrigin(0, 0) - if np.sum(np.isnan(galImg.array)) > 0: - # ERROR happens - return 2, pos_shear - for order in ["A","B","C","D","E"]: - galImg_List.append(galImg) - - # starImg = gal.drawImage( - # wcs=chip_wcs_local, offset=offset, method='real_space') - - origin_star = [y_nominal - (galImg.center.y - galImg.ymin), - x_nominal - (galImg.center.x - galImg.xmin)] - galImg.setOrigin(0, 0) + starImg = gal.drawImage( + wcs=chip_wcs_local, offset=offset, method='real_space') + + origin_star = [y_nominal - (starImg.center.y - starImg.ymin), + x_nominal - (starImg.center.x - starImg.xmin)] + starImg.setOrigin(0, 0) gal_origin = [origin_star[0], origin_star[1]] - gal_end = [origin_star[0] + galImg.array.shape[0] - - 1, origin_star[1] + galImg.array.shape[1] - 1] + gal_end = [origin_star[0] + starImg.array.shape[0] - + 1, origin_star[1] + starImg.array.shape[1] - 1] if gal_origin[1] < grating_split_pos_chip < gal_end[1]: subSlitPos = int(grating_split_pos_chip - gal_origin[1] + 1) # part img disperse - star_p1s=[] - for galImg in galImg_List: - - subImg_p1 = galImg.array[:, 0:subSlitPos] - star_p1 = galsim.Image(subImg_p1) - star_p1.setOrigin(0, 0) - star_p1s.append(star_p1) + subImg_p1 = starImg.array[:, 0:subSlitPos] + star_p1 = galsim.Image(subImg_p1) + star_p1.setOrigin(0, 0) origin_p1 = origin_star xcenter_p1 = min(x_nominal, grating_split_pos_chip-1) - 0 ycenter_p1 = y_nominal-0 - sdp_p1 = SpecDisperser(orig_img=star_p1s, xcenter=xcenter_p1, + sdp_p1 = SpecDisperser(orig_img=star_p1, xcenter=xcenter_p1, ycenter=ycenter_p1, origin=origin_p1, tar_spec=normalSED, band_start=brange[0], band_end=brange[1], @@ -397,25 +352,21 @@ class Galaxy(MockObject): isAlongY=0, flat_cube=flat_cube) - self.addSLStoChipImage(sdp=sdp_p1, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local) - # pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp_p1, chip=chip, pos_img_local=[xcenter_p1, ycenter_p1], - # psf_model=psf_model, bandNo=i + 1, - # grating_split_pos=grating_split_pos, - # local_wcs=chip_wcs_local, pos_img=pos_img) - - star_p2s=[] - for galImg in galImg_List: + # self.addSLStoChipImage(sdp=sdp_p1, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local) + pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp_p1, chip=chip, pos_img_local=[xcenter_p1, ycenter_p1], + psf_model=psf_model, bandNo=i + 1, + grating_split_pos=grating_split_pos, + local_wcs=chip_wcs_local, pos_img=pos_img) - subImg_p2 = galImg.array[:, - subSlitPos + 1:galImg.array.shape[1]] - star_p2 = galsim.Image(subImg_p2) - star_p2.setOrigin(0, 0) - star_p2s.append(star_p2) + subImg_p2 = starImg.array[:, + subSlitPos+1:starImg.array.shape[1]] + star_p2 = galsim.Image(subImg_p2) + star_p2.setOrigin(0, 0) origin_p2 = [origin_star[0], grating_split_pos_chip] xcenter_p2 = max(x_nominal, grating_split_pos_chip - 1) - 0 ycenter_p2 = y_nominal - 0 - sdp_p2 = SpecDisperser(orig_img=star_p2s, xcenter=xcenter_p2, + sdp_p2 = SpecDisperser(orig_img=star_p2, xcenter=xcenter_p2, ycenter=ycenter_p2, origin=origin_p2, tar_spec=normalSED, band_start=brange[0], band_end=brange[1], @@ -423,41 +374,41 @@ class Galaxy(MockObject): isAlongY=0, flat_cube=flat_cube) - self.addSLStoChipImage(sdp=sdp_p2, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local) - # pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp_p2, chip=chip, pos_img_local=[xcenter_p2, ycenter_p2], - # psf_model=psf_model, bandNo=i + 1, - # grating_split_pos=grating_split_pos, - # local_wcs=chip_wcs_local, pos_img=pos_img) + # self.addSLStoChipImage(sdp=sdp_p2, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local) + pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp_p2, chip=chip, pos_img_local=[xcenter_p2, ycenter_p2], + psf_model=psf_model, bandNo=i + 1, + grating_split_pos=grating_split_pos, + local_wcs=chip_wcs_local, pos_img=pos_img) del sdp_p1 del sdp_p2 elif grating_split_pos_chip <= gal_origin[1]: - sdp = SpecDisperser(orig_img=galImg_List, xcenter=x_nominal - 0, + sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0, ycenter=y_nominal - 0, origin=origin_star, tar_spec=normalSED, band_start=brange[0], band_end=brange[1], conf=chip.sls_conf[1], isAlongY=0, flat_cube=flat_cube) - self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local) - # pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal], - # psf_model=psf_model, bandNo=i + 1, - # grating_split_pos=grating_split_pos, - # local_wcs=chip_wcs_local, pos_img=pos_img) + # self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local) + pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal], + psf_model=psf_model, bandNo=i + 1, + grating_split_pos=grating_split_pos, + local_wcs=chip_wcs_local, pos_img=pos_img) del sdp elif grating_split_pos_chip >= gal_end[1]: - sdp = SpecDisperser(orig_img=galImg_List, xcenter=x_nominal - 0, + sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0, ycenter=y_nominal - 0, origin=origin_star, tar_spec=normalSED, band_start=brange[0], band_end=brange[1], conf=chip.sls_conf[0], isAlongY=0, flat_cube=flat_cube) - self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local) - # pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal], - # psf_model=psf_model, bandNo=i + 1, - # grating_split_pos=grating_split_pos, - # local_wcs=chip_wcs_local, pos_img=pos_img) + # self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local) + pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal], + psf_model=psf_model, bandNo=i + 1, + grating_split_pos=grating_split_pos, + local_wcs=chip_wcs_local, pos_img=pos_img) del sdp # print(self.y_nominal, starImg.center.y, starImg.ymin) diff --git a/observation_sim/mock_objects/MockObject.py b/observation_sim/mock_objects/MockObject.py index afbbc37..12e3612 100755 --- a/observation_sim/mock_objects/MockObject.py +++ b/observation_sim/mock_objects/MockObject.py @@ -11,8 +11,6 @@ from observation_sim.mock_objects._util import integrate_sed_bandpass, getNormFa getABMAG from observation_sim.mock_objects.SpecDisperser import SpecDisperser -from observation_sim.instruments.chip import chip_utils - class MockObject(object): def __init__(self, param, logger=None): @@ -61,13 +59,13 @@ class MockObject(object): flux = self.getFluxFilter(filt) return flux * tel.pupil_area * exptime - def getPosWorld(self, ra_offset=0., dec_offset=0.): - ra = self.param["ra"] + ra_offset - dec = self.param["dec"] + dec_offset + def getPosWorld(self): + ra = self.param["ra"] + dec = self.param["dec"] return galsim.CelestialCoord(ra=ra * galsim.degrees, dec=dec * galsim.degrees) - def getPosImg_Offset_WCS(self, img, fdmodel=None, chip=None, verbose=True, chip_wcs=None, img_header=None, ra_offset=0., dec_offset=0.): - self.posImg = img.wcs.toImage(self.getPosWorld(ra_offset, dec_offset)) + def getPosImg_Offset_WCS(self, img, fdmodel=None, chip=None, verbose=True, chip_wcs=None, img_header=None): + self.posImg = img.wcs.toImage(self.getPosWorld()) self.localWCS = img.wcs.local(self.posImg) # Apply field distortion model if (fdmodel is not None) and (chip is not None): @@ -120,6 +118,7 @@ class MockObject(object): if self.logger: self.logger.error(e) return 2, None + # Set Galsim Parameters if self.getMagFilter(filt) <= 15: folding_threshold = 5.e-4 @@ -196,26 +195,33 @@ class MockObject(object): # 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]) + 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]) - orig_off = 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]) + 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]) ######################################################### - stamp = galsim.ImageF(img_s) + specImg = galsim.ImageF(img_s) + photons = galsim.PhotonArray.makeFromImage(specImg) + photons.x += origin_order_x + photons.y += origin_order_y + + xlen_imf = int(specImg.xmax - specImg.xmin + 1) + ylen_imf = int(specImg.ymax - specImg.ymin + 1) + stamp = galsim.ImageF(xlen_imf, ylen_imf) stamp.wcs = local_wcs stamp.setOrigin(origin_order_x, origin_order_y) @@ -224,140 +230,72 @@ class MockObject(object): if bounds.area() == 0: continue chip.img.setOrigin(0, 0) - chip.img[bounds] = chip.img[bounds]+stamp[bounds] + stamp[bounds] = chip.img[bounds] + chip.sensor.accumulate(photons, stamp) + chip.img[bounds] = stamp[bounds] chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin) del stamp del spec_orders def addSLStoChipImageWithPSF(self, sdp=None, chip=None, pos_img_local=[1, 1], psf_model=None, bandNo=1, grating_split_pos=3685, local_wcs=None, pos_img=None): spec_orders = sdp.compute_spec_orders() - pos_shear = galsim.Shear(e=0., beta=(np.pi/2)*galsim.radians) - if chip.slsPSFOptim: - for k, v in spec_orders.items(): - img_s = v[0] - - 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]) - ######################################################### - # img_s, orig_off = convolveImg(img_s, psf_img_m) - orig_off = [0,0] - origin_order_x = v[1] - orig_off[0] - origin_order_y = v[2] - orig_off[1] - - specImg = galsim.ImageF(img_s) - - specImg.wcs = local_wcs - specImg.setOrigin(origin_order_x, origin_order_y) - - bounds = specImg.bounds & galsim.BoundsI( - 0, chip.npix_x - 1, 0, chip.npix_y - 1) - if bounds.area() == 0: - continue - - # orders = {'A': 'order1', 'B': 'order0', 'C': 'order2', 'D': 'order-1', 'E': 'order-2'} - orders = {'A': 'order1', 'B': 'order0', 'C': 'order0', 'D': 'order0', 'E': 'order0'} - gratingN = chip_utils.getChipSLSGratingID(chip.chipID)[1] - if pos_img_local[0] < grating_split_pos: - gratingN = chip_utils.getChipSLSGratingID(chip.chipID)[0] - - - chip.img_stack[gratingN][orders[k]]['w' + str(bandNo)].setOrigin(0, 0) - chip.img_stack[gratingN][orders[k]]['w' + str(bandNo)][bounds] = chip.img_stack[gratingN][orders[k]]['w' + str(bandNo)][bounds] + specImg[bounds] - chip.img_stack[gratingN][orders[k]]['w' + str(bandNo)].setOrigin(chip.bound.xmin, chip.bound.ymin) + for k, v in spec_orders.items(): + img_s = v[0] + # print(bandNo,k) + try: + psf, pos_shear = psf_model.get_PSF( + chip, pos_img_local=pos_img_local, bandNo=bandNo, galsimGSObject=True, g_order=k, grating_split_pos=grating_split_pos) + except: + psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img) - else: - for k, v in spec_orders.items(): - # img_s = v[0] - # # print(bandNo,k) - # try: - # psf, pos_shear = psf_model.get_PSF( - # chip, pos_img_local=pos_img_local, bandNo=bandNo, galsimGSObject=True, g_order=k, grating_split_pos=grating_split_pos) - # except: - # psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img) - - # psf_img = psf.drawImage(nx=100, ny=100, wcs=local_wcs) - - # psf_img_m = psf_img.array - - # ######################################################### - # # DEBUG - # ######################################################### - # # ids_p = psf_img_m < 0 - # # psf_img_m[ids_p] = 0 - - # # from astropy.io import fits - # # fits.writeto(str(bandNo) + '_' + str(k) + '_psf.fits', psf_img_m) - - # # 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]) - # ######################################################### - # img_s, orig_off = convolveImg(img_s, psf_img_m) - # origin_order_x = v[1] - orig_off[0] - # origin_order_y = v[2] - orig_off[1] - - # specImg = galsim.ImageF(img_s) - # # photons = galsim.PhotonArray.makeFromImage(specImg) - # # photons.x += origin_order_x - # # photons.y += origin_order_y - - # # xlen_imf = int(specImg.xmax - specImg.xmin + 1) - # # ylen_imf = int(specImg.ymax - specImg.ymin + 1) - # # stamp = galsim.ImageF(xlen_imf, ylen_imf) - # # stamp.wcs = local_wcs - # # stamp.setOrigin(origin_order_x, origin_order_y) - - # specImg.wcs = local_wcs - # specImg.setOrigin(origin_order_x, origin_order_y) - - # print('DEBUG: BEGIN -----------',bandNo,k) - - img_s = v[0] - - 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]) - ######################################################### - origin_order_x = v[1] - origin_order_y = v[2] - specImg = galsim.ImageF(img_s) - - specImg.wcs = local_wcs - specImg.setOrigin(origin_order_x, origin_order_y) - try: - specImg = psf_model.get_PSF_AND_convolve_withsubImg(chip, cutImg=specImg, pos_img_local=pos_img_local, bandNo=bandNo, g_order=k, grating_split_pos=grating_split_pos) - except: - psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img) - - psf_img = psf.drawImage(nx=100, ny=100, wcs=local_wcs) - - psf_img_m = psf_img.array - - img_s, orig_off = convolveImg(img_s, psf_img_m) - origin_order_x = v[1] - orig_off[0] - origin_order_y = v[2] - orig_off[1] - - specImg = galsim.ImageF(img_s) - - specImg.wcs = local_wcs - specImg.setOrigin(origin_order_x, origin_order_y) - - bounds = specImg.bounds & galsim.BoundsI( - 0, chip.npix_x - 1, 0, chip.npix_y - 1) - if bounds.area() == 0: - continue - chip.img.setOrigin(0, 0) - chip.img[bounds] = chip.img[bounds] + specImg[bounds] - # stamp[bounds] = chip.img[bounds] - # # chip.sensor.accumulate(photons, stamp) - # chip.img[bounds] = stamp[bounds] - chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin) - # del stamp + psf_img = psf.drawImage(nx=100, ny=100, wcs=local_wcs) + + psf_img_m = psf_img.array + + ######################################################### + # DEBUG + ######################################################### + # ids_p = psf_img_m < 0 + # psf_img_m[ids_p] = 0 + + # from astropy.io import fits + # fits.writeto(str(bandNo) + '_' + str(k) + '_psf.fits', psf_img_m) + + # 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]) + ######################################################### + img_s, orig_off = convolveImg(img_s, psf_img_m) + origin_order_x = v[1] - orig_off[0] + origin_order_y = v[2] - orig_off[1] + + specImg = galsim.ImageF(img_s) + # photons = galsim.PhotonArray.makeFromImage(specImg) + # photons.x += origin_order_x + # photons.y += origin_order_y + + # xlen_imf = int(specImg.xmax - specImg.xmin + 1) + # ylen_imf = int(specImg.ymax - specImg.ymin + 1) + # stamp = galsim.ImageF(xlen_imf, ylen_imf) + # stamp.wcs = local_wcs + # stamp.setOrigin(origin_order_x, origin_order_y) + + specImg.wcs = local_wcs + specImg.setOrigin(origin_order_x, origin_order_y) + + bounds = specImg.bounds & galsim.BoundsI( + 0, chip.npix_x - 1, 0, chip.npix_y - 1) + if bounds.area() == 0: + continue + chip.img.setOrigin(0, 0) + chip.img[bounds] = chip.img[bounds] + specImg[bounds] + # stamp[bounds] = chip.img[bounds] + # # chip.sensor.accumulate(photons, stamp) + # chip.img[bounds] = stamp[bounds] + chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin) + # del stamp del spec_orders return pos_shear @@ -421,76 +359,32 @@ class MockObject(object): # psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, # folding_threshold=folding_threshold) + star = galsim.DeltaFunction(gsparams=gsp) + star = star.withFlux(tel.pupil_area * exptime) + psf_tmp = galsim.Gaussian(sigma=0.002) + star = galsim.Convolve(psf_tmp, star) - # star = galsim.DeltaFunction(gsparams=gsp) - # star = star.withFlux(tel.pupil_area * exptime) - - #psf list :["A","B","C","D","E"] - starImg_List = [] - try: - pos_img_local = [0,0] - x_start = chip.x_cen/chip.pix_size - chip.npix_x / 2. - y_start = chip.y_cen/chip.pix_size - chip.npix_y / 2. - pos_img_local[0] = pos_img.x - x_start - pos_img_local[1] = pos_img.y - y_start - nnx = 0 - nny = 0 - for order in ["A","B"]: - psf, pos_shear = psf_model.get_PSF( - chip, pos_img_local=pos_img_local, bandNo=i+1, galsimGSObject=True, g_order=order, grating_split_pos=grating_split_pos) - # star_p = galsim.Convolve(psf, star) - star_p = psf.withFlux(tel.pupil_area * exptime) - if nnx == 0: - starImg = star_p.drawImage(wcs=chip_wcs_local, offset=offset) - nnx = starImg.xmax - starImg.xmin + 1 - nny = starImg.ymax - starImg.ymin + 1 - else: - starImg = star_p.drawImage(nx = nnx, ny = nny, wcs=chip_wcs_local, offset=offset) - # n1 = np.sum(np.isinf(starImg.array)) - # n2 = np.sum(np.isnan(starImg.array)) - # if n1>0 or n2 > 0: - # print("DEBUG: MockObject, inf:%d, nan:%d"%(n1, n2)) - starImg.setOrigin(0, 0) - starImg_List.append(starImg) - for order in ["C","D","E"]: - starImg_List.append(starImg) - except: - psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img) - # star_p = galsim.Convolve(psf, star) - star_p = psf.withFlux(tel.pupil_area * exptime) - starImg = star_p.drawImage(wcs=chip_wcs_local, offset=offset) - starImg.setOrigin(0, 0) - for order in ["A","B","C","D","E"]: - starImg_List.append(starImg) - - - # psf_tmp = galsim.Gaussian(sigma=0.002) - # star = galsim.Convolve(psf_tmp, star) - - # starImg = star.drawImage( - # nx=60, ny=60, wcs=chip_wcs_local, offset=offset) + starImg = star.drawImage( + nx=60, ny=60, wcs=chip_wcs_local, offset=offset) origin_star = [y_nominal - (starImg.center.y - starImg.ymin), x_nominal - (starImg.center.x - starImg.xmin)] - + starImg.setOrigin(0, 0) gal_origin = [origin_star[0], origin_star[1]] gal_end = [origin_star[0] + starImg.array.shape[0] - 1, origin_star[1] + starImg.array.shape[1] - 1] if gal_origin[1] < grating_split_pos_chip < gal_end[1]: subSlitPos = int(grating_split_pos_chip - gal_origin[1] + 1) # part img disperse - star_p1s=[] - for starImg in starImg_List: - subImg_p1 = starImg.array[:, 0:subSlitPos] - star_p1 = galsim.Image(subImg_p1) - star_p1.setOrigin(0, 0) - star_p1s.append(star_p1) + subImg_p1 = starImg.array[:, 0:subSlitPos] + star_p1 = galsim.Image(subImg_p1) origin_p1 = origin_star + star_p1.setOrigin(0, 0) xcenter_p1 = min(x_nominal, grating_split_pos_chip - 1) - 0 ycenter_p1 = y_nominal - 0 - sdp_p1 = SpecDisperser(orig_img=star_p1s, xcenter=xcenter_p1, + sdp_p1 = SpecDisperser(orig_img=star_p1, xcenter=xcenter_p1, ycenter=ycenter_p1, origin=origin_p1, tar_spec=normalSED, band_start=brange[0], band_end=brange[1], @@ -498,25 +392,20 @@ class MockObject(object): isAlongY=0, flat_cube=flat_cube) - self.addSLStoChipImage(sdp=sdp_p1, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local) - # pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp_p1, chip=chip, pos_img_local=[xcenter_p1, ycenter_p1], - # psf_model=psf_model, bandNo=i+1, grating_split_pos=grating_split_pos, - # local_wcs=chip_wcs_local, pos_img=pos_img) - - - star_p2s=[] - for starImg in starImg_List: + # self.addSLStoChipImage(sdp=sdp_p1, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local) + pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp_p1, chip=chip, pos_img_local=[xcenter_p1, ycenter_p1], + psf_model=psf_model, bandNo=i+1, grating_split_pos=grating_split_pos, + local_wcs=chip_wcs_local, pos_img=pos_img) - subImg_p2 = starImg.array[:, + subImg_p2 = starImg.array[:, subSlitPos + 1:starImg.array.shape[1]] - star_p2 = galsim.Image(subImg_p2) - star_p2.setOrigin(0, 0) - star_p2s.append(star_p2) + star_p2 = galsim.Image(subImg_p2) + star_p2.setOrigin(0, 0) origin_p2 = [origin_star[0], grating_split_pos_chip] xcenter_p2 = max(x_nominal, grating_split_pos_chip - 1) - 0 ycenter_p2 = y_nominal - 0 - sdp_p2 = SpecDisperser(orig_img=star_p2s, xcenter=xcenter_p2, + sdp_p2 = SpecDisperser(orig_img=star_p2, xcenter=xcenter_p2, ycenter=ycenter_p2, origin=origin_p2, tar_spec=normalSED, band_start=brange[0], band_end=brange[1], @@ -524,38 +413,38 @@ class MockObject(object): isAlongY=0, flat_cube=flat_cube) - self.addSLStoChipImage(sdp=sdp_p2, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local) - # pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp_p2, chip=chip, pos_img_local=[xcenter_p2, ycenter_p2], - # psf_model=psf_model, bandNo=i + 1, grating_split_pos=grating_split_pos, - # local_wcs=chip_wcs_local, pos_img=pos_img) + # self.addSLStoChipImage(sdp=sdp_p2, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local) + pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp_p2, chip=chip, pos_img_local=[xcenter_p2, ycenter_p2], + psf_model=psf_model, bandNo=i + 1, grating_split_pos=grating_split_pos, + local_wcs=chip_wcs_local, pos_img=pos_img) del sdp_p1 del sdp_p2 elif grating_split_pos_chip <= gal_origin[1]: - sdp = SpecDisperser(orig_img=starImg_List, xcenter=x_nominal - 0, + sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0, ycenter=y_nominal - 0, origin=origin_star, tar_spec=normalSED, band_start=brange[0], band_end=brange[1], conf=chip.sls_conf[1], isAlongY=0, flat_cube=flat_cube) - self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local) - # pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal], - # psf_model=psf_model, bandNo=i + 1, grating_split_pos=grating_split_pos, - # local_wcs=chip_wcs_local, pos_img=pos_img) + # self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local) + pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal], + psf_model=psf_model, bandNo=i + 1, grating_split_pos=grating_split_pos, + local_wcs=chip_wcs_local, pos_img=pos_img) del sdp elif grating_split_pos_chip >= gal_end[1]: - sdp = SpecDisperser(orig_img=starImg_List, xcenter=x_nominal - 0, + sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0, ycenter=y_nominal - 0, origin=origin_star, tar_spec=normalSED, band_start=brange[0], band_end=brange[1], conf=chip.sls_conf[0], isAlongY=0, flat_cube=flat_cube) - self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local) - # pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal], - # psf_model=psf_model, bandNo=i + 1, grating_split_pos=grating_split_pos, - # local_wcs=chip_wcs_local, pos_img=pos_img) + # self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local) + pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal], + psf_model=psf_model, bandNo=i + 1, grating_split_pos=grating_split_pos, + local_wcs=chip_wcs_local, pos_img=pos_img) del sdp # del psf return 1, pos_shear diff --git a/observation_sim/mock_objects/SpecDisperser/SpecDisperser.py b/observation_sim/mock_objects/SpecDisperser/SpecDisperser.py index 2b2d7c2..6a1e3ef 100644 --- a/observation_sim/mock_objects/SpecDisperser/SpecDisperser.py +++ b/observation_sim/mock_objects/SpecDisperser/SpecDisperser.py @@ -65,30 +65,10 @@ class SpecDisperser(object): # self.img_x = orig_img.shape[1] # self.img_y = orig_img.shape[0] - # 5 orders, A, B , - orderName=["A","B","C","D","E"] - self.orig_img_orders = OrderedDict() - if isinstance(orig_img, list): - orig_img_list = orig_img - list_len = len(orig_img_list) - if list_len < 5: - for i in np.arange(5-list_len): - orig_img_list.append(orig_img_list[list_len-1]) - for i, k in enumerate(orig_img_list): - self.orig_img_orders[orderName[i]] = k - - - if isinstance(orig_img, galsim.Image): - for i in np.arange(5): - self.orig_img_orders[orderName[i]] = orig_img - - - orig_img_one = self.orig_img_orders["A"] - - self.thumb_img = np.abs(orig_img_one.array) - self.thumb_x = orig_img_one.center.x - self.thumb_y = orig_img_one.center.y - self.img_sh = orig_img_one.array.shape + self.thumb_img = np.abs(orig_img.array) + self.thumb_x = orig_img.center.x + self.thumb_y = orig_img.center.y + self.img_sh = orig_img.array.shape self.id = gid @@ -98,13 +78,10 @@ class SpecDisperser(object): self.isAlongY = isAlongY self.flat_cube = flat_cube if self.isAlongY == 1: - for order in orderName: - self.orig_img_orders[order], self.thumb_x, self.thumb_y = rotate90(array_orig=self.orig_img_orders[order], xc=orig_img_one.center.x, - yc=orig_img_one.center.y, isClockwise=1) - # self.thumb_img, self.thumb_x, self.thumb_y = rotate90(array_orig=self.thumb_img, xc=orig_img_one.center.x, - # yc=orig_img_one.center.y, isClockwise=1) + self.thumb_img, self.thumb_x, self.thumb_y = rotate90(array_orig=self.thumb_img, xc=orig_img.center.x, + yc=orig_img.center.y, isClockwise=1) - self.img_sh = self.orig_img_orders[order].array.T.shape + self.img_sh = orig_img.array.T.shape self.xcenter = ycenter self.ycenter = xcenter @@ -134,16 +111,10 @@ class SpecDisperser(object): def compute_spec(self, beam): - # if beam == "B": - # return self.thumb_img, self.origin[1], self.origin[0], None, None, None - from .disperse_c import interp from .disperse_c import disperse # from MockObject.disperse_c import disperse - self.thumb_img = np.abs(self.orig_img_orders[beam].array) - self.thumb_x = self.orig_img_orders[beam].center.x - self.thumb_y = self.orig_img_orders[beam].center.y - self.img_sh = self.orig_img_orders[beam].array.shape + dx = self.grating_conf.dxlam[beam] xoff = 0 ytrace_beam, lam_beam = self.grating_conf.get_beam_trace(x=self.xcenter, y=self.ycenter, dx=(dx + xoff), @@ -198,8 +169,7 @@ class SpecDisperser(object): dyc = cast[int](np.floor(ytrace_beam+0.5)) - # dypix = cast[int](np.floor(ytrace_beam - dyc[0] + x0[0] + 0.5)) - dypix = dyc - dyc[0] + x0[0] + dypix = cast[int](np.floor(ytrace_beam - dyc[0] + x0[0] + 0.5)) frac_ids = yfrac_beam < 0 @@ -278,8 +248,7 @@ class SpecDisperser(object): # beam_flat[k] = self.flat_cube[:, originOut_y + i, originOut_x + j] status = disperse.disperse_grism_object(self.thumb_img.astype(np.float32), - flat_index[nonz], - yfrac_beam[nonz], + flat_index[nonz], yfrac_beam[nonz], sensitivity_beam[nonz], modelf, x0, array(self.img_sh, @@ -289,24 +258,11 @@ class SpecDisperser(object): lam_beam[lam_index][nonz]) model = modelf.reshape(beam_sh) - # n1 = np.sum(np.isinf(model)) - # n2 = np.sum(np.isnan(model)) - # n3 = np.sum(np.isinf(modelf)) - # n4 = np.sum(np.isnan(modelf)) - # if n1>0 or n2 > 0: - # print("DEBUG: SpecDisperser, inf:%d, nan:%d--------%d,%d"%(n1, n2, n3, n4)) - # print(dypix) - # n1 = np.sum(np.isinf(self.thumb_img.astype(np.float32))) - # n2 = np.sum(np.isnan(self.thumb_img.astype(np.float32))) - # n3 = np.sum(np.isinf(yfrac_beam)) - # n4 = np.sum(np.isnan(yfrac_beam)) - # n5 = np.sum(np.isinf(sensitivity_beam)) - # n6 = np.sum(np.isnan(sensitivity_beam)) - # print("DEBUG: SpecDisperser, innput ---inf:%d, nan:%d, yfrac_beam:%d/%d, sensitivity_beam:%d/%d"%(n1, n2, n3, n4, n5, n6)) self.beam_flux[beam] = sum(modelf) if self.isAlongY == 1: model, _, _ = rotate90(array_orig=model, isClockwise=0) + return model, originOut_x, originOut_y, dxpix, dypix, lam_beam, ysens def writerSensitivityFile(self, conffile='', beam='', w=None, sens=None): diff --git a/observation_sim/mock_objects/SpecDisperser/disperse_c/disperse.pyx b/observation_sim/mock_objects/SpecDisperser/disperse_c/disperse.pyx index 707c670..41dbf14 100644 --- a/observation_sim/mock_objects/SpecDisperser/disperse_c/disperse.pyx +++ b/observation_sim/mock_objects/SpecDisperser/disperse_c/disperse.pyx @@ -20,30 +20,7 @@ import cython cdef extern from "math.h": double sqrt(double x) double exp(double x) - - -def check_nan2D(np.ndarray[FTYPE_t, ndim=2] arr): - cdef int i, j - cdef int nrows = arr.shape[0] - cdef int ncols = arr.shape[1] - # 遍历数组的每个元素并检查是否存在 NaN - for i in range(nrows): - for j in range(ncols): - if np.isnan(arr[i, j]) | np.isinf(arr[i, j]): - return True - return False - -def check_nan1d(np.ndarray[DTYPE_t, ndim=1] arr): - cdef int i - cdef int n = arr.shape[0] - - # 遍历数组的每个元素并检查是否存在 NaN - for i in range(n): - if np.isnan(arr[i]) | np.isinf(arr[i]): - return True - return False - @cython.boundscheck(False) @cython.wraparound(False) @cython.embedsignature(True) @@ -76,18 +53,6 @@ def disperse_grism_object(np.ndarray[FTYPE_t, ndim=2] flam, nk = len(idxl) nl = len(full) - - - #if check_nan2D(flam): - # print("DEBUG: disperse, input Array 'flam' contains NaN.") - #if check_nan1d(ysens): - # print("DEBUG: disperse, input Array 'ysens' contains NaN.") - #if check_nan1d(yfrac): - # print("DEBUG: disperse, input Array 'yfrac' contains NaN.") - #if check_nan1d(full): - # print("DEBUG: disperse, input Array 'full' contains NaN before processing.") - - if (flat is not None): nlamb = len(wlambda) @@ -130,15 +95,14 @@ def disperse_grism_object(np.ndarray[FTYPE_t, ndim=2] flam, else: for i in range(0-x0[1], x0[1]): - x_pos = x0[1]+i - if (x_pos < 0) | (x_pos >= shd[1]): + if (x0[1]+i < 0) | (x0[1]+i >= shd[1]): continue for j in range(0-x0[0], x0[0]): - y_pos = x0[0]+j - if (y_pos < 0) | (y_pos >= shd[0]): + if (x0[0]+j < 0) | (x0[0]+j >= shd[0]): continue - fl_ij = flam[y_pos, x_pos] #/1.e-17 + + fl_ij = flam[x0[0]+j, x0[1]+i] #/1.e-17 if (fl_ij == 0): continue @@ -146,13 +110,10 @@ def disperse_grism_object(np.ndarray[FTYPE_t, ndim=2] flam, k1 = idxl[k]+j*shg[1]+i if (k1 >= 0) & (k1 < nl): full[k1] += ysens[k]*fl_ij*(1-yfrac[k]) + k2 = idxl[k]+(j+1)*shg[1]+i if (k2 >= 0) & (k2 < nl): full[k2] += ysens[k]*fl_ij*yfrac[k] - - #if (check_nan1d(full)): - # print("DEBUG: disperse, output Array 'full' contains NaN after processing.+++++++++++++++++++++++++++") - return True diff --git a/observation_sim/mock_objects/__init__.py b/observation_sim/mock_objects/__init__.py index c6c567f..0a655e1 100755 --- a/observation_sim/mock_objects/__init__.py +++ b/observation_sim/mock_objects/__init__.py @@ -5,4 +5,3 @@ from .Quasar import Quasar from .Star import Star from .Stamp import Stamp from .FlatLED import FlatLED -from .ExtinctionMW import ExtinctionMW diff --git a/observation_sim/psf/PSFGauss.py b/observation_sim/psf/PSFGauss.py index e3036d9..a21b0c8 100644 --- a/observation_sim/psf/PSFGauss.py +++ b/observation_sim/psf/PSFGauss.py @@ -17,7 +17,7 @@ class PSFGauss(PSFModel): self.fwhm = self.fwhmGauss(r=psfRa) self.sigGauss = psfRa # 80% light radius self.sigSpin = sigSpin - self.psf = galsim.Gaussian(flux=1.0, fwhm=self.fwhm) + self.psf = galsim.Gaussian(flux=1.0, fwhm=fwhm) def perfGauss(self, r, sig): """ @@ -122,4 +122,4 @@ class PSFGauss(PSFModel): # return ell, beta, qr PSFshear = galsim.Shear(e=ell, beta=beta*galsim.radians) - return self.psf.shear(PSFshear), PSFshear \ No newline at end of file + return self.psf.shear(PSFshear), PSFshear diff --git a/observation_sim/psf/PSFInterpSLS.py b/observation_sim/psf/PSFInterpSLS.py index 7a68b3b..7adfedc 100644 --- a/observation_sim/psf/PSFInterpSLS.py +++ b/observation_sim/psf/PSFInterpSLS.py @@ -20,10 +20,8 @@ import os from astropy.io import fits from astropy.modeling.models import Gaussian2D -from scipy import signal, interpolate -import datetime -import gc -from jax import numpy as jnp +from scipy import signal + LOG_DEBUG = False # ***# NPSF = 900 # ***# 30*30 @@ -435,16 +433,7 @@ class PSFInterpSLS(PSFModel): # PSF_int_trans[ids_szero] = 0 # print(PSF_int_trans[ids_szero].shape[0],PSF_int_trans.shape) PSF_int_trans = PSF_int_trans/np.sum(PSF_int_trans) - ###DEBGU - ids_szero = PSF_int_trans<0 - n01 = PSF_int_trans[ids_szero].shape[0] - - n1 = np.sum(np.isinf(PSF_int_trans)) - n2 = np.sum(np.isnan(PSF_int_trans)) - if n1>0 or n2>0: - print("DEBUG: PSFInterpSLS, inf:%d, nan:%d, 0 num:%d"%(n1, n2, n01)) - #### # from astropy.io import fits # fits.writeto(str(bandNo) + '_' + g_order+ '_psf_o.fits', PSF_int_trans) @@ -490,215 +479,6 @@ class PSFInterpSLS(PSFModel): return PSF_int_trans, PSF_int - def get_PSF_AND_convolve_withsubImg(self, chip, cutImg=None, pos_img_local=[1000, 1000], bandNo=1, g_order='A', grating_split_pos=3685): - """ - Get the PSF at a given image position - - Parameters: - chip: A 'Chip' object representing the chip we want to extract PSF from. - pos_img: A 'galsim.Position' object representing the image position. - bandpass: A 'galsim.Bandpass' object representing the wavelength range. - pixSize: The pixels size of psf matrix - findNeighMode: 'treeFind' or 'hoclistFind' - Returns: - PSF: A 'galsim.GSObject'. - """ - order_IDs = {'A': '1', 'B': '0', 'C': '0', 'D': '0', 'E': '0'} - contam_order_sigma = {'C': 0.28032344707964174, - 'D': 0.39900182912061344, 'E': 1.1988309797685412} # arcsec - x_start = chip.x_cen/chip.pix_size - chip.npix_x / 2. - y_start = chip.y_cen/chip.pix_size - chip.npix_y / 2. - # print(pos_img.x - x_start) - # pos_img_x = pos_img_local[0] + x_start - # pos_img_y = pos_img_local[1] + y_start - # pos_img = galsim.PositionD(pos_img_x, pos_img_y) - # centerPos_local = cutImg.ncol/2. - if pos_img_local[0] < grating_split_pos: - psf_data = self.grating1_data - else: - psf_data = self.grating2_data - - grating_order = order_IDs[g_order] - # if grating_order in ['-2','-1','2']: - # grating_order = '1' - - # if grating_order in ['0', '1']: - psf_order = psf_data['order'+grating_order] - psf_order_b = psf_order['band'+str(bandNo)] - psf_b_dat = psf_order_b['band_data'] - # pos_p = psf_b_dat[1].data - pos_p = psf_b_dat[1].data/chip.pix_size - np.array([y_start, x_start]) - - pc_coeff = psf_b_dat[2].data - pcs = psf_b_dat[0].data - - npc = 10 - m_size = int(pcs.shape[0]**0.5) - sumImg = np.sum(cutImg.array) - tmp_img = cutImg*0 - for j in np.arange(npc): - X_ = jnp.hstack((pos_p[:,1].flatten()[:, None], pos_p[:,0].flatten()[:, None]),dtype=np.float32) - Z_ = (pc_coeff[j].astype(np.float32)).flatten() - # print(pc_coeff[j].shape[0], pos_p[:,1].shape[0], pos_p[:,0].shape[0]) - cx_len = int(chip.npix_x) - cy_len = int(chip.npix_y) - n_x = jnp.arange(0, cx_len, 1, dtype = int) - n_y = jnp.arange(0, cy_len, 1, dtype = int) - M, N = jnp.meshgrid(n_x, n_y) - # t1=datetime.datetime.now() - # U = interpolate.griddata(X_, Z_, (M[0:cy_len, 0:cx_len],N[0:cy_len, 0:cx_len]), - # method='nearest',fill_value=1.0) - b_img = galsim.Image(cx_len, cy_len) - b_img.setOrigin(0,0) - bounds = cutImg.bounds & b_img.bounds - if bounds.area() == 0: - continue - - # ys = cutImg.ymin - # if ys < 0: - # ys = 0 - # ye = cutImg.ymin+cutImg.nrow - # if ye >= cy_len-1: - # ye = cy_len-1 - # if ye - ys <=0: - # continue - # xs = cutImg.xmin - # if xs < 0: - # xs = 0 - # xe = cutImg.xmin+cutImg.ncol - # if xe >= cx_len-1: - # xe = cx_len-1 - # if xe - xs <=0: - # continue - ys = bounds.ymin - ye = bounds.ymax+1 - xs = bounds.xmin - xe = bounds.xmax+1 - U = interpolate.griddata(X_, Z_, (M[ys:ye, xs:xe],N[ys:ye, xs:xe]), - method='nearest',fill_value=1.0) - # t2=datetime.datetime.now() - - # print("time interpolate:", t2-t1) - - # if U.shape != cutImg.array.shape: - # print('DEBUG:SHAPE',cutImg.ncol,cutImg.nrow,cutImg.xmin, cutImg.ymin) - # continue - img_tmp = cutImg - img_tmp[bounds] = img_tmp[bounds]*U - psf = pcs[:, j].reshape(m_size, m_size) - tmp_img = tmp_img + signal.fftconvolve(img_tmp.array, psf, mode='same', axes=None) - - # t3=datetime.datetime.now() - # print("time convole:", t3-t2) - del U - del img_tmp - if np.sum(tmp_img.array)==0: - tmp_img = cutImg - else: - tmp_img = tmp_img/np.sum(tmp_img.array)*sumImg - return tmp_img - - - def convolveFullImgWithPCAPSF(self, chip, folding_threshold=5.e-3): - keys_L1= chip_utils.getChipSLSGratingID(chip.chipID) - # keys_L2 = ['order-2','order-1','order0','order1','order2'] - keys_L2 = ['order0','order1'] - keys_L3 = ['w1','w2','w3','w4'] - - npca = 10 - - x_start = chip.x_cen/chip.pix_size - chip.npix_x / 2. - y_start = chip.y_cen/chip.pix_size - chip.npix_y / 2. - - for i,gt in enumerate(keys_L1): - psfCo = self.grating1_data - if i > 0: - psfCo = self.grating2_data - for od in keys_L2: - psfCo_L2 = psfCo['order1'] - if od in ['order-2','order-1','order0','order2']: - psfCo_L2 = psfCo['order0'] - for w in keys_L3: - img = chip.img_stack[gt][od][w] - pcs = psfCo_L2['band'+w[1]]['band_data'][0].data - pos_p = psfCo_L2['band'+w[1]]['band_data'][1].data/chip.pix_size - np.array([y_start, x_start]) - pc_coeff = psfCo_L2['band'+w[1]]['band_data'][2].data - # print("DEBUG-----------",np.max(pos_p[:,1]),np.min(pos_p[:,1]), np.max(pos_p[:,0]),np.min(pos_p[:,0])) - sum_img = np.sum(img.array) - - - # coeff_mat = np.zeros([npca, chip.npix_y, chip.npix_x]) - # for m in np.arange(chip.npix_y): - # for n in np.arange(chip.npix_x): - # px = n - # py = m - - # dist2 = (pos_p[:, 1] - px)*(pos_p[:, 1] - px) + (pos_p[:, 0] - py)*(pos_p[:, 0] - py) - # temp_sort_dist = np.zeros([dist2.shape[0], 2]) - # temp_sort_dist[:, 0] = np.arange(0, dist2.shape[0], 1) - # temp_sort_dist[:, 1] = dist2 - # # print(temp_sort_dist) - # dits2_sortlist = sorted(temp_sort_dist, key=lambda x: x[1]) - # # print(dits2_sortlist) - # nearest4p = np.zeros([4, 3]) - # pc_coeff_4p = np.zeros([npca, 4]) - - # for i in np.arange(4): - # smaller_ids = int(dits2_sortlist[i][0]) - # nearest4p[i, 0] = pos_p[smaller_ids, 1] - # nearest4p[i, 1] = pos_p[smaller_ids, 0] - # # print(pos_p[smaller_ids, 1],pos_p[smaller_ids, 0]) - # nearest4p[i, 2] = dits2_sortlist[i][1] - # pc_coeff_4p[:, i] = pc_coeff[npca, smaller_ids] - # # idw_dist = 1/(np.sqrt((px-nearest4p[:, 0]) * (px-nearest4p[:, 0]) + ( - # # py-nearest4p[:, 1]) * (py-nearest4p[:, 1]))) - # idw_dist = 1/(np.sqrt(nearest4p[:, 2])) - - # coeff_int = np.zeros(npca) - # for i in np.arange(4): - # coeff_int = coeff_int + pc_coeff_4p[:, i]*idw_dist[i] - # coeff_mat[:, m, n] = coeff_int - - m_size = int(pcs.shape[0]**0.5) - tmp_img = np.zeros_like(img.array,dtype=np.float32) - for j in np.arange(npca): - print(gt, od, w, j) - X_ = jnp.hstack((pos_p[:,1].flatten()[:, None], pos_p[:,0].flatten()[:, None]),dtype=np.float32) - Z_ = (pc_coeff[j].astype(np.float32)).flatten() - # print(pc_coeff[j].shape[0], pos_p[:,1].shape[0], pos_p[:,0].shape[0]) - sub_size = 4 - cx_len = int(chip.npix_x/sub_size) - cy_len = int(chip.npix_y/sub_size) - n_x = jnp.arange(0, chip.npix_x, sub_size, dtype = int) - n_y = jnp.arange(0, chip.npix_y, sub_size, dtype = int) - - M, N = jnp.meshgrid(n_x, n_y) - t1=datetime.datetime.now() - # U = interpolate.griddata(X_, Z_, (M[0:cy_len, 0:cx_len],N[0:cy_len, 0:cx_len]), - # method='nearest',fill_value=1.0) - U1 = interpolate.griddata(X_, Z_, (M, N), - method='nearest',fill_value=1.0) - U = np.zeros_like(chip.img.array, dtype=np.float32) - for mi in np.arange(cy_len): - for mj in np.arange(cx_len): - U[mi*sub_size:(mi+1)*sub_size, mj*sub_size:(mj+1)*sub_size]=U1[mi,mj] - t2=datetime.datetime.now() - - print("time interpolate:", t2-t1) - - img_tmp = img.array*U - psf = pcs[:, j].reshape(m_size, m_size) - tmp_img = tmp_img + signal.fftconvolve(img_tmp, psf, mode='same', axes=None) - - t3=datetime.datetime.now() - print("time convole:", t3-t2) - del U - del U1 - - chip.img = chip.img + tmp_img*sum_img/np.sum(tmp_img) - del tmp_img - gc.collect() - # pixSize = np.rad2deg(self.pixsize*1e-3/28)*3600 #set psf pixsize # # # assert self.iccd == int(chip.getChipLabel(chipID=chip.chipID)), 'ERROR: self.iccd != chip.chipID' diff --git a/observation_sim/sim_steps/__init__.py b/observation_sim/sim_steps/__init__.py index e8e2661..65a8e21 100644 --- a/observation_sim/sim_steps/__init__.py +++ b/observation_sim/sim_steps/__init__.py @@ -1,13 +1,10 @@ import os - class SimSteps: - def __init__(self, overall_config, chip_output, all_filters, ra_offset=0., dec_offset=0.): + def __init__(self, overall_config, chip_output, all_filters): self.overall_config = overall_config self.chip_output = chip_output self.all_filters = all_filters - self.ra_offset = ra_offset - self.dec_offset = dec_offset from .prepare_headers import prepare_headers, updateHeaderInfo from .add_sky_background import add_sky_background_sci, add_sky_flat_calibration, add_sky_background @@ -18,7 +15,6 @@ class SimSteps: from .readout_output import add_prescan_overscan, add_readout_noise, apply_gain, quantization_and_output from .add_LED_flat import add_LED_Flat - SIM_STEP_TYPES = { "scie_obs": "add_objects", "sky_background": "add_sky_background", @@ -35,6 +31,6 @@ SIM_STEP_TYPES = { "readout_noise": "add_readout_noise", "gain": "apply_gain", "quantization_and_output": "quantization_and_output", - "led_calib_model": "add_LED_Flat", - "sky_flatField": "add_sky_flat_calibration", -} + "led_calib_model":"add_LED_Flat", + "sky_flatField":"add_sky_flat_calibration", +} \ No newline at end of file diff --git a/observation_sim/sim_steps/add_objects.py b/observation_sim/sim_steps/add_objects.py index 6c9c7c1..9f41d36 100644 --- a/observation_sim/sim_steps/add_objects.py +++ b/observation_sim/sim_steps/add_objects.py @@ -145,7 +145,7 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param): # Get position of object on the focal plane pos_img, _, _, _, fd_shear = obj.getPosImg_Offset_WCS( - img=chip.img, fdmodel=fd_model, chip=chip, verbose=False, chip_wcs=chip_wcs, img_header=self.h_ext, ra_offset=self.ra_offset, dec_offset=self.dec_offset) + img=chip.img, fdmodel=fd_model, chip=chip, verbose=False, chip_wcs=chip_wcs, img_header=self.h_ext) # [TODO] For now, only consider objects which their centers (after field distortion) are projected within the focal plane # Otherwise they will be considered missed objects @@ -194,8 +194,7 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param): if isUpdated == 1: # TODO: add up stats - self.chip_output.cat_add_obj( - obj, pos_img, pos_shear, ra_offset=self.ra_offset, dec_offset=self.dec_offset) + self.chip_output.cat_add_obj(obj, pos_img, pos_shear) pass elif isUpdated == 0: missed_obj += 1 @@ -217,27 +216,7 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param): # Unload SED: obj.unload_SED() del obj - # gc.collect() - - if chip.survey_type == "spectroscopic" and not self.overall_config["run_option"]["out_cat_only"] and chip.slsPSFOptim: - # from observation_sim.instruments.chip import chip_utils as chip_utils - # gn = chip_utils.getChipSLSGratingID(chip.chipID)[0] - # img1 = np.zeros([2,chip.img.array.shape[0],chip.img.array.shape[1]]) - - # for id1 in np.arange(2): - # gn = chip_utils.getChipSLSGratingID(chip.chipID)[id1] - # img_i = 0 - # for id2 in ['0','1']: - # o_n = "order"+id2 - # for id3 in ['1','2','3','4']: - # w_n = "w"+id3 - # img1[img_i] = img1[img_i] + chip.img_stack[gn][o_n][w_n].array - # img_i = img_i + 1 - # from astropy.io import fits - # fits.writeto('order0.fits',img1[0],overwrite=True) - # fits.writeto('order1.fits',img1[1],overwrite=True) - - psf_model.convolveFullImgWithPCAPSF(chip) + gc.collect() del psf_model gc.collect() diff --git a/observation_sim/sim_steps/add_pattern_noise.py b/observation_sim/sim_steps/add_pattern_noise.py index f92cf36..fc091c5 100644 --- a/observation_sim/sim_steps/add_pattern_noise.py +++ b/observation_sim/sim_steps/add_pattern_noise.py @@ -27,7 +27,7 @@ def add_poisson_and_dark(self, chip, filt, tel, pointing, catalog, obs_param): InputDark=None) else: chip.img, _ = chip_utils.add_poisson(img=chip.img, - chip=chip, + chip=self, exptime=exptime, poisson_noise=chip.poisson_noise, dark_noise=0.) @@ -81,5 +81,5 @@ def add_bias(self, chip, filt, tel, pointing, catalog, obs_param): nsecx=chip.nsecx, seed=self.overall_config["random_seeds"]["seed_biasNonUniform"]+chip.chipID) elif obs_param["bias_16channel"] == False: - chip.img += chip.bias_level + chip.img += self.bias_level return chip, filt, tel, pointing diff --git a/observation_sim/sim_steps/readout_output.py b/observation_sim/sim_steps/readout_output.py index c2cb8a1..8ef872f 100644 --- a/observation_sim/sim_steps/readout_output.py +++ b/observation_sim/sim_steps/readout_output.py @@ -86,10 +86,10 @@ def quantization_and_output(self, chip, filt, tel, pointing, catalog, obs_param) fname = os.path.join(self.chip_output.subdir, self.h_prim['FILENAME'] + '.fits') - # f_name_size = 68 - # if (len(self.h_prim['FILENAME']) > f_name_size): - # self.updateHeaderInfo(header_flag='prim', keys=['FILENAME'], values=[ - # self.h_prim['FILENAME'][0:f_name_size]]) + f_name_size = 68 + if (len(self.h_prim['FILENAME']) > f_name_size): + self.updateHeaderInfo(header_flag='prim', keys=['FILENAME'], values=[ + self.h_prim['FILENAME'][0:f_name_size]]) hdu1 = fits.PrimaryHDU(header=self.h_prim) diff --git a/submit_jobs.slurm b/submit_jobs.slurm index 69d019f..34867ae 100755 --- a/submit_jobs.slurm +++ b/submit_jobs.slurm @@ -1,10 +1,10 @@ #!/bin/bash #SBATCH -J CSSTSim -#SBATCH -N 2 +#SBATCH -N 1 #SBATCH --ntasks-per-node=6 #SBATCH -p debug -#SBATCH --mem=96G +#SBATCH --mem=60G module load mpi/hpcx/2.4.1/gcc-7.3.1 date @@ -12,4 +12,4 @@ date #限定单节点任务数 srun hostname -s | sort -n | awk -F"-" '{print $2}' | uniq > pnodes -mpirun -mca pml ucx -x UCX_NET_DEVICES=mlx5_0:1 -machinefile pnodes -np 12 --map-by node python3 /public/home/fangyuedong/project/csst_msc_sim/run_sim.py --config_file config_overall.yaml --catalog C9_Catalog -c /public/home/fangyuedong/project/csst_msc_sim/config \ No newline at end of file +mpirun -mca pml ucx -x UCX_NET_DEVICES=mlx5_0:1 -machinefile pnodes -np 6 --map-by node python3 /public/home/fangyuedong/project/csst_msc_sim/run_sim.py --config_file config_overall.yaml --catalog C9_Catalog -c /public/home/fangyuedong/project/csst_msc_sim/config \ No newline at end of file diff --git a/tools/get_pointing.py b/tools/get_pointing.py index 39216a1..073d160 100755 --- a/tools/get_pointing.py +++ b/tools/get_pointing.py @@ -63,15 +63,10 @@ class Chip(object): ycen = self.cen_pix_y if pix_scale == None: pix_scale = self.pix_scale - # dudx = -np.cos(img_rot.rad) * pix_scale - # dudy = -np.sin(img_rot.rad) * pix_scale - # dvdx = -np.sin(img_rot.rad) * pix_scale - # dvdy = +np.cos(img_rot.rad) * pix_scale - - dudx = -np.cos(img_rot.rad) * pix_scale - dudy = +np.sin(img_rot.rad) * pix_scale - dvdx = -np.sin(img_rot.rad) * pix_scale - dvdy = -np.cos(img_rot.rad) * pix_scale + dudx = -np.cos(img_rot.rad) * pix_scale + dudy = -np.sin(img_rot.rad) * pix_scale + dvdx = -np.sin(img_rot.rad) * pix_scale + dvdy = +np.cos(img_rot.rad) * pix_scale # dudx = +np.sin(img_rot.rad) * pix_scale # dudy = +np.cos(img_rot.rad) * pix_scale @@ -144,11 +139,12 @@ def getobsPA(ra, dec): angle = math.acos(np.dot(l1l2cross,pdl2cross)/(np.linalg.norm(l1l2cross)*np.linalg.norm(pdl2cross))) angle = (angle)/math.pi*180 - angle = angle + 90 - if (ra<90 or ra> 270): - angle=-angle + + # if (ra>90 and ra< 270): + # angle=-angle return angle + # @jit() def getSelectPointingList(center = [60,-40], radius = 2): points = np.loadtxt('sky.dat') @@ -265,7 +261,7 @@ def findPointingbyChipID(chipID = 8, ra = 60., dec = -40.): if __name__ == "__main__": - tchip, tra, tdec = 13, 60., -40. + tchip, tra, tdec = 8, 60., -40. pointing = findPointingbyChipID(chipID=tchip, ra=tra, dec=tdec) print("[ra_center, dec_center, image_rot]: ", pointing) diff --git a/tools/get_pointing_accuracy.py b/tools/get_pointing_accuracy.py deleted file mode 100644 index 33d8575..0000000 --- a/tools/get_pointing_accuracy.py +++ /dev/null @@ -1,165 +0,0 @@ - -from pylab import * -import math, sys, numpy as np -import astropy.coordinates as coord -from astropy.coordinates import SkyCoord -from astropy import wcs, units as u -from observation_sim.config.header import ImageHeader -from observation_sim.instruments import Telescope, Chip, FilterParam, Filter, FocalPlane - - -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 ecl2radec(lon_ecl, lat_ecl): - ## convert from ecliptic coordinates to equatorial coordinates - c_ecl = SkyCoord( - lon=lon_ecl * u.degree, lat=lat_ecl * u.degree, frame="barycentrictrueecliptic" - ) - c_eq = c_ecl.transform_to("icrs") - ra_eq, dec_eq = c_eq.ra.degree, c_eq.dec.degree - return ra_eq, dec_eq - - -def radec2ecl(ra, dec): - ## convert from equatorial coordinates to ecliptic coordinates - c_eq = SkyCoord(ra=ra * u.degree, dec=dec * u.degree, frame="icrs") - c_ecl = c_eq.transform_to("barycentrictrueecliptic") - lon_ecl, lat_ecl = c_ecl.lon.degree, c_ecl.lat.degree - return lon_ecl, lat_ecl - -def cal_FoVcenter_1P_equatorial(ra_FieldCenter, dec_FieldCenter, chipID = 1, pa = -23.5): - - ### [ra_FieldCenter, dec_FieldCenter] is the center ra, dec of calibration fileds, such as: NEP, NGC 6397, etc. - ### [ra_ChipCenter, dec_ChipCenter] is the center ra, dec of the Chip center. - ### [ra_PointCenter, dec_PointCenter] is the telescope pointing center. - ## Calculate PA angle - chip = Chip(chipID) - - h_ext = ImageHeader.generateExtensionHeader( - chip=chip, - xlen=chip.npix_x, - ylen=chip.npix_y, - ra=(ra_FieldCenter * u.degree).value, - dec=(dec_FieldCenter * u.degree).value, - pa=pa, - gain=chip.gain, - readout=chip.read_noise, - dark=chip.dark_noise, - saturation=90000, - pixel_scale=chip.pix_scale, - pixel_size=chip.pix_size, - xcen=chip.x_cen, - ycen=chip.y_cen, - ) - - # assume that the telescope point to the sky center; then abtain the chip center under this situation; then calculate the differenc between the assumed telescope center and chip center; then add this difference to the sky center - - wcs_h = wcs.WCS(h_ext) - pixs = np.array([[4608, 4616]]) - world_point = wcs_h.wcs_pix2world(pixs, 0) - ra_ChipCenter, dec_ChipCenter = world_point[0][0], world_point[0][1] - - d_ra = ra_FieldCenter - ra_ChipCenter - d_dec = dec_FieldCenter - dec_ChipCenter - - ra_PointCenter = ra_FieldCenter + d_ra - dec_PointCenter = dec_FieldCenter + d_dec - - lon_ecl_PointCenter, lat_ecl_PointCenter = radec2ecl( - ra_PointCenter, dec_PointCenter - ) - - return ra_PointCenter, dec_PointCenter, lon_ecl_PointCenter, lat_ecl_PointCenter - -def cal_FoVcenter_1P_ecliptic(lon_ecl_FieldCenter, lat_ecl_FieldCenter, chipID = 1, pa = -23.5): - - ### [ra_FieldCenter, dec_FieldCenter] is the center ra, dec of calibration fileds, such as: NEP, NGC 6397, etc. - ### [ra_ChipCenter, dec_ChipCenter] is the center ra, dec of the Chip center. - ### [ra_PointCenter, dec_PointCenter] is the telescope pointing center. - - ra_FieldCenter, dec_FieldCenter = ecl2radec( - lon_ecl_FieldCenter, lat_ecl_FieldCenter - ) - - ## Calculate PA angle - chip = Chip(chipID) - - h_ext = ImageHeader.generateExtensionHeader( - chip=chip, - xlen=chip.npix_x, - ylen=chip.npix_y, - ra=(ra_FieldCenter * u.degree).value, - dec=(dec_FieldCenter * u.degree).value, - pa=pa, - gain=chip.gain, - readout=chip.read_noise, - dark=chip.dark_noise, - saturation=90000, - pixel_scale=chip.pix_scale, - pixel_size=chip.pix_size, - xcen=chip.x_cen, - ycen=chip.y_cen, - ) - - # assume that the telescope point to the sky center; then abtain the chip center under this situation; then calculate the differenc between the assumed telescope center and chip center; then add this difference to the sky center - - wcs_h = wcs.WCS(h_ext) - pixs = np.array([[4608, 4616]]) - world_point = wcs_h.wcs_pix2world(pixs, 0) - ra_ChipCenter, dec_ChipCenter = world_point[0][0], world_point[0][1] - - d_ra = ra_FieldCenter - ra_ChipCenter - d_dec = dec_FieldCenter - dec_ChipCenter - - ra_PointCenter = ra_FieldCenter + d_ra - dec_PointCenter = dec_FieldCenter + d_dec - - lon_ecl_PointCenter, lat_ecl_PointCenter = radec2ecl( - ra_PointCenter, dec_PointCenter - ) - - return ra_PointCenter, dec_PointCenter, lon_ecl_PointCenter, lat_ecl_PointCenter - -def getChipCenterRaDec(chipID = 1, p_ra = 60., p_dec = -40.): - chip = Chip(chipID) - - h_ext = ImageHeader.generateExtensionHeader( - chip=chip, - xlen=chip.npix_x, - ylen=chip.npix_y, - ra=(p_ra * u.degree).value, - dec=(p_dec * u.degree).value, - pa=pa, - gain=chip.gain, - readout=chip.read_noise, - dark=chip.dark_noise, - saturation=90000, - pixel_scale=chip.pix_scale, - pixel_size=chip.pix_size, - xcen=chip.x_cen, - ycen=chip.y_cen, - ) - wcs_h = wcs.WCS(h_ext) - pixs = np.array([[4608, 4616]]) - world_point = wcs_h.wcs_pix2world(pixs, 0) - RA_chip, Dec_chip = world_point[0][0], world_point[0][1] - return RA_chip, Dec_chip - -if __name__ == '__main__': - ra_input, dec_input = 270.00000, 66.56000 # NEP - pa = 23.5 - # chipid = 2 - - for chipid in np.arange(1,31,1): - ra, dec, lon_ecl, lat_ecl = cal_FoVcenter_1P_equatorial( - ra_input, dec_input,chipID=chipid, pa=pa) - - print("chip id is %d, chip center [ra,dec] is [%f, %f], pointing center calculated [ra,dec] is [%f, %f]"%(chipid, ra_input, dec_input, ra, dec)) - #for check the result - # testRA, testDec = getChipCenterRaDec(chipID = chipid, p_ra = ra, p_dec = dec) - # print(ra_input-testRA, dec_input-testDec) - diff --git a/tools/target_location_check.py b/tools/target_location_check.py index 11f22ba..34293be 100644 --- a/tools/target_location_check.py +++ b/tools/target_location_check.py @@ -119,15 +119,10 @@ def getTanWCS(ra, dec, img_rot, pix_scale=0.074): """ xcen, ycen = 0, 0 img_rot = img_rot * galsim.degrees - # dudx = -np.cos(img_rot.rad) * pix_scale - # dudy = -np.sin(img_rot.rad) * pix_scale - # dvdx = -np.sin(img_rot.rad) * pix_scale - # dvdy = +np.cos(img_rot.rad) * pix_scale - dudx = -np.cos(img_rot.rad) * pix_scale - dudy = +np.sin(img_rot.rad) * pix_scale + dudy = -np.sin(img_rot.rad) * pix_scale dvdx = -np.sin(img_rot.rad) * pix_scale - dvdy = -np.cos(img_rot.rad) * pix_scale + dvdy = +np.cos(img_rot.rad) * pix_scale moscen = galsim.PositionD(x=xcen, y=ycen) sky_center = galsim.CelestialCoord( -- GitLab