add orbit position inputs, add astrometric modeling

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]

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"]

ObservationSim/Astrometry/Astrometry_util.py

+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)

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:

ObservationSim/Config/Pointing.py

+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

ObservationSim/Config/__init__.py

-from .Config import ConfigDir, ReadConfig, ParseConfig
-from .ChipOutput import ChipOutput
+from .Config import *
+from .ChipOutput import ChipOutput
\ No newline at end of file
+from .Pointing import Pointing
\ No newline at end of file

ObservationSim/Instrument/FilterParam.py

 import galsim
-import pylab as pl
 import numpy as np
 
 class FilterParam(object):

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
     

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.")

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:

ObservationSim/ObservationSim.py

-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)

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:

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
 ###############################################

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

run_sim.py

 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)

tests/shao_test.py

+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()
+