diff --git a/Catalog/C3Catalog.py b/Catalog/C3Catalog.py index f28a1930780d2de83a3509a6ed8718557beba647..384599b620e839f2d0cf63f1f4ef74380ad89e33 100644 --- a/Catalog/C3Catalog.py +++ b/Catalog/C3Catalog.py @@ -8,14 +8,16 @@ import astropy.constants as cons from astropy.coordinates import spherical_to_cartesian from astropy.table import Table from scipy import interpolate +from datetime import datetime from ObservationSim.MockObject import CatalogBase, Star, Galaxy, Quasar from ObservationSim.MockObject._util import seds, sed_assign, extAv, tag_sed, getObservedSED +from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position NSIDE = 128 class C3Catalog(CatalogBase): - def __init__(self, config, chip, **kwargs): + def __init__(self, config, chip, pointing, **kwargs): super().__init__() self.cat_dir = os.path.join(config["data_dir"], config["input_path"]["cat_dir"]) self.seed_Av = config["random_seeds"]["seed_Av"] @@ -24,7 +26,9 @@ class C3Catalog(CatalogBase): self.normF_star = Table.read(os.path.join(self.normalize_dir, 'SLOAN_SDSS.g.fits')) self.normF_galaxy = Table.read(os.path.join(self.normalize_dir, 'lsst_throuput_g.fits')) + self.config = config self.chip = chip + self.pointing = pointing if "star_cat" in config["input_path"] and config["input_path"]["star_cat"]: star_file = config["input_path"]["star_cat"] @@ -77,6 +81,11 @@ class C3Catalog(CatalogBase): param = self.initialize_param() param['ra'] = gals['ra_true'][igals] param['dec'] = gals['dec_true'][igals] + if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200): + continue + param['mag_use_normal'] = gals['mag_true_g_lsst'][igals] + if param['mag_use_normal'] >= 26.5: + continue param['z'] = gals['redshift_true'][igals] param['model_tag'] = 'None' param['gamma1'] = 0 @@ -84,10 +93,10 @@ class C3Catalog(CatalogBase): param['kappa'] = 0 param['delta_ra'] = 0 param['delta_dec'] = 0 - sersicB = gals['sersic_bulge'][igals] + # sersicB = gals['sersic_bulge'][igals] hlrMajB = gals['size_bulge_true'][igals] hlrMinB = gals['size_minor_bulge_true'][igals] - sersicD = gals['sersic_disk'][igals] + # sersicD = gals['sersic_disk'][igals] hlrMajD = gals['size_disk_true'][igals] hlrMinD = gals['size_minor_disk_true'][igals] aGal = gals['size_true'][igals] @@ -112,11 +121,6 @@ class C3Catalog(CatalogBase): param['av'] = 0.6 * param['av'] / 3.0 # for quasar, av=[0, 0.2], 3.0=av.max-av.im param['star'] = 2 # Quasar - if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200): - continue - param['mag_use_normal'] = gals['mag_true_g_lsst'][igals] - if param['mag_use_normal'] >= 26.5: - continue self.ids += 1 param['id'] = self.ids @@ -129,10 +133,43 @@ class C3Catalog(CatalogBase): def _load_stars(self, stars, pix_id=None): nstars = len(stars['sourceID']) + # Apply astrometric modeling + ra_arr = stars["RA"][:] + dec_arr = stars["Dec"][:] + if "astrometric_lib" in self.config["obs_setting"] and self.config["obs_setting"]["enable_astrometric_model"]: + ra_list = ra_arr.tolist() + dec_list = dec_arr.tolist() + pmra_list = np.zeros(nstars).tolist() + pmdec_list = np.zeros(nstars).tolist() + rv_list = np.zeros(nstars).tolist() + parallax_list = [1e-9] * nstars + lib_path = os.path.join(self.config["data_dir"], self.config["obs_setting"]["astrometric_lib"]) + dt = datetime.fromtimestamp(self.pointing.timestamp) + date_str = dt.date().isoformat() + time_str = dt.time().isoformat() + ra_arr, dec_arr = on_orbit_obs_position( + input_ra_list=ra_list, + input_dec_list=dec_list, + input_pmra_list=pmra_list, + input_pmdec_list=pmdec_list, + input_rv_list=rv_list, + input_parallax_list=parallax_list, + input_nstars=nstars, + input_x=self.pointing.sat_x, + input_y=self.pointing.sat_y, + input_z=self.pointing.sat_z, + input_vx=self.pointing.sat_vx, + input_vy=self.pointing.sat_vy, + input_vz=self.pointing.sat_vz, + input_epoch="J2015.5", + input_date_str=date_str, + input_time_str=time_str, + lib_path=lib_path + ) for istars in range(nstars): param = self.initialize_param() - param['ra'] = stars['RA'][istars] - param['dec'] = stars['Dec'][istars] + param['ra'] = ra_arr[istars] + param['dec'] = dec_arr[istars] if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200): continue param['mag_use_normal'] = stars['app_sdss_g'][istars] diff --git a/Catalog/Catalog_example.py b/Catalog/Catalog_example.py index 3c58780d5da8f5b93d65bfbb91ccde30a7251180..b27a0393f65ae9df6d61234151c06d2188e18a8b 100644 --- a/Catalog/Catalog_example.py +++ b/Catalog/Catalog_example.py @@ -33,7 +33,7 @@ class Catalog_example(CatalogBase): load the filter throughput for the input catalog's photometric system. """ - def __init__(self, config, chip, **kwargs): + def __init__(self, config, chip, pointing, **kwargs): """Constructor method. Parameters @@ -54,6 +54,7 @@ class Catalog_example(CatalogBase): super().__init__() self.cat_dir = os.path.join(config["data_dir"], config["input_path"]["cat_dir"]) self.chip = chip + self.pointing = pointing if "star_cat" in config["input_path"] and config["input_path"]["star_cat"]: star_file = config["input_path"]["star_cat"] star_SED_file = config["SED_templates_path"]["star_SED"] diff --git a/ObservationSim/Astrometry/Astrometry_util.py b/ObservationSim/Astrometry/Astrometry_util.py new file mode 100644 index 0000000000000000000000000000000000000000..4667a17d695dbd38ddb7ff62faaaf4b81a54b20d --- /dev/null +++ b/ObservationSim/Astrometry/Astrometry_util.py @@ -0,0 +1,124 @@ +from ctypes import * +import numpy as np + +def checkInputList(input_list, n): + if not isinstance(input_list, list): + raise TypeError("Input type is not list!", input_list) + for i in input_list: + if type(i) != type(1.1): + if type(i) != type(1): + raise TypeError("Input list's element is not float or int!", input_list) + if len(input_list) != n: + raise RuntimeError("Length of input list is not equal to stars' number!", input_list) + +def on_orbit_obs_position(input_ra_list, input_dec_list, input_pmra_list, input_pmdec_list, input_rv_list, input_parallax_list, input_nstars, input_x, input_y, input_z, input_vx, input_vy, input_vz, input_epoch, input_date_str, input_time_str, lib_path): + #Check input parameters + if not isinstance(input_nstars, int): + raise TypeError("Parameter 7 is not int!", input_nstars) + + checkInputList(input_ra_list, input_nstars) + checkInputList(input_dec_list, input_nstars) + checkInputList(input_pmra_list, input_nstars) + checkInputList(input_pmdec_list, input_nstars) + checkInputList(input_rv_list, input_nstars) + checkInputList(input_parallax_list, input_nstars) + + if not isinstance(input_x, float): + raise TypeError("Parameter 8 is not double!", input_x) + if not isinstance(input_y, float): + raise TypeError("Parameter 9 is not double!", input_y) + if not isinstance(input_z, float): + raise TypeError("Parameter 10 is not double!", input_z) + if not isinstance(input_vx, float): + raise TypeError("Parameter 11 is not double!", input_vx) + if not isinstance(input_vy, float): + raise TypeError("Parameter 12 is not double!", input_vy) + if not isinstance(input_vz, float): + raise TypeError("Parameter 13 is not double!", input_vz) + #Convert km -> m + input_x = input_x*1000.0 + input_y = input_y*1000.0 + input_z = input_z*1000.0 + input_vx = input_vx*1000.0 + input_vy = input_vy*1000.0 + input_vz = input_vz*1000.0 + + if not isinstance(input_date_str, str): + raise TypeError("Parameter 15 is not string!", input_date_str) + else: + input_date_str = input_date_str.strip() + if not (input_date_str[4]=="-" and input_date_str[7]=="-"): + raise TypeError("Parameter 15 format error (1)!", input_date_str) + else: + tmp = input_date_str.split("-") + if len(tmp) != 3: + raise TypeError("Parameter 15 format error (2)!", input_date_str) + input_year = int(tmp[0]) + input_month = int(tmp[1]) + input_day = int(tmp[2]) + if not (input_year>=1900 and input_year<=2100): + raise TypeError("Parameter 15 year range error [1900 ~ 2100]!", input_year) + if not (input_month>=1 and input_month<=12): + raise TypeError("Parameter 15 month range error [1 ~ 12]!", input_month) + if not (input_day>=1 and input_day<=31): + raise TypeError("Parameter 15 day range error [1 ~ 31]!", input_day) + + if not isinstance(input_time_str, str): + raise TypeError("Parameter 16 is not string!", input_time_str) + else: + input_time_str = input_time_str.strip() + if not (input_time_str[2]==":" and input_time_str[5]==":"): + raise TypeError("Parameter 16 format error (1)!", input_time_str) + else: + tmp = input_time_str.split(":") + if len(tmp) != 3: + raise TypeError("Parameter 16 format error (2)!", input_time_str) + input_hour = int(tmp[0]) + input_minute = int(tmp[1]) + input_second = float(tmp[2]) + if not (input_hour>=0 and input_hour<=23): + raise TypeError("Parameter 16 hour range error [0 ~ 23]!", input_hour) + if not (input_minute>=0 and input_minute<=59): + raise TypeError("Parameter 16 minute range error [0 ~ 59]!", input_minute) + if not (input_second>=0 and input_second<60.0): + raise TypeError("Parameter 16 second range error [0 ~ 60)!", input_second) + #Inital dynamic lib + shao = cdll.LoadLibrary(lib_path) + shao.onOrbitObs.restype = c_int + + d3 = c_double * 3 + shao.onOrbitObs.argtypes = [c_double, c_double, c_double, c_double, c_double, c_double, \ + c_int, c_int, c_int, c_int, c_int, c_double, \ + c_double, POINTER(d3), POINTER(d3), \ + c_int, c_int, c_int, c_int, c_int, c_double, \ + POINTER(c_double), POINTER(c_double) ] + output_ra_list = list() + output_dec_list = list() + for i in range(input_nstars): + input_ra = c_double(input_ra_list[i]) + input_dec = c_double(input_dec_list[i]) + input_pmra = c_double(input_pmra_list[i]) + input_pmdec = c_double(input_pmdec_list[i]) + input_rv = c_double(input_rv_list[i]) + input_parallax = c_double(input_parallax_list[i]) + p3 = d3(input_x, input_y, input_z) + v3 = d3(input_vx, input_vy, input_vz) + input_year_c = c_int(input_year) + input_month_c = c_int(input_month) + input_day_c = c_int(input_day) + input_hour_c = c_int(input_hour) + input_minute_c = c_int(input_minute) + input_second_c = c_double(input_second) + DAT = c_double(37.0) + output_ra = c_double(0.0) + output_dec = c_double(0.0) + rs = shao.onOrbitObs(input_ra, input_dec, input_pmra, input_pmdec, input_rv, input_parallax, \ + input_year_c, input_month_c, input_day_c, input_hour_c, input_minute_c, input_second_c, \ + DAT, byref(p3), byref(v3), \ + input_year_c, input_month_c, input_day_c, input_hour_c, input_minute_c, input_second_c, \ + byref(output_ra), byref(output_dec)) + if rs != 0: + raise RuntimeError("Calculate error!") + output_ra_list.append(output_ra.value) + output_dec_list.append(output_dec.value) + return np.array(output_ra_list), np.array(output_dec_list) diff --git a/ObservationSim/Astrometry/__init__.py b/ObservationSim/Astrometry/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/ObservationSim/Config/Config.py b/ObservationSim/Config/Config.py index 0f6dcf891f03d033f6cf3e462d46e64b97487cf3..b9478b3495b913924a206ee097b817a5991880ec 100755 --- a/ObservationSim/Config/Config.py +++ b/ObservationSim/Config/Config.py @@ -2,7 +2,7 @@ import galsim import os from astropy.time import Time as asTime -def ConfigDir(config, work_dir=None, data_dir=None): +def config_dir(config, work_dir=None, data_dir=None): path_dict = {} # Working directory if work_dir == None: @@ -19,7 +19,6 @@ def ConfigDir(config, work_dir=None, data_dir=None): path_dict["data_dir"] = data_dir # Data sub-catalogs # Object catalog direcotry - # path_dict["cat_dir"] = os.path.join(path_dict["data_dir"], "catalog_points_7degree2/", cat_dir) path_dict["cat_dir"] = os.path.join(path_dict["data_dir"], config["input_path"]["cat_dir"]) # PSF data directory path_dict["psf_dir"] = os.path.join(path_dict["data_dir"], config["psf_setting"]["psf_dir"]) @@ -27,8 +26,6 @@ def ConfigDir(config, work_dir=None, data_dir=None): # SED catalog directory # TODO: SED_dir is deprecated path_dict["SED_dir"] = os.path.join(path_dict["data_dir"], "imageSims/Catalog/SEDObject") - # path_dict["template_dir"] = path_dict["data_dir"] + "Templates/" - # path_dict["template_dir"] = os.path.join(path_dict["data_dir"], config["SED_templates_path"]["galaxy_SED"]) # Directories/files for instrument parameters, e.g. efficiency curves. path_dict["filter_dir"] = os.path.join(path_dict["data_dir"], config["Efficiency_curve_path"]["filter_eff"]) path_dict["ccd_dir"] = os.path.join(path_dict["data_dir"], config["Efficiency_curve_path"]["ccd_eff"]) @@ -44,7 +41,7 @@ def ConfigDir(config, work_dir=None, data_dir=None): return path_dict -def ReadConfig(config_filename): +def read_config(config_filename): """Read in a configuration file and return the corresponding dict(s). Parameters: @@ -73,10 +70,10 @@ def ReadConfig(config_filename): else: print("!! Something is wrong with parameter '%s'."%row[0]) return - config = ParseConfig(config) + config = parse_config(config) return config -def ParseConfig(config): +def parse_config(config): """Parse the config values to the right type Parameters: diff --git a/ObservationSim/Config/Pointing.py b/ObservationSim/Config/Pointing.py new file mode 100644 index 0000000000000000000000000000000000000000..3d6557e0db0c9bd84721af987e7faf51d209e233 --- /dev/null +++ b/ObservationSim/Config/Pointing.py @@ -0,0 +1,36 @@ +import numpy as np +import galsim +from astropy.time import Time + +class Pointing(object): + def __init__(self, id=0, ra=0., dec=0., img_pa=0., timestamp=1621915200, sat_x=0., sat_y=0, sat_z=0, sat_vx=0, sat_vy=0, sat_vz=0, exp_time=150., pointing_type='MS'): + self.id = id + self.ra = ra + self.dec = dec + self.img_pa = img_pa * galsim.degrees + self.timestamp = timestamp + self.sat_x, self.sat_y, self.sat_z = sat_x, sat_y, sat_z + self.sat_vx, self.sat_vy, self.sat_vz = sat_vx, sat_vy, sat_vz + self.exp_time = exp_time + self.pointing_type = pointing_type + + def read_pointing_columns(self, columns, id=0, t=1621915200, pointing_type='MS'): + self.id = id + col_len = len(columns) + self.ra = float(columns[0]) + self.dec = float(columns[1]) + self.img_pa = float(columns[4]) * galsim.degrees + self.pointing_type = pointing_type + if col_len > 5: + jdt = np.double(columns[5]) + t_temp = Time(jdt, format='jd') + self.timestamp = t_temp.unix + self.sat_x = float(columns[6]) + self.sat_y = float(columns[7]) + self.sat_z = float(columns[8]) + self.sat_vx = float(columns[15]) + self.sat_vy = float(columns[16]) + self.sat_vz = float(columns[17]) + self.exp_time = float(columns[18]) + else: + self.timestamp = t diff --git a/ObservationSim/Config/__init__.py b/ObservationSim/Config/__init__.py index 0d8b52ef73e40d3132de0c8548b93211b530a243..96ca93e25999a633ba3697782b2d72e7e66b17fa 100644 --- a/ObservationSim/Config/__init__.py +++ b/ObservationSim/Config/__init__.py @@ -1,2 +1,3 @@ -from .Config import ConfigDir, ReadConfig, ParseConfig -from .ChipOutput import ChipOutput \ No newline at end of file +from .Config import * +from .ChipOutput import ChipOutput +from .Pointing import Pointing \ No newline at end of file diff --git a/ObservationSim/Instrument/FilterParam.py b/ObservationSim/Instrument/FilterParam.py index c7a8e876270aea6780987e29403dadce60b6af2e..0437ed545acbc5a8cc489730dcfcda46e397e37b 100755 --- a/ObservationSim/Instrument/FilterParam.py +++ b/ObservationSim/Instrument/FilterParam.py @@ -1,5 +1,4 @@ import galsim -import pylab as pl import numpy as np class FilterParam(object): diff --git a/ObservationSim/MockObject/CatalogBase.py b/ObservationSim/MockObject/CatalogBase.py index bcb287e4a19232b98689090438c4d4dcd980e205..08e82ca110ad837bb105849cced1dbc73f918b96 100644 --- a/ObservationSim/MockObject/CatalogBase.py +++ b/ObservationSim/MockObject/CatalogBase.py @@ -46,7 +46,13 @@ class CatalogBase(metaclass=ABCMeta): "ell_tot":0, "teff":0, "logg":0, - "feh":0 + "feh":0, + "g1":0, + "g2":0, + "pmra":0, + "pmdec":0, + "rv":0, + "parallax":1e-9 } return param diff --git a/ObservationSim/MockObject/Galaxy.py b/ObservationSim/MockObject/Galaxy.py index 8462114c424aeca4853617995b4657071f3115c1..fb0d4c3b68d54a75960bf53f3a4ccb7a4ee51354 100755 --- a/ObservationSim/MockObject/Galaxy.py +++ b/ObservationSim/MockObject/Galaxy.py @@ -28,90 +28,11 @@ class Galaxy(MockObject): if rotation is not None: self.rotateEllipticity(rotation) - # def load_SED(self, survey_type, sed_path=None, cosids=None, objtypes=None, sed_templates=None, normFilter=None, target_filt=None): - # if survey_type == "photometric": - # norm_thr_rang_ids = normFilter['SENSITIVITY'] > 0.001 - # if sed_templates is None: - # # Read SED data directly - # itype = objtypes[cosids==self.sed_type][0] - # sed_file = os.path.join(sed_path, itype + "_ID%s.sed"%(self.sed_type)) - # if not os.path.exists(sed_file): - # raise ValueError("!!! No SED found.") - # sed_data = Table.read(sed_file, format="ascii") - # wave, flux = sed_data["observedLambda"].data, sed_data["observedFlux"].data - # else: - # # Load SED from templates - # sed_data = sed_templates[self.sed_type] - # # redshift, intrinsic extinction - # sed_data = getObservedSED( - # sedCat=sed_data, - # redshift=self.z, - # av=self.param['av'], - # redden=self.param['redden']) - # wave, flux = sed_data[0], sed_data[1] - # flux_photon = flux * (wave / (cons.h.value * cons.c.value)) * 1e-13 - # sed_photon = Table(np.array([wave, flux_photon]).T, names=('WAVELENGTH', 'FLUX')) - # # Get scaling factor for SED - # sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=self.param['mag_use_normal'], - # spectrum=sed_photon, - # norm_thr=normFilter, - # sWave=np.floor(normFilter[norm_thr_rang_ids][0][0]), - # eWave=np.ceil(normFilter[norm_thr_rang_ids][-1][0])) - # sed_photon = np.array([sed_photon['WAVELENGTH'], sed_photon['FLUX']*sedNormFactor]).T - # # Convert to galsim.SED object - # spec = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(sed_photon[:, 1]), interpolant='nearest') - # self.sed = galsim.SED(spec, wave_type='A', flux_type='1', fast=False) - # # Get magnitude - # interFlux = integrate_sed_bandpass(sed=self.sed, bandpass=target_filt.bandpass_full) - # self.param['mag_%s'%target_filt.filter_type] = getABMAG( - # interFlux=interFlux, - # bandpass=target_filt.bandpass_full) - # # print('mag_use_normal = ', self.param['mag_use_normal']) - # # print('mag_%s = '%target_filt.filter_type, self.param['mag_%s'%target_filt.filter_type]) - # # print('redshift = %.3f'%(self.z)) - # # print('sed_type = %d, av = %.2f, redden = %d'%(self.sed_type, self.param['av'], self.param['redden'])) - - # elif survey_type == "spectroscopic": - # if sed_templates is None: - # self.sedPhotons(sed_path=sed_path, cosids=cosids, objtypes=objtypes) - # else: - # sed_data = sed_templates[self.sed_type] - # sed_data = getObservedSED( - # sedCat=sed_data, - # redshift=self.z, - # av=self.param['av'], - # redden=self.param['redden']) - # speci = interpolate.interp1d(sed_data[0], sed_data[1]) - # lamb = np.arange(2500, 10001 + 0.5, 0.5) - # y = speci(lamb) - # # erg/s/cm2/A --> photo/s/m2/A - # all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13 - # self.sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX')) - - def unload_SED(self): """(Test) free up SED memory """ del self.sed - # def sedPhotons(self, sed_path, cosids, objtypes): - - # itype = objtypes[cosids == self.sed_type][0] - # sed_file = os.path.join(sed_path, itype + "_ID%s.sed" % (self.sed_type)) - # if not os.path.exists(sed_file): - # raise ValueError("!!! No SED found.") - # sed = Table.read(sed_file, format="ascii") - # spec_data = {} - # f_orig = sed["observedFlux"].data - # w_orig = sed["observedLambda"].data - - # speci = interpolate.interp1d(w_orig, f_orig) - # lamb = np.arange(2500, 10001 + 0.5, 0.5) - # y = speci(lamb) - # # erg/s/cm2/A --> photo/s/m2/A - # all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13 - # self.sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX')) - def getGSObj_multiband(self, tel, psf_list, bandpass_list, filt, nphotons_tot=None, g1=0, g2=0, exptime=150.): if len(psf_list) != len(bandpass_list): raise ValueError("!!!The number of PSF profiles and the number of bandpasses must be equal.") diff --git a/ObservationSim/MockObject/Star.py b/ObservationSim/MockObject/Star.py index 31cba083ddcf255eff3d323248861f1a7d26be2c..c80334a0cd329edfcdebc0fa3334262d9766ea9d 100755 --- a/ObservationSim/MockObject/Star.py +++ b/ObservationSim/MockObject/Star.py @@ -17,61 +17,6 @@ class Star(MockObject): """ del self.sed - # def load_SED(self, survey_type, normFilter=None, target_filt=None, sed_lib=None, sed_path=None): - # if survey_type == "photometric": - # norm_thr_rang_ids = normFilter['SENSITIVITY'] > 0.001 - # # spec_lambda = Table.read(sed_path, path=f"/SED/wave_{self.model_tag}") - # # spec_flux = Table.read(sed_path, path=f"/SED/{self.sed_type}") - # # wave, flux = spec_lambda['col0'].data, spec_flux['col0'].data - # _, wave, flux = tag_sed( - # h5file=sed_lib, - # model_tag=self.param['model_tag'], - # teff=self.param['teff'], - # logg=self.param['logg'], - # feh=self.param['feh']) - # flux_photon = flux * (wave / (cons.h.value * cons.c.value)) * 1e-13 - # sed_photon = Table(np.array([wave, flux_photon]).T, names=('WAVELENGTH', 'FLUX')) - # # Get scaling factor for SED - # sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=self.param['mag_use_normal'], - # spectrum=sed_photon, - # norm_thr=normFilter, - # sWave=np.floor(normFilter[norm_thr_rang_ids][0][0]), - # eWave=np.ceil(normFilter[norm_thr_rang_ids][-1][0])) - # sed_photon = np.array([sed_photon['WAVELENGTH'], sed_photon['FLUX']*sedNormFactor]).T - # # Convert to galsim.SED object - # spec = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(sed_photon[:, 1]), interpolant='nearest') - # self.sed = galsim.SED(spec, wave_type='A', flux_type='1', fast=False) - # # Get magnitude - # interFlux = integrate_sed_bandpass(sed=self.sed, bandpass=target_filt.bandpass_full) - # self.param['mag_%s'%target_filt.filter_type] = getABMAG( - # interFlux=interFlux, - # bandpass=target_filt.bandpass_full) - # # print('mag_use_normal = ', self.param['mag_use_normal']) - # # print('mag_%s = '%target_filt.filter_type, self.param['mag_%s'%target_filt.filter_type]) - - # elif survey_type == "spectroscopic": - # # self.sedPhotons(sed_path=sed_path) - # self.sedPhotons(sed_lib=sed_lib) - - # def sedPhotons(self, sed_path=None, sed_lib=None): - # # spec_lambda = Table.read(sed_path, path=f"/SED/wave_{self.model_tag}") - # # spec_flux = Table.read(sed_path, path=f"/SED/{self.sed_type}") - # _, w_orig, f_orig = tag_sed( - # h5file=sed_lib, - # model_tag=self.param['model_tag'], - # teff=self.param['teff'], - # logg=self.param['logg'], - # feh=self.param['feh']) - # # spec_data = {} - # # f_orig = spec_flux["col0"].data - # # w_orig = spec_lambda["col0"].data - # speci = interpolate.interp1d(w_orig, f_orig) - # lamb = np.arange(2500, 10001 + 0.5, 0.5) - # y = speci(lamb) - # # erg/s/cm2/A --> photo/s/m2/A - # all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13 - # self.sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX')) - def getGSObj(self, psf, g1=0, g2=0, flux=None, filt=None, tel=None, exptime=150.): if flux == None: flux = self.getElectronFluxFilt(filt, tel, exptime) @@ -96,6 +41,7 @@ class Star(MockObject): for i in range(len(bandpass_list)): bandpass = bandpass_list[i] + psf = psf_list[i] try: sub = integrate_sed_bandpass(sed=self.sed, bandpass=bandpass) except Exception as e: diff --git a/ObservationSim/ObservationSim.py b/ObservationSim/ObservationSim.py index ece5dd7bee91451bc478b502d92c9aca9313be42..85b667e227c3ba829ff8931e977a67d636adf026 100755 --- a/ObservationSim/ObservationSim.py +++ b/ObservationSim/ObservationSim.py @@ -1,29 +1,26 @@ -import sys import os - -from ObservationSim.Config import ConfigDir, ChipOutput -from ObservationSim.Config.Header import generatePrimaryHeader, generateExtensionHeader -from ObservationSim.Instrument import Telescope, Filter, FilterParam, FocalPlane, Chip -from ObservationSim.MockObject import calculateSkyMap_split_g -from ObservationSim.PSF import PSFGauss, FieldDistortion, PSFInterp -from ObservationSim._util import getShearFiled, makeSubDir_PointingList - -from astropy.io import fits -from datetime import datetime import numpy as np import mpi4py.MPI as MPI import galsim import logging import psutil +from astropy.io import fits +from datetime import datetime + +from ObservationSim.Config import config_dir, ChipOutput +from ObservationSim.Config.Header import generatePrimaryHeader, generateExtensionHeader +from ObservationSim.Instrument import Telescope, Filter, FilterParam, FocalPlane, Chip +from ObservationSim.MockObject import calculateSkyMap_split_g +from ObservationSim.PSF import PSFGauss, FieldDistortion, PSFInterp +from ObservationSim._util import get_shear_field, makeSubDir_PointingList class Observation(object): def __init__(self, config, Catalog, work_dir=None, data_dir=None): - self.path_dict = ConfigDir(config=config, work_dir=work_dir, data_dir=data_dir) - # self.config = ReadConfig(self.path_dict["config_file"]) + self.path_dict = config_dir(config=config, work_dir=work_dir, data_dir=data_dir) self.config = config - self.tel = Telescope(optEffCurve_path=self.path_dict["mirror_file"]) # Currently the default values are hard coded in - self.focal_plane = FocalPlane(survey_type=self.config["obs_setting"]["survey_type"]) # Currently the default values are hard coded in - self.filter_param = FilterParam(filter_dir=self.path_dict["filter_dir"]) # Currently the default values are hard coded in + self.tel = Telescope(optEffCurve_path=self.path_dict["mirror_file"]) + self.focal_plane = FocalPlane(survey_type=self.config["obs_setting"]["survey_type"]) + self.filter_param = FilterParam(filter_dir=self.path_dict["filter_dir"]) self.chip_list = [] self.filter_list = [] self.Catalog = Catalog @@ -42,23 +39,34 @@ class Observation(object): continue # Make Chip & Filter lists - chip = Chip(chipID, ccdEffCurve_dir=self.path_dict["ccd_dir"], CRdata_dir=self.path_dict["CRdata_dir"], normalize_dir=self.path_dict["normalize_dir"], sls_dir=self.path_dict["sls_dir"], config=self.config) # currently there is no config file for chips + chip = Chip( + chipID=chipID, + ccdEffCurve_dir=self.path_dict["ccd_dir"], + CRdata_dir=self.path_dict["CRdata_dir"], + normalize_dir=self.path_dict["normalize_dir"], + sls_dir=self.path_dict["sls_dir"], config=self.config) filter_id, filter_type = chip.getChipFilter() - filt = Filter(filter_id=filter_id, filter_type=filter_type, filter_param=self.filter_param, ccd_bandpass=chip.effCurve) + filt = Filter(filter_id=filter_id, + filter_type=filter_type, + filter_param=self.filter_param, + ccd_bandpass=chip.effCurve) self.chip_list.append(chip) self.filter_list.append(filt) # Read catalog and shear(s) - self.g1_field, self.g2_field, self.nshear = getShearFiled(config=self.config) + self.g1_field, self.g2_field, self.nshear = get_shear_field(config=self.config) + def run_one_chip(self, chip, filt, pointing, chip_output, wcs_fp=None, psf_model=None, shear_cat_file=None, cat_dir=None, sed_dir=None): - def runOneChip(self, chip, filt, chip_output, wcs_fp=None, psf_model=None, pointing_ID=0, ra_cen=None, dec_cen=None, img_rot=None, exptime=150., timestamp_obs=1621915200, pointing_type="MS", shear_cat_file=None, cat_dir=None, sed_dir=None): - - if (ra_cen is None) or (dec_cen is None): - ra_cen = self.config["obs_setting"]["ra_center"] - dec_cen = self.config["obs_setting"]["dec_center"] - if img_rot is None: - img_rot = self.config["obs_setting"]["image_rot"] + print(':::::::::::::::::::Current Pointing Information::::::::::::::::::') + print("RA: %f, DEC; %f" % (pointing.ra, pointing.dec)) + print("Time: %s" % datetime.fromtimestamp(pointing.timestamp).isoformat()) + print("Exposure time: %f" % pointing.exp_time) + print("Satellite Position (x, y, z): (%f, %f, %f)" % (pointing.sat_x, pointing.sat_y, pointing.sat_z)) + print("Satellite Velocity (x, y, z): (%f, %f, %f)" % (pointing.sat_vx, pointing.sat_vy, pointing.sat_vz)) + print("Position Angle: %f" % pointing.img_pa.deg) + print('Chip : %d' % chip.chipID) + print(':::::::::::::::::::::::::::END:::::::::::::::::::::::::::::::::::') if self.config["psf_setting"]["psf_model"] == "Gauss": psf_model = PSFGauss(chip=chip) @@ -69,11 +77,11 @@ class Observation(object): # Get (extra) shear fields if shear_cat_file is not None: - self.g1_field, self.g2_field, self.nshear = getShearFiled(config=self.config, shear_cat_file=shear_cat_file) + self.g1_field, self.g2_field, self.nshear = get_shear_field(config=self.config, shear_cat_file=shear_cat_file) # Get WCS for the focal plane if wcs_fp == None: - wcs_fp = self.focal_plane.getTanWCS(ra_cen, dec_cen, img_rot, chip.pix_scale) + wcs_fp = self.focal_plane.getTanWCS(pointing.ra, pointing.dec, pointing.img_pa, chip.pix_scale) # Create chip Image chip.img = galsim.ImageF(chip.npix_x, chip.npix_y) @@ -84,9 +92,10 @@ class Observation(object): elif chip.survey_type == "spectroscopic": sky_map = calculateSkyMap_split_g(xLen=chip.npix_x, yLen=chip.npix_y, blueLimit=filt.blue_limit, redLimit=filt.red_limit, skyfn=self.path_dict["sky_file"], conf=chip.sls_conf, pixelSize=chip.pix_scale, isAlongY=0) - if pointing_type == 'MS': + # if pointing_type == 'MS': + if pointing.pointing_type == 'MS': # Load catalogues and templates - self.cat = self.Catalog(config=self.config, chip=chip, cat_dir=cat_dir, sed_dir=sed_dir) + self.cat = self.Catalog(config=self.config, chip=chip, pointing=pointing, cat_dir=cat_dir, sed_dir=sed_dir) self.nobj = len(self.cat.objs) # Loop over objects @@ -134,13 +143,13 @@ class Observation(object): if self.config["shear_setting"]["shear_type"] == "constant": if obj.type == 'star': - g1, g2 = 0, 0 + obj.param["g1"], obj.param["g2"] = 0, 0 else: - g1, g2 = self.g1_field, self.g2_field + obj.param["g1"], obj.param["g2"] = self.g1_field, self.g2_field elif self.config["shear_setting"]["shear_type"] == "extra": try: # TODO: every object with individual shear from input catalog(s) - g1, g2 = self.g1_field[j], self.g2_field[j] + obj.param["g1"], obj.param["g2"] = self.g1_field[j], self.g2_field[j] except: print("failed to load external shear.") pass @@ -169,9 +178,10 @@ class Observation(object): bandpass_list=filt.bandpass_sub_list, filt=filt, chip=chip, - g1=g1, - g2=g2, - exptime=exptime) + g1=obj.param["g1"], + g2=obj.param["g2"], + exptime=pointing.exp_time + ) elif chip.survey_type == "spectroscopic" and not self.config["out_cat_only"]: isUpdated, pos_shear = obj.drawObj_slitless( tel=self.tel, @@ -180,15 +190,14 @@ class Observation(object): bandpass_list=filt.bandpass_sub_list, filt=filt, chip=chip, - g1=g1, - g2=g2, - exptime=exptime, - # normFilter=normF, + g1=obj.param["g1"], + g2=obj.param["g2"], + exptime=pointing.exp_time, normFilter=norm_filt, ) if isUpdated: # TODO: add up stats - chip_output.cat_add_obj(obj, pos_img, pos_shear, g1, g2) + chip_output.cat_add_obj(obj, pos_img, pos_shear, obj.param["g1"], obj.param["g2"]) pass else: # print("object omitted", flush=True) @@ -203,49 +212,46 @@ class Observation(object): del psf_model del self.cat - print("check running:1: pointing-{:} chip-{:} pid-{:} memory-{:6.2}GB".format(pointing_ID, chip.chipID, os.getpid(), (psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024 / 1024) ), flush=True) + print("check running:1: pointing-{:} chip-{:} pid-{:} memory-{:6.2}GB".format(pointing.id, chip.chipID, os.getpid(), (psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024 / 1024) ), flush=True) # Detector Effects # =========================================================== - # whether to output zero, dark, flat calibration images. chip.img = chip.addEffects( config=self.config, img=chip.img, chip_output=chip_output, filt=filt, - ra_cen=ra_cen, - dec_cen=dec_cen, - img_rot=img_rot, - pointing_ID=pointing_ID, - timestamp_obs=timestamp_obs, - pointing_type=pointing_type, + ra_cen=pointing.ra, + dec_cen=pointing.dec, + img_rot=pointing.img_pa, + pointing_ID=pointing.id, + timestamp_obs=pointing.timestamp, + pointing_type=pointing.pointing_type, sky_map=sky_map, tel = self.tel) - if pointing_type == 'MS': - datetime_obs = datetime.fromtimestamp(timestamp_obs) + if pointing.pointing_type == 'MS': + datetime_obs = datetime.fromtimestamp(pointing.timestamp) date_obs = datetime_obs.strftime("%y%m%d") time_obs = datetime_obs.strftime("%H%M%S") h_prim = generatePrimaryHeader( xlen=chip.npix_x, ylen=chip.npix_y, - pointNum = str(pointing_ID), - ra=ra_cen, - dec=dec_cen, + pointNum = str(pointing.id), + ra=pointing.ra, + dec=pointing.dec, psize=chip.pix_scale, row_num=chip.rowID, col_num=chip.colID, - # date=self.config["date_obs"], - # time_obs=self.config["time_obs"], date=date_obs, time_obs=time_obs, im_type='MS') h_ext = generateExtensionHeader( xlen=chip.npix_x, ylen=chip.npix_y, - ra=ra_cen, - dec=dec_cen, - pa=img_rot.deg, + ra=pointing.ra, + dec=pointing.dec, + pa=pointing.img_pa.deg, gain=chip.gain, readout=chip.read_noise, dark=chip.dark_noise, @@ -255,7 +261,6 @@ class Observation(object): col_num=chip.colID, extName='raw') chip.img = galsim.Image(chip.img.array, dtype=np.uint16) - # chip.img = galsim.Image(chip.img.array, dtype=np.uint32) hdu1 = fits.PrimaryHDU(header=h_prim) hdu2 = fits.ImageHDU(chip.img.array, header=h_ext) hdu1 = fits.HDUList([hdu1, hdu2]) @@ -266,9 +271,9 @@ class Observation(object): print("# objects that are missed %d out of %d"%(missed_obj, self.nobj)) del chip.img - print("check running:2: pointing-{:} chip-{:} pid-{:} memory-{:6.2}GB".format(pointing_ID, chip.chipID, os.getpid(), (psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024 / 1024) ), flush=True) + print("check running:2: pointing-{:} chip-{:} pid-{:} memory-{:6.2}GB".format(pointing.id, chip.chipID, os.getpid(), (psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024 / 1024) ), flush=True) - def runExposure_MPI_PointingList(self, ra_cen=None, dec_cen=None, pRange=None, timestamp_obs=np.array([1621915200]), pointing_type=np.array(['MS']), img_rot=None, exptime=150., shear_cat_file=None, chips=None, use_mpi=False): + def runExposure_MPI_PointingList(self, pointing_list, shear_cat_file=None, chips=None, use_mpi=False): if use_mpi: comm = MPI.COMM_WORLD ind_thread = comm.Get_rank() @@ -290,30 +295,11 @@ class Observation(object): run_chips.append(chip) run_filts.append(filt) - # TEMP - if len(timestamp_obs) == 1: - timestamp_obs = np.tile(timestamp_obs, len(ra_cen)) - pointing_type = np.tile(pointing_type, len(ra_cen)) - - if pRange is not None: - timestamp_obs = timestamp_obs[pRange] - pointing_type = pointing_type[pRange] - ra_cen = ra_cen[pRange] - dec_cen = dec_cen[pRange] - - # The Starting pointing ID - if pRange is not None: - pStart = pRange[0] - else: - pStart = 0 - - for ipoint in range(len(ra_cen)): + for ipoint in range(len(pointing_list)): for ichip in range(nchips_per_fp): i = ipoint*nchips_per_fp + ichip - if pRange is None: - pointing_ID = pStart + ipoint - else: - pointing_ID = pRange[ipoint] + pointing = pointing_list[ipoint] + pointing_ID = pointing.id if use_mpi: if i % num_thread != ind_thread: continue @@ -322,8 +308,6 @@ class Observation(object): sub_img_dir, prefix = makeSubDir_PointingList(path_dict=self.path_dict, config=self.config, pointing_ID=pointing_ID) - # chip = self.chip_list[ichip] - # filt = self.filter_list[ichip] chip = run_chips[ichip] filt = run_filts[ichip] print("running pointing#%d, chip#%d, at PID#%d..."%(pointing_ID, chip.chipID, pid), flush=True) @@ -332,21 +316,15 @@ class Observation(object): focal_plane=self.focal_plane, chip=chip, filt=filt, - exptime=exptime, - pointing_type=pointing_type[ipoint], + exptime=pointing.exp_time, + pointing_type=pointing.pointing_type, pointing_ID=pointing_ID, subdir=sub_img_dir, prefix=prefix) - self.runOneChip( + self.run_one_chip( chip=chip, filt=filt, chip_output=chip_output, - pointing_ID = pointing_ID, - ra_cen=ra_cen[ipoint], - dec_cen=dec_cen[ipoint], - img_rot=img_rot, - exptime=exptime, - timestamp_obs=timestamp_obs[ipoint], - pointing_type=pointing_type[ipoint], + pointing=pointing, cat_dir=self.path_dict["cat_dir"]) print("finished running chip#%d..."%(chip.chipID), flush=True) diff --git a/ObservationSim/_util.py b/ObservationSim/_util.py index 2e69c3c053c5a1e0b85a2d56c3991fc9cb803d2e..b5514ef50f97f33be69757ab6ffcbb55cc85a628 100755 --- a/ObservationSim/_util.py +++ b/ObservationSim/_util.py @@ -2,6 +2,9 @@ import numpy as np import os from datetime import datetime import argparse +from astropy.time import Time + +from ObservationSim.Config import Pointing def parse_args(): ''' @@ -15,9 +18,32 @@ def parse_args(): parser.add_argument('-w', '--work_dir', help='The path for output.') return parser.parse_args() -def generate_pointings(config, pointing_filename=None, data_dir=None): - pRA = [] - pDEC = [] +def generate_pointing_list(config, pointing_filename=None, data_dir=None): + pointing_list = [] + + # Only valid when the pointing list does not contain time stamp column + t0 = datetime(2021, 5, 25, 12, 0, 0) + delta_t = 10. # Time elapsed between exposures (minutes) + + # Calculate starting time(s) for CAL exposures + # NOTE: temporary implementation + t = datetime.timestamp(t0) + ncal = config['obs_setting']['np_cal'] + ipoint = 0 + for i in range(ncal): + pointing = Pointing( + id = ipoint, + ra=config["obs_setting"]["ra_center"], + dec=config["obs_setting"]["dec_center"], + img_pa=config["obs_setting"]["image_rot"], + timestamp=t, + pointing_type='CAL') + t += 3 * delta_t * 60. # 3 calibration exposures for each pointing + pointing_list.append(pointing) + ipoint += 1 + + run_pointings = config['obs_setting']['run_pointings'] + if pointing_filename and data_dir: pointing_file = os.path.join(data_dir, pointing_filename) f = open(pointing_file, 'r') @@ -25,49 +51,33 @@ def generate_pointings(config, pointing_filename=None, data_dir=None): header = f.readline() iline = 0 for line in f: + if run_pointings and iline not in run_pointings: + iline += 1 + continue line = line.strip() columns = line.split() - pRA.append(float(columns[0])) - pDEC.append(float(columns[1])) + pointing = Pointing() + pointing.read_pointing_columns(columns=columns, id=ipoint, t=t) + t += delta_t * 60. + pointing_list.append(pointing) + iline += 1 + ipoint += 1 f.close() else: - pRA.append(config["obs_setting"]["ra_center"]) - pDEC.append(config["obs_setting"]["dec_center"]) - pRA = np.array(pRA) - pDEC = np.array(pDEC) + pointing = Pointing( + id=ipoint, + ra=config["obs_setting"]["ra_center"], + dec=config["obs_setting"]["dec_center"], + img_pa=config["obs_setting"]["image_rot"], + timestamp=t, + pointing_type='MS' + ) + t += delta_t * 60. + pointing_list.append(pointing) + ipoint += 1 + return pointing_list - # Create calibration pointings - # NOTE: temporary implementation - ncal = config['obs_setting']['np_cal'] - pointing_type = ['MS']*len(pRA) - pRA = np.append([pRA[0]]*ncal, pRA) - pDEC = np.append([pDEC[0]]*ncal, pDEC) - pointing_type = ['CAL']*ncal + pointing_type - - # Calculate starting time(s) - # NOTE: temporary implementation - t0 = datetime(2021, 5, 25, 12, 0, 0) - t = datetime.timestamp(t0) - timestamp_obs = [] - delta_t = 10 # Time elapsed between exposures (minutes) - for i in range(len(pointing_type)): - timestamp_obs.append(t) - if pointing_type[i] == 'CAL': - t += 3 * delta_t * 60 # 3 calibration exposures for each pointing - elif pointing_type[i] == 'MS': - t += delta_t * 60 - timestamp_obs = np.array(timestamp_obs) - pointing_type = np.array(pointing_type) - - if config['obs_setting']['run_pointings'] is None: - pRange = list(range(len(pRA))) - else: - ncal = config['obs_setting']['np_cal'] - plist = config['obs_setting']['run_pointings'] - pRange = list(range(ncal)) + [x + ncal for x in plist] - return pRA, pDEC, timestamp_obs, pointing_type, pRange - -def make_run_dirs(work_dir, run_name, nPointings, pRange=None): +def make_run_dirs(work_dir, run_name, pointing_list): if not os.path.exists(work_dir): try: os.makedirs(work_dir, exist_ok=True) @@ -80,11 +90,8 @@ def make_run_dirs(work_dir, run_name, nPointings, pRange=None): except OSError: pass prefix = "MSC_" - for pointing_ID in range(nPointings): - if pRange is not None: - if pointing_ID not in pRange: - continue - fname=prefix + str(pointing_ID).rjust(7, '0') + for pointing in pointing_list: + fname=prefix + str(pointing.id).rjust(7, '0') subImgDir = os.path.join(imgDir, fname) if not os.path.exists(subImgDir): try: @@ -129,14 +136,13 @@ def makeSubDir_PointingList(path_dict, config, pointing_ID=0): pass return subImgdir, prefix -def getShearFiled(config, shear_cat_file=None): +def get_shear_field(config, shear_cat_file=None): if not config["shear_setting"]["shear_type"] in ["constant", "extra"]: raise ValueError("Please set a right 'shear_method' parameter.") if config["shear_setting"]["shear_type"] == "constant": g1 = config["shear_setting"]["reduced_g1"] g2 = config["shear_setting"]["reduced_g2"] - reduced_shear = np.sqrt(g1**2 + g2**2) nshear = 1 # TODO logging else: diff --git a/config/config_C3.yaml b/config/config_C3.yaml index f4bd7fe300e0267b324d27f2b89e456607e6771c..edbf5334aa75d88b38ad8690615296086c198cf9 100644 --- a/config/config_C3.yaml +++ b/config/config_C3.yaml @@ -12,14 +12,15 @@ # work_dir: "/public/home/fangyuedong/sim_code_release/CSST/test/" work_dir: "/public/home/fangyuedong/20211203/CSST/workplace/" data_dir: "/data/simudata/CSSOSDataProductsSims/data/" -run_name: "C3_20211216" +run_name: "TEST_astrometry" # (Optional) a file of point list # if you just want to run default pointing: # - pointing_dir: null # - pointing_file: null pointing_dir: "/data/simudata/CSSOSDataProductsSims/data/" -pointing_file: "pointing10_20210202.dat" +# pointing_file: "pointing10_20210202.dat" +pointing_file: "pointing_test1214.dat" # Whether to use MPI run_option: @@ -27,14 +28,14 @@ run_option: # NOTE: "n_threads" paramters is currently not used in the backend # simulation codes. It should be implemented later in the web frontend # in order to config the number of threads to request from NAOC cluster - n_threads: 40 + n_threads: 80 # Output catalog only? # If yes, no imaging simulation will run out_cat_only: NO # Only simulate stars? -star_only: NO +star_only: YES # Only simulate galaxies? galaxy_only: NO @@ -68,13 +69,13 @@ obs_setting: # Number of calibration pointings # Note: only valid when a pointing list is specified - np_cal: 3 + np_cal: 0 # Run specific pointing(s): # - 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, 2, 3, 4, 5] + run_pointings: [1, 2, 3] # Run specific chip(s): # - give a list of indexes of chips: [ip_1, ip_2...] @@ -82,6 +83,10 @@ obs_setting: # Note: for all pointings run_chips: [1, 26, 29, 6, 7, 22, 23, 19, 20] + # Whether to enable astrometric modeling + astrometric_lib: "libshao.so" + enable_astrometric_model: True + ############################################### # Input path setting ############################################### diff --git a/run_C3.pbs b/run_C3.pbs index d421f6ab63863cab061294fcecee27654d597626..025f24e62e0aab029987a1590939621fdee6959c 100755 --- a/run_C3.pbs +++ b/run_C3.pbs @@ -7,7 +7,7 @@ ##mpdboot -n 10 -f ./mpd.hosts ##PBS -j oe -#PBS -l nodes=comput108:ppn=80 +#PBS -l nodes=comput112:ppn=80 #####PBS -q longq #PBS -q batch #PBS -u fangyuedong diff --git a/run_sim.py b/run_sim.py index 64447998d30c677443171a294408c8c1ba4fc286..9e5a9ce14d84686a2c605543e859c091e5747cf6 100755 --- a/run_sim.py +++ b/run_sim.py @@ -1,8 +1,6 @@ from ObservationSim.ObservationSim import Observation -from ObservationSim._util import parse_args, generate_pointings, make_run_dirs -from Catalog.Catalog_example import Catalog_example +from ObservationSim._util import parse_args, make_run_dirs, generate_pointing_list import os -import galsim import yaml import gc @@ -11,6 +9,15 @@ gc.enable() def run_sim(Catalog): """ Main method for simulation call. + + Parameters + ---------- + Catalog : Class + a catalog class which is inherited from ObservationSim.MockObject.CatalogBase + + Returns + ---------- + None """ args = parse_args() if args.config_dir is None: @@ -24,7 +31,6 @@ def run_sim(Catalog): print (key + " : " + str(value)) except yaml.YAMLError as exc: print(exc) - config["obs_setting"]["image_rot"] = float(config["obs_setting"]["image_rot"])*galsim.degrees # Overwrite the data and working directories # if they are specified by the command line inputs @@ -33,33 +39,30 @@ def run_sim(Catalog): if args.work_dir is not None: config['work_dir'] = args.work_dir - # Generate lists of RA, DEC, time_stamps, pointing_type (MS or CAL), and pRange - # (selections of pointings to run) based on the input pointing list (or default + # Generate lists pointings based on the input pointing list (or default # pointing RA, DEC) and "config["obs_setting"]["run_pointings"]". # "config['obs_setting']['np_cal']"" is the number of CAL pointings which will be - # appended to the front of all above lists. + # appended to the front. # NOTE: the implementation of gerenating time_stamps is temporary. - pRA, pDEC, timestamp_obs, pointing_type, pRange = generate_pointings(config=config, pointing_filename=config['pointing_file'], data_dir=config['pointing_dir']) + pointing_list = generate_pointing_list(config=config, pointing_filename=config['pointing_file'], data_dir=config['pointing_dir']) # Make the main output directories - make_run_dirs(work_dir=config['work_dir'], run_name=config['run_name'], nPointings=len(pRA), pRange=pRange) + make_run_dirs(work_dir=config['work_dir'], run_name=config['run_name'], pointing_list=pointing_list) # Initialize the simulation obs = Observation(config=config, Catalog=Catalog, work_dir=config['work_dir'], data_dir=config['data_dir']) # Run simulation obs.runExposure_MPI_PointingList( - ra_cen=pRA, - dec_cen=pDEC, - pRange=pRange, - timestamp_obs=timestamp_obs, - pointing_type=pointing_type, - exptime=config["obs_setting"]["exp_time"], + pointing_list=pointing_list, use_mpi=config["run_option"]["use_mpi"], chips=config["obs_setting"]["run_chips"]) if __name__=='__main__': - run_sim(Catalog=Catalog_example) + # To run with the example input catalog + # from Catalog.Catalog_example import Catalog_example + # run_sim(Catalog=Catalog_example) + # To run cycle-3 simulation - # from Catalog.C3Catalog import C3Catalog - # run_sim(Catalog=C3Catalog) + from Catalog.C3Catalog import C3Catalog + run_sim(Catalog=C3Catalog) diff --git a/tests/libshao.so b/tests/libshao.so new file mode 100644 index 0000000000000000000000000000000000000000..47d97a44a74c8f011914dc07f52d8957c81eff39 Binary files /dev/null and b/tests/libshao.so differ diff --git a/tests/shao_test.py b/tests/shao_test.py new file mode 100644 index 0000000000000000000000000000000000000000..a4dab231537136780a79181e36928db1181fc787 --- /dev/null +++ b/tests/shao_test.py @@ -0,0 +1,75 @@ +import os +import sys + +from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position + +def readFits2List(fn): + from astropy.io import fits + if not os.path.exists(fn): + print("Can not find "+fn); return False + + hdul=fits.open(fn) + data = hdul[1].data + ra_list = data['RA'].tolist() + dec_list = data['Dec'].tolist() + pmra_list = data['pmra'].tolist() + pmdec_list = data['pmdec'].tolist() + parallax_list = data['parallax'].tolist() + rv_list = [0.0 for i in range(len(ra_list))] + hdul.close() + return ra_list, dec_list, pmra_list, pmdec_list, rv_list, parallax_list + +def usageExample(targets_fn, results_fn, x, y, z, vx, vy, vz, date_str, time_str): + if not os.path.exists(targets_fn): + print("Can not find " + targets_fn) + sys.exit(1) + ra_list, dec_list, pmra_list, pmdec_list, rv_list, parallax_list = readFits2List(targets_fn) + output_ra_list, output_dec_list = on_orbit_obs_position(ra_list, dec_list, pmra_list, \ + pmdec_list, rv_list, parallax_list, len(ra_list), \ + x, y, z, vx, vy, vz, "J2015.5", date_str, time_str, lib_path='./libshao.so') + f = open(results_fn, "w") + ll = list() + ll.append("n,date,time,ra,dec,pm_ra,pm_dec,rv,parallax,obs_ra,obs_dec\n") + for i in range(len(output_ra_list)): + ra_str = str(ra_list[i]) + dec_str = str(dec_list[i]) + pmra_str = str(pmra_list[i]) + pmdec_str = str(pmdec_list[i]) + rv_str = str(rv_list[i]) + parallax_str = str(parallax_list[i]) + ora_str = str(output_ra_list[i]) + odec_str = str(output_dec_list[i]) + l = str(i)+date_str+","+time_str+","+ra_str+","+dec_str+","+pmra_str+","+pmdec_str+","+rv_str+","+parallax_str+","+ora_str+","+odec_str+"\n" + ll.append(l) + f.writelines(ll) + f.close() + print("Process finished. Results save to "+results_fn) + +def usageTips(): + print("Usage 1: python3 ./shao_test.py") + print("Usage 2: python3 ./shao_test.py yyyy-MM-dd HH:mm:ss.ms Positon1_KM Positon2_KM Positon3_KM Velocity1_MK/S Velocity2_MK/S Velocity3_MK/S $PATH1/MMW_Gaia_Cluster_D20_SS_astrometry.fits $PATH2/results.csv") + print("Caution: Do no include space in path; Unit is KM or KM/S") + print("Example: python3 ./shao_test.py 2025-03-05 22:20:15.12 2347.766100 5132.421392 3726.591334 5.282357 4.644825 -3.074722 ./MMW_Gaia_Cluster_D20_SS_astrometry.fits ./results.csv") + +if __name__ == "__main__": + args_value = sys.argv + if len(args_value) == 1: + usageExample("./MMW_Gaia_Cluster_D20_SS_astrometry.fits", "./results.csv", 2347.766100, 5132.421392, 3726.591334, 5.282357, 4.644825, -3.074722, "2025-03-05", "22:20:15.12"); sys.exit(1) + elif len(args_value) == 11: + date_str = args_value[1] + time_str = args_value[2] + x = float(args_value[3]) + y = float(args_value[4]) + z = float(args_value[5]) + vx = float(args_value[6]) + vy = float(args_value[7]) + vz = float(args_value[8]) + targets_fn = args_value[9] + results_fn = args_value[10] + if not os.path.exists(targets_fn): + print("Can not find " + targets_fn) + sys.exit(1) + usageExample(targets_fn, results_fn, x, y, z, vx, vy, vz, date_str, time_str) + else: + usageTips() +