From 5bfe67607c59cf36e6f999261693801ae515893c Mon Sep 17 00:00:00 2001 From: zhangxin Date: Thu, 11 Jul 2024 10:30:16 +0800 Subject: [PATCH 01/13] add calculate pointing center tools by defining chip center ra dec --- tools/get_pointing_accurcy.py | 165 ++++++++++++++++++++++++++++++++++ 1 file changed, 165 insertions(+) create mode 100644 tools/get_pointing_accurcy.py diff --git a/tools/get_pointing_accurcy.py b/tools/get_pointing_accurcy.py new file mode 100644 index 0000000..33d8575 --- /dev/null +++ b/tools/get_pointing_accurcy.py @@ -0,0 +1,165 @@ + +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) + -- GitLab From f5ae25788dafcf68e140f001796e5017e0f19134 Mon Sep 17 00:00:00 2001 From: fangyuedong Date: Mon, 15 Jul 2024 23:32:37 +0800 Subject: [PATCH 02/13] add radec_offsets (due to aberration compensation by FGS) to wcs projection --- observation_sim/ObservationSim.py | 7 ++++++- observation_sim/mock_objects/MockObject.py | 10 +++++----- observation_sim/sim_steps/__init__.py | 12 ++++++++---- observation_sim/sim_steps/add_objects.py | 2 +- 4 files changed, 20 insertions(+), 11 deletions(-) diff --git a/observation_sim/ObservationSim.py b/observation_sim/ObservationSim.py index 1e82700..247636f 100755 --- a/observation_sim/ObservationSim.py +++ b/observation_sim/ObservationSim.py @@ -92,16 +92,21 @@ class Observation(object): input_time_str=time_str ) ra_cen, dec_cen = ra_cen[0], dec_cen[0] + ra_offset, dec_offset = pointing.ra - ra_cen, pointing.dec - dec_cen else: ra_cen = pointing.ra dec_cen = pointing.dec + ra_offset, dec_offset = 0., 0. # Prepare necessary chip properties for simulation 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) + chip_output=chip_output, + all_filters=self.all_filters, + ra_offset=ra_offset, + dec_offset=dec_offset) for step in pointing.obs_param["call_sequence"]: if self.config["run_option"]["out_cat_only"]: diff --git a/observation_sim/mock_objects/MockObject.py b/observation_sim/mock_objects/MockObject.py index 12e3612..ef660ca 100755 --- a/observation_sim/mock_objects/MockObject.py +++ b/observation_sim/mock_objects/MockObject.py @@ -59,13 +59,13 @@ class MockObject(object): flux = self.getFluxFilter(filt) return flux * tel.pupil_area * exptime - def getPosWorld(self): - ra = self.param["ra"] - dec = self.param["dec"] + def getPosWorld(self, ra_offset=0., dec_offset=0.): + ra = self.param["ra"] + ra_offset + dec = self.param["dec"] + dec_offset 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): - self.posImg = img.wcs.toImage(self.getPosWorld()) + 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)) self.localWCS = img.wcs.local(self.posImg) # Apply field distortion model if (fdmodel is not None) and (chip is not None): diff --git a/observation_sim/sim_steps/__init__.py b/observation_sim/sim_steps/__init__.py index 65a8e21..e8e2661 100644 --- a/observation_sim/sim_steps/__init__.py +++ b/observation_sim/sim_steps/__init__.py @@ -1,10 +1,13 @@ import os + class SimSteps: - def __init__(self, overall_config, chip_output, all_filters): + def __init__(self, overall_config, chip_output, all_filters, ra_offset=0., dec_offset=0.): 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 @@ -15,6 +18,7 @@ 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", @@ -31,6 +35,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", -} \ No newline at end of file + "led_calib_model": "add_LED_Flat", + "sky_flatField": "add_sky_flat_calibration", +} diff --git a/observation_sim/sim_steps/add_objects.py b/observation_sim/sim_steps/add_objects.py index 9f41d36..b50c7af 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) + 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) # [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 -- GitLab From f02a4de08692699023733bf0f486e829e229e839 Mon Sep 17 00:00:00 2001 From: chengliang Date: Tue, 16 Jul 2024 19:40:11 +0800 Subject: [PATCH 03/13] optimize gc.collect --- observation_sim/ObservationSim.py | 3 +++ observation_sim/sim_steps/add_objects.py | 2 +- 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/observation_sim/ObservationSim.py b/observation_sim/ObservationSim.py index 247636f..f80123a 100755 --- a/observation_sim/ObservationSim.py +++ b/observation_sim/ObservationSim.py @@ -135,6 +135,9 @@ 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/sim_steps/add_objects.py b/observation_sim/sim_steps/add_objects.py index b50c7af..741d963 100644 --- a/observation_sim/sim_steps/add_objects.py +++ b/observation_sim/sim_steps/add_objects.py @@ -216,7 +216,7 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param): # Unload SED: obj.unload_SED() del obj - gc.collect() + # gc.collect() del psf_model gc.collect() -- GitLab From 5e8554bd4b6f4a991c13d22fe72f670231614458 Mon Sep 17 00:00:00 2001 From: fangyuedong Date: Mon, 22 Jul 2024 23:11:09 +0800 Subject: [PATCH 04/13] add milky way extinction for galaxy --- catalog/C9_Catalog.py | 19 +++- config/config_overall.yaml | 6 +- config/obs_config_SCI.yaml | 2 +- observation_sim/mock_objects/ExtinctionMW.py | 95 ++++++++++++++++++++ observation_sim/mock_objects/__init__.py | 1 + submit_jobs.slurm | 6 +- 6 files changed, 121 insertions(+), 8 deletions(-) create mode 100644 observation_sim/mock_objects/ExtinctionMW.py diff --git a/catalog/C9_Catalog.py b/catalog/C9_Catalog.py index 948fc96..978d74d 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 +from observation_sim.mock_objects import CatalogBase, Star, Galaxy, Quasar, ExtinctionMW 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,6 +110,12 @@ class Catalog(CatalogBase): self.max_size = 0. + # [TODO] Milky Way extinction + self.mw_ext = ExtinctionMW() + self.mw_ext.init_ext_model(model_name="odonnell") + self.mw_ext.load_Planck_ext( + file_path="/public/home/fangyuedong/project/ext_maps/planck/HFI_CompMap_ThermalDustModel_2048_R1.20.fits") + 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( @@ -248,6 +254,9 @@ class Catalog(CatalogBase): input_time_str=time_str ) + # [TODO] get Milky Way extinction AVs + MW_Av_arr = self.mw_ext.Av_from_Planck(ra=ra_arr, dec=dec_arr) + for igals in range(ngals): # # (TEST) # if igals > 100: @@ -259,6 +268,9 @@ 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 @@ -530,6 +542,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': + 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')) diff --git a/config/config_overall.yaml b/config/config_overall.yaml index 8624cc0..5027e49 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: "no_nonlinearity_test" +run_name: "ext_test" # Project cycle and run counter are used to name the outputs project_cycle: 9 @@ -44,7 +44,7 @@ catalog_options: star_only: NO # Only simulate galaxies? - galaxy_only: NO + galaxy_only: YES # rotate galaxy ellipticity rotateEll: 0. # [degree] @@ -65,7 +65,7 @@ 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: [1] + run_pointings: [0, 1, 2, 3, 4] # Whether to enable astrometric modeling enable_astrometric_model: True diff --git a/config/obs_config_SCI.yaml b/config/obs_config_SCI.yaml index 405950f..ef82bcd 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: [1] +run_chips: [17, 22] # Define observation sequence call_sequence: diff --git a/observation_sim/mock_objects/ExtinctionMW.py b/observation_sim/mock_objects/ExtinctionMW.py new file mode 100644 index 0000000..a98106a --- /dev/null +++ b/observation_sim/mock_objects/ExtinctionMW.py @@ -0,0 +1,95 @@ +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, 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) + 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/__init__.py b/observation_sim/mock_objects/__init__.py index 0a655e1..c6c567f 100755 --- a/observation_sim/mock_objects/__init__.py +++ b/observation_sim/mock_objects/__init__.py @@ -5,3 +5,4 @@ from .Quasar import Quasar from .Star import Star from .Stamp import Stamp from .FlatLED import FlatLED +from .ExtinctionMW import ExtinctionMW diff --git a/submit_jobs.slurm b/submit_jobs.slurm index 34867ae..69d019f 100755 --- a/submit_jobs.slurm +++ b/submit_jobs.slurm @@ -1,10 +1,10 @@ #!/bin/bash #SBATCH -J CSSTSim -#SBATCH -N 1 +#SBATCH -N 2 #SBATCH --ntasks-per-node=6 #SBATCH -p debug -#SBATCH --mem=60G +#SBATCH --mem=96G 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 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 +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 -- GitLab From 1aa4f5339045dae11b99fd09e67aedf9db731741 Mon Sep 17 00:00:00 2001 From: fangyuedong Date: Tue, 23 Jul 2024 08:01:37 +0800 Subject: [PATCH 05/13] 1. add option for galaxy milky way extinction in overall_config 2. bug fix: planck map nested ordering --- catalog/C9_Catalog.py | 21 +++++++++++++------- config/config_overall.yaml | 10 +++++++--- observation_sim/mock_objects/ExtinctionMW.py | 5 +++-- 3 files changed, 24 insertions(+), 12 deletions(-) diff --git a/catalog/C9_Catalog.py b/catalog/C9_Catalog.py index 978d74d..b028bdf 100644 --- a/catalog/C9_Catalog.py +++ b/catalog/C9_Catalog.py @@ -111,10 +111,14 @@ class Catalog(CatalogBase): self.max_size = 0. # [TODO] Milky Way extinction - self.mw_ext = ExtinctionMW() - self.mw_ext.init_ext_model(model_name="odonnell") - self.mw_ext.load_Planck_ext( - file_path="/public/home/fangyuedong/project/ext_maps/planck/HFI_CompMap_ThermalDustModel_2048_R1.20.fits") + 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 @@ -255,7 +259,10 @@ class Catalog(CatalogBase): ) # [TODO] get Milky Way extinction AVs - MW_Av_arr = self.mw_ext.Av_from_Planck(ra=ra_arr, dec=dec_arr) + 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.ones(len(ra_arr)) for igals in range(ngals): # # (TEST) @@ -544,8 +551,8 @@ class Catalog(CatalogBase): y = speci(lamb) # [TODO] Apply Milky Way extinction - if obj.type != 'star': - self.mw_ext.apply_extinction(y, Av=obj.mw_Av) + 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 diff --git a/config/config_overall.yaml b/config/config_overall.yaml index 5027e49..d22dc68 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_test" +run_name: "ext_on" # Project cycle and run counter are used to name the outputs project_cycle: 9 @@ -44,11 +44,15 @@ catalog_options: star_only: NO # Only simulate galaxies? - galaxy_only: YES + galaxy_only: NO # 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 ############################################### @@ -68,7 +72,7 @@ obs_setting: run_pointings: [0, 1, 2, 3, 4] # Whether to enable astrometric modeling - enable_astrometric_model: True + enable_astrometric_model: YES # Cut by saturation magnitude in which band? cut_in_band: "z" diff --git a/observation_sim/mock_objects/ExtinctionMW.py b/observation_sim/mock_objects/ExtinctionMW.py index a98106a..a754deb 100644 --- a/observation_sim/mock_objects/ExtinctionMW.py +++ b/observation_sim/mock_objects/ExtinctionMW.py @@ -13,7 +13,7 @@ class ExtinctionMW(object): @staticmethod def radec2pix(ra, dec, NSIDE=2048): - return hp.pixelfunc.ang2pix(nside=NSIDE, theta=ra, phi=dec, lonlat=True) + return hp.pixelfunc.ang2pix(nside=NSIDE, theta=ra, phi=dec, nested=True, lonlat=True) def init_ext_model(self, model_name="odonnell", Av=1.0, Rv=3.1, lamb=None): self.model_name = model_name @@ -79,7 +79,8 @@ class ExtinctionMW(object): 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) + pix = hp.pixelfunc.ang2pix( + nside=NSIDE, theta=np.pi/2. - b, phi=l, nested=True) return self.ebv_planck[pix] * self.Rv def apply_extinction(self, spec, Av=1.0): -- GitLab From 1ba8e79022cfd21d8ba9ce3c5b771adcdc46ed33 Mon Sep 17 00:00:00 2001 From: zhangxin Date: Wed, 24 Jul 2024 15:10:55 +0800 Subject: [PATCH 06/13] fix file name to get_pointing_accuracy.py --- tools/{get_pointing_accurcy.py => get_pointing_accuracy.py} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename tools/{get_pointing_accurcy.py => get_pointing_accuracy.py} (100%) diff --git a/tools/get_pointing_accurcy.py b/tools/get_pointing_accuracy.py similarity index 100% rename from tools/get_pointing_accurcy.py rename to tools/get_pointing_accuracy.py -- GitLab From b8d1624727f2e38f5eee85db01fc80cb241ff0c9 Mon Sep 17 00:00:00 2001 From: fangyuedong Date: Wed, 24 Jul 2024 21:25:02 +0800 Subject: [PATCH 07/13] fix a typo in add_poisson_and_dark --- observation_sim/sim_steps/add_pattern_noise.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/observation_sim/sim_steps/add_pattern_noise.py b/observation_sim/sim_steps/add_pattern_noise.py index dfe90d0..f92cf36 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=self, + chip=chip, exptime=exptime, poisson_noise=chip.poisson_noise, dark_noise=0.) -- GitLab From e66531a393106988a678adc7bb945599807e83e2 Mon Sep 17 00:00:00 2001 From: fangyuedong Date: Mon, 29 Jul 2024 06:54:42 +0800 Subject: [PATCH 08/13] bug fix C9_Catalog.py and ExtinctionMW.py --- catalog/C9_Catalog.py | 2 +- observation_sim/mock_objects/ExtinctionMW.py | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/catalog/C9_Catalog.py b/catalog/C9_Catalog.py index b028bdf..590dda6 100644 --- a/catalog/C9_Catalog.py +++ b/catalog/C9_Catalog.py @@ -262,7 +262,7 @@ class Catalog(CatalogBase): 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.ones(len(ra_arr)) + MW_Av_arr = np.zeros(len(ra_arr)) for igals in range(ngals): # # (TEST) diff --git a/observation_sim/mock_objects/ExtinctionMW.py b/observation_sim/mock_objects/ExtinctionMW.py index a754deb..85ae692 100644 --- a/observation_sim/mock_objects/ExtinctionMW.py +++ b/observation_sim/mock_objects/ExtinctionMW.py @@ -13,7 +13,7 @@ class ExtinctionMW(object): @staticmethod def radec2pix(ra, dec, NSIDE=2048): - return hp.pixelfunc.ang2pix(nside=NSIDE, theta=ra, phi=dec, nested=True, lonlat=True) + 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 @@ -80,7 +80,7 @@ class ExtinctionMW(object): 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, nested=True) + 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): @@ -91,6 +91,6 @@ class ExtinctionMW(object): raise ValueError( "Need to initialize the extinction model (init_ext_model) first") scale = 10**(-.4 * self.ext * Av) - print("scale = ", scale) + # print("scale = ", scale) spec *= scale return spec -- GitLab From 36106f04b9468e914a43550f49dc7bc2f608ede6 Mon Sep 17 00:00:00 2001 From: fangyuedong Date: Tue, 30 Jul 2024 03:15:54 +0800 Subject: [PATCH 09/13] change bluest and reddest end points (to 2000 and 11000 respectively) for sub bandpasses. --- observation_sim/instruments/data/filters/fgs_sub.list | 2 +- observation_sim/instruments/data/filters/g_sub.list | 4 ++-- observation_sim/instruments/data/filters/i_sub.list | 4 ++-- observation_sim/instruments/data/filters/nuv_sub.list | 4 ++-- observation_sim/instruments/data/filters/r_sub.list | 4 ++-- observation_sim/instruments/data/filters/u_sub.list | 4 ++-- observation_sim/instruments/data/filters/y_sub.list | 4 ++-- observation_sim/instruments/data/filters/z_sub.list | 4 ++-- 8 files changed, 15 insertions(+), 15 deletions(-) diff --git a/observation_sim/instruments/data/filters/fgs_sub.list b/observation_sim/instruments/data/filters/fgs_sub.list index 5a826ee..c02a44f 100644 --- a/observation_sim/instruments/data/filters/fgs_sub.list +++ b/observation_sim/instruments/data/filters/fgs_sub.list @@ -1,4 +1,4 @@ -3000 +2000 4500 4750 5000 diff --git a/observation_sim/instruments/data/filters/g_sub.list b/observation_sim/instruments/data/filters/g_sub.list index f5dda57..1384ade 100644 --- a/observation_sim/instruments/data/filters/g_sub.list +++ b/observation_sim/instruments/data/filters/g_sub.list @@ -1,4 +1,4 @@ -3800 +2000 4217 4432 4631 @@ -6,4 +6,4 @@ 5002 5179 5354 -5799 \ No newline at end of file +11000 \ 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 715f74a..0449403 100644 --- a/observation_sim/instruments/data/filters/i_sub.list +++ b/observation_sim/instruments/data/filters/i_sub.list @@ -1,4 +1,4 @@ -6600 +2000 7061 7255 7448 @@ -6,4 +6,4 @@ 7833 8027 8226 -8999 \ No newline at end of file +11000 \ 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 c89b5dc..b6c60c5 100644 --- a/observation_sim/instruments/data/filters/nuv_sub.list +++ b/observation_sim/instruments/data/filters/nuv_sub.list @@ -1,4 +1,4 @@ -2513 +2000 2621 2716 2805 @@ -6,4 +6,4 @@ 2969 3050 3132 -3499 \ No newline at end of file +11000 \ 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 0abf191..ec0852e 100644 --- a/observation_sim/instruments/data/filters/r_sub.list +++ b/observation_sim/instruments/data/filters/r_sub.list @@ -1,4 +1,4 @@ -5100 +2000 5642 5821 6001 @@ -6,4 +6,4 @@ 6363 6547 6735 -7199 \ No newline at end of file +11000 \ 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 f3a2b72..3ec5a58 100644 --- a/observation_sim/instruments/data/filters/u_sub.list +++ b/observation_sim/instruments/data/filters/u_sub.list @@ -1,4 +1,4 @@ -3000 +2000 3277 3380 3485 @@ -6,4 +6,4 @@ 3703 3813 3918 -4499 \ No newline at end of file +11000 \ 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 81dc243..69a05ee 100644 --- a/observation_sim/instruments/data/filters/y_sub.list +++ b/observation_sim/instruments/data/filters/y_sub.list @@ -1,4 +1,4 @@ -9000 +2000 9322 9405 9489 @@ -6,4 +6,4 @@ 9695 9832 10024 -10590 \ No newline at end of file +11000 \ 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 b0191a5..f105711 100644 --- a/observation_sim/instruments/data/filters/z_sub.list +++ b/observation_sim/instruments/data/filters/z_sub.list @@ -1,4 +1,4 @@ -7800 +2000 8494 8638 8790 @@ -6,4 +6,4 @@ 9141 9363 9663 -10590 \ No newline at end of file +11000 \ No newline at end of file -- GitLab From dfb71c23e1d3a99f280f5393b675daee4903d3ba Mon Sep 17 00:00:00 2001 From: fangyuedong Date: Tue, 30 Jul 2024 20:06:33 +0800 Subject: [PATCH 10/13] magnitude ('mag_use_normal') calculations for galaxies and QSOs now is only based on their SEDs --- catalog/C9_Catalog.py | 18 ++++++++++++++++-- 1 file changed, 16 insertions(+), 2 deletions(-) diff --git a/catalog/C9_Catalog.py b/catalog/C9_Catalog.py index 590dda6..146205d 100644 --- a/catalog/C9_Catalog.py +++ b/catalog/C9_Catalog.py @@ -557,8 +557,22 @@ class Catalog(CatalogBase): # 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 + # 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': 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') -- GitLab From b9723717902e49a0615a6c7f7048da81ae8454bb Mon Sep 17 00:00:00 2001 From: Fang Yuedong Date: Tue, 30 Jul 2024 15:23:01 +0000 Subject: [PATCH 11/13] Add Milky Way extinction and extend the range for SED integration --- catalog/C9_Catalog.py | 39 ++++++++++++++----- config/config_overall.yaml | 10 +++-- .../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 | 7 ++-- .../sim_steps/add_pattern_noise.py | 2 +- 12 files changed, 57 insertions(+), 31 deletions(-) diff --git a/catalog/C9_Catalog.py b/catalog/C9_Catalog.py index 978d74d..146205d 100644 --- a/catalog/C9_Catalog.py +++ b/catalog/C9_Catalog.py @@ -111,10 +111,14 @@ class Catalog(CatalogBase): self.max_size = 0. # [TODO] Milky Way extinction - self.mw_ext = ExtinctionMW() - self.mw_ext.init_ext_model(model_name="odonnell") - self.mw_ext.load_Planck_ext( - file_path="/public/home/fangyuedong/project/ext_maps/planck/HFI_CompMap_ThermalDustModel_2048_R1.20.fits") + 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 @@ -255,7 +259,10 @@ class Catalog(CatalogBase): ) # [TODO] get Milky Way extinction AVs - MW_Av_arr = self.mw_ext.Av_from_Planck(ra=ra_arr, dec=dec_arr) + 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) @@ -544,14 +551,28 @@ class Catalog(CatalogBase): y = speci(lamb) # [TODO] Apply Milky Way extinction - if obj.type != 'star': - self.mw_ext.apply_extinction(y, Av=obj.mw_Av) + 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 + # 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': 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 5027e49..d22dc68 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_test" +run_name: "ext_on" # Project cycle and run counter are used to name the outputs project_cycle: 9 @@ -44,11 +44,15 @@ catalog_options: star_only: NO # Only simulate galaxies? - galaxy_only: YES + galaxy_only: NO # 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 ############################################### @@ -68,7 +72,7 @@ obs_setting: run_pointings: [0, 1, 2, 3, 4] # Whether to enable astrometric modeling - enable_astrometric_model: True + enable_astrometric_model: YES # Cut by saturation magnitude in which band? cut_in_band: "z" diff --git a/observation_sim/instruments/data/filters/fgs_sub.list b/observation_sim/instruments/data/filters/fgs_sub.list index 5a826ee..c02a44f 100644 --- a/observation_sim/instruments/data/filters/fgs_sub.list +++ b/observation_sim/instruments/data/filters/fgs_sub.list @@ -1,4 +1,4 @@ -3000 +2000 4500 4750 5000 diff --git a/observation_sim/instruments/data/filters/g_sub.list b/observation_sim/instruments/data/filters/g_sub.list index f5dda57..1384ade 100644 --- a/observation_sim/instruments/data/filters/g_sub.list +++ b/observation_sim/instruments/data/filters/g_sub.list @@ -1,4 +1,4 @@ -3800 +2000 4217 4432 4631 @@ -6,4 +6,4 @@ 5002 5179 5354 -5799 \ No newline at end of file +11000 \ 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 715f74a..0449403 100644 --- a/observation_sim/instruments/data/filters/i_sub.list +++ b/observation_sim/instruments/data/filters/i_sub.list @@ -1,4 +1,4 @@ -6600 +2000 7061 7255 7448 @@ -6,4 +6,4 @@ 7833 8027 8226 -8999 \ No newline at end of file +11000 \ 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 c89b5dc..b6c60c5 100644 --- a/observation_sim/instruments/data/filters/nuv_sub.list +++ b/observation_sim/instruments/data/filters/nuv_sub.list @@ -1,4 +1,4 @@ -2513 +2000 2621 2716 2805 @@ -6,4 +6,4 @@ 2969 3050 3132 -3499 \ No newline at end of file +11000 \ 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 0abf191..ec0852e 100644 --- a/observation_sim/instruments/data/filters/r_sub.list +++ b/observation_sim/instruments/data/filters/r_sub.list @@ -1,4 +1,4 @@ -5100 +2000 5642 5821 6001 @@ -6,4 +6,4 @@ 6363 6547 6735 -7199 \ No newline at end of file +11000 \ 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 f3a2b72..3ec5a58 100644 --- a/observation_sim/instruments/data/filters/u_sub.list +++ b/observation_sim/instruments/data/filters/u_sub.list @@ -1,4 +1,4 @@ -3000 +2000 3277 3380 3485 @@ -6,4 +6,4 @@ 3703 3813 3918 -4499 \ No newline at end of file +11000 \ 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 81dc243..69a05ee 100644 --- a/observation_sim/instruments/data/filters/y_sub.list +++ b/observation_sim/instruments/data/filters/y_sub.list @@ -1,4 +1,4 @@ -9000 +2000 9322 9405 9489 @@ -6,4 +6,4 @@ 9695 9832 10024 -10590 \ No newline at end of file +11000 \ 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 b0191a5..f105711 100644 --- a/observation_sim/instruments/data/filters/z_sub.list +++ b/observation_sim/instruments/data/filters/z_sub.list @@ -1,4 +1,4 @@ -7800 +2000 8494 8638 8790 @@ -6,4 +6,4 @@ 9141 9363 9663 -10590 \ No newline at end of file +11000 \ No newline at end of file diff --git a/observation_sim/mock_objects/ExtinctionMW.py b/observation_sim/mock_objects/ExtinctionMW.py index a98106a..85ae692 100644 --- a/observation_sim/mock_objects/ExtinctionMW.py +++ b/observation_sim/mock_objects/ExtinctionMW.py @@ -13,7 +13,7 @@ class ExtinctionMW(object): @staticmethod def radec2pix(ra, dec, NSIDE=2048): - return hp.pixelfunc.ang2pix(nside=NSIDE, theta=ra, phi=dec, lonlat=True) + 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 @@ -79,7 +79,8 @@ class ExtinctionMW(object): 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) + 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): @@ -90,6 +91,6 @@ class ExtinctionMW(object): raise ValueError( "Need to initialize the extinction model (init_ext_model) first") scale = 10**(-.4 * self.ext * Av) - print("scale = ", scale) + # print("scale = ", scale) spec *= scale return spec diff --git a/observation_sim/sim_steps/add_pattern_noise.py b/observation_sim/sim_steps/add_pattern_noise.py index dfe90d0..f92cf36 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=self, + chip=chip, exptime=exptime, poisson_noise=chip.poisson_noise, dark_noise=0.) -- GitLab From 615e951fbdd9dafa6fbfb564588d7fda26498381 Mon Sep 17 00:00:00 2001 From: fangyuedong Date: Thu, 1 Aug 2024 09:46:23 +0800 Subject: [PATCH 12/13] bug fix for guidance abberation correction --- observation_sim/ObservationSim.py | 26 +++++++++++++----------- observation_sim/config/ChipOutput.py | 4 ++-- observation_sim/sim_steps/add_objects.py | 5 +++-- 3 files changed, 19 insertions(+), 16 deletions(-) diff --git a/observation_sim/ObservationSim.py b/observation_sim/ObservationSim.py index 4fa42d4..cb060ab 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, slsPSFOptim=False): # Get WCS for the focal plane if wcs_fp == None: wcs_fp = self.focal_plane.getTanWCS( @@ -35,21 +35,22 @@ class Observation(object): 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: + 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']: + for id2 in ['0', '1']: o_n = "order"+id2 allbands = {} - for id3 in ['1','2','3','4']: + 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].setOrigin( + chip.bound.xmin, chip.bound.ymin) allbands[w_n].wcs = wcs_fp - orders[o_n] = allbands + orders[o_n] = allbands chip.img_stack[gn] = orders else: chip.img_stack = {} @@ -111,16 +112,17 @@ class Observation(object): input_date_str=date_str, input_time_str=time_str ) - ra_cen, dec_cen = ra_cen[0], dec_cen[0] - ra_offset, dec_offset = pointing.ra - ra_cen, pointing.dec - dec_cen + ra_offset, dec_offset = pointing.ra - \ + ra_cen[0], pointing.dec - dec_cen[0] else: - ra_cen = pointing.ra - dec_cen = pointing.dec ra_offset, dec_offset = 0., 0. - + 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, slsPSFOptim=slsPSFOpt) # Initialize SimSteps sim_steps = SimSteps(overall_config=self.config, diff --git a/observation_sim/config/ChipOutput.py b/observation_sim/config/ChipOutput.py index 5efd061..7db5e12 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): + def cat_add_obj(self, obj, pos_img, pos_shear, ra_offset=0., dec_offset=0.): 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, obj.dec, obj.ra_orig, obj.dec_orig, obj.z, obj.getMagFilter( + 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( self.filt), obj.type, obj.pmra, obj.pmdec, obj.rv, obj.parallax) line += obj.additional_output_str diff --git a/observation_sim/sim_steps/add_objects.py b/observation_sim/sim_steps/add_objects.py index 9ce6c0a..6c9c7c1 100644 --- a/observation_sim/sim_steps/add_objects.py +++ b/observation_sim/sim_steps/add_objects.py @@ -194,7 +194,8 @@ 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) + self.chip_output.cat_add_obj( + obj, pos_img, pos_shear, ra_offset=self.ra_offset, dec_offset=self.dec_offset) pass elif isUpdated == 0: missed_obj += 1 @@ -236,7 +237,7 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param): # fits.writeto('order0.fits',img1[0],overwrite=True) # fits.writeto('order1.fits',img1[1],overwrite=True) - psf_model.convolveFullImgWithPCAPSF(chip) + psf_model.convolveFullImgWithPCAPSF(chip) del psf_model gc.collect() -- GitLab From edffea7bf814cce27b9322aa0906357480d21822 Mon Sep 17 00:00:00 2001 From: zhangxin Date: Thu, 1 Aug 2024 10:37:55 +0800 Subject: [PATCH 13/13] delete jax lib --- observation_sim/psf/PSFInterpSLS.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/observation_sim/psf/PSFInterpSLS.py b/observation_sim/psf/PSFInterpSLS.py index 7a68b3b..90ac027 100644 --- a/observation_sim/psf/PSFInterpSLS.py +++ b/observation_sim/psf/PSFInterpSLS.py @@ -23,7 +23,7 @@ from astropy.modeling.models import Gaussian2D from scipy import signal, interpolate import datetime import gc -from jax import numpy as jnp +# from jax import numpy as jnp LOG_DEBUG = False # ***# NPSF = 900 # ***# 30*30 @@ -537,14 +537,14 @@ class PSFInterpSLS(PSFModel): 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) + X_ = np.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) + n_x = np.arange(0, cx_len, 1, dtype = int) + n_y = np.arange(0, cy_len, 1, dtype = int) + M, N = np.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) @@ -663,16 +663,16 @@ class PSFInterpSLS(PSFModel): 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) + X_ = np.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) + n_x = np.arange(0, chip.npix_x, sub_size, dtype = int) + n_y = np.arange(0, chip.npix_y, sub_size, dtype = int) - M, N = jnp.meshgrid(n_x, n_y) + M, N = np.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) -- GitLab