Merge branch 'develop' into 'master'

release version 2.1.0

See merge request csst_sim/csst-simulation!17

Catalog/C6_50sqdeg.py

+import os
+import galsim
+import random
+import numpy as np
+import h5py as h5
+import healpy as hp
+import astropy.constants as cons
+import traceback
+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 tag_sed, getObservedSED, getABMAG, integrate_sed_bandpass, comoving_dist
+from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position
+
+# (TEST)
+from astropy.cosmology import FlatLambdaCDM
+from astropy import constants
+from astropy import units as U
+from astropy.coordinates import SkyCoord
+from astropy.io import fits
+
+try:
+    import importlib.resources as pkg_resources
+except ImportError:
+    # Try backported to PY<37 'importlib_resources'
+    import importlib_resources as pkg_resources
+
+NSIDE = 128
+
+bundle_file_list = ['galaxies_C6_bundle000199.h5','galaxies_C6_bundle000200.h5','galaxies_C6_bundle000241.h5','galaxies_C6_bundle000242.h5','galaxies_C6_bundle000287.h5','galaxies_C6_bundle000288.h5','galaxies_C6_bundle000714.h5','galaxies_C6_bundle000715.h5','galaxies_C6_bundle000778.h5','galaxies_C6_bundle000779.h5','galaxies_C6_bundle000842.h5','galaxies_C6_bundle000843.h5','galaxies_C6_bundle002046.h5','galaxies_C6_bundle002110.h5','galaxies_C6_bundle002111.h5','galaxies_C6_bundle002173.h5','galaxies_C6_bundle002174.h5','galaxies_C6_bundle002238.h5','galaxies_C6_bundle002596.h5','galaxies_C6_bundle002597.h5','galaxies_C6_bundle002656.h5','galaxies_C6_bundle002657.h5','galaxies_C6_bundle002711.h5','galaxies_C6_bundle002712.h5','galaxies_C6_bundle002844.h5','galaxies_C6_bundle002845.h5','galaxies_C6_bundle002884.h5','galaxies_C6_bundle002885.h5','galaxies_C6_bundle002921.h5','galaxies_C6_bundle002922.h5']
+
+qsosed_file_list = ['quickspeclib_interp1d_run1.fits','quickspeclib_interp1d_run2.fits','quickspeclib_interp1d_run3.fits','quickspeclib_interp1d_run4.fits','quickspeclib_interp1d_run5.fits','quickspeclib_interp1d_run6.fits','quickspeclib_interp1d_run7.fits','quickspeclib_interp1d_run8.fits','quickspeclib_interp1d_run9.fits','quickspeclib_interp1d_run10.fits','quickspeclib_interp1d_run11.fits','quickspeclib_interp1d_run12.fits','quickspeclib_interp1d_run13.fits','quickspeclib_interp1d_run14.fits','quickspeclib_interp1d_run15.fits','quickspeclib_interp1d_run16.fits','quickspeclib_interp1d_run17.fits','quickspeclib_interp1d_run18.fits','quickspeclib_interp1d_run19.fits','quickspeclib_interp1d_run20.fits','quickspeclib_interp1d_run21.fits','quickspeclib_interp1d_run22.fits','quickspeclib_interp1d_run23.fits','quickspeclib_interp1d_run24.fits','quickspeclib_interp1d_run25.fits','quickspeclib_interp1d_run26.fits','quickspeclib_interp1d_run27.fits','quickspeclib_interp1d_run28.fits','quickspeclib_interp1d_run29.fits','quickspeclib_interp1d_run30.fits']
+
+star_file_list = ['C7_Gaia_Galaxia_RA170DECm23_healpix.hdf5', 'C7_Gaia_Galaxia_RA180DECp60_healpix.hdf5', 'C7_Gaia_Galaxia_RA240DECp30_healpix.hdf5', 'C7_Gaia_Galaxia_RA300DECm60_healpix.hdf5', 'C7_Gaia_Galaxia_RA30DECm48_healpix.hdf5']
+star_center_list = [(170., -23.), (180., 60.), (240., 30.), (300., -60.), (30., -48.)]
+
+def get_bundleIndex(healpixID_ring, bundleOrder=4, healpixOrder=7):
+    assert NSIDE == 2**healpixOrder
+    shift = healpixOrder - bundleOrder
+    shift = 2*shift
+
+    nside_bundle = 2**bundleOrder
+    nside_healpix= 2**healpixOrder
+
+    healpixID_nest= hp.ring2nest(nside_healpix, healpixID_ring)
+    bundleID_nest = (healpixID_nest >> shift)
+    bundleID_ring = hp.nest2ring(nside_bundle, bundleID_nest)
+
+    return bundleID_ring
+
+def get_agnsed_file(bundle_file_name):
+    return qsosed_file_list[bundle_file_list.index(bundle_file_name)]
+
+def get_star_cat(ra_pointing, dec_pointing):
+    pointing_c = SkyCoord(ra=ra_pointing*U.deg, dec=dec_pointing*U.deg)
+    max_dist = 10
+    return_star_path = None
+    for star_file, center in zip(star_file_list, star_center_list):
+        center_c = SkyCoord(ra=center[0]*U.deg, dec=center[1]*U.deg)
+        dist = pointing_c.separation(center_c).to(U.deg).value
+        if dist < max_dist:
+            return_star_path = star_file
+            max_dist = dist
+    return return_star_path
+
+class Catalog(CatalogBase):
+    def __init__(self, config, chip, pointing, chip_output, filt, **kwargs):
+        super().__init__()
+        self.cat_dir = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["cat_dir"])
+        self.seed_Av = config["catalog_options"]["seed_Av"]
+
+        self.cosmo = FlatLambdaCDM(H0=67.66, Om0=0.3111)
+
+        self.chip_output = chip_output
+        self.filt = filt
+        self.logger = chip_output.logger
+
+        with pkg_resources.path('Catalog.data', 'SLOAN_SDSS.g.fits') as filter_path:
+            self.normF_star = Table.read(str(filter_path))
+        
+        self.config = config
+        self.chip = chip
+        self.pointing = pointing
+
+        self.max_size = 0.
+
+        if "star_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["star_cat"] and not config["catalog_options"]["galaxy_only"]:
+            # Get the cloest star catalog file
+            star_file_name = get_star_cat(ra_pointing=self.pointing.ra, dec_pointing=self.pointing.dec)
+            star_path = os.path.join(config["catalog_options"]["input_path"]["star_cat"], star_file_name)
+            star_SED_file = config["catalog_options"]["SED_templates_path"]["star_SED"]
+            self.star_path = os.path.join(self.cat_dir, star_path)
+            self.star_SED_path = os.path.join(config["data_dir"], star_SED_file)
+            self._load_SED_lib_star()
+        
+        if "galaxy_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["galaxy_cat"] and not config["catalog_options"]["star_only"]:
+            galaxy_dir = config["catalog_options"]["input_path"]["galaxy_cat"]
+            self.galaxy_path = os.path.join(self.cat_dir, galaxy_dir)
+            self.galaxy_SED_path = os.path.join(config["data_dir"], config["catalog_options"]["SED_templates_path"]["galaxy_SED"])
+            self._load_SED_lib_gals()
+            self.agn_seds = {}
+
+        if "AGN_SED" in config["catalog_options"]["SED_templates_path"] and not config["catalog_options"]["star_only"]:
+            self.AGN_SED_path = os.path.join(config["data_dir"], config["catalog_options"]["SED_templates_path"]["AGN_SED"])
+
+        if "rotateEll" in config["catalog_options"]:
+            self.rotation = np.radians(float(config["catalog_options"]["rotateEll"]))
+        else:
+            self.rotation = 0.
+
+        # Update output .cat header with catalog specific output columns
+        self._add_output_columns_header()
+
+        self._get_healpix_list()
+        self._load()
+    
+    def _add_output_columns_header(self):
+        self.add_hdr = " model_tag teff logg feh"
+        self.add_hdr += " bulgemass diskmass detA e1 e2 kappa g1 g2 size galType veldisp "
+
+        self.add_fmt = " %10s %8.4f %8.4f %8.4f"
+        self.add_fmt += " %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %4d %8.4f "
+        self.chip_output.update_ouptut_header(additional_column_names=self.add_hdr)
+
+    def _get_healpix_list(self):
+        self.sky_coverage = self.chip.getSkyCoverageEnlarged(self.chip.img.wcs, margin=0.2)
+        ra_min, ra_max, dec_min, dec_max = self.sky_coverage.xmin, self.sky_coverage.xmax, self.sky_coverage.ymin, self.sky_coverage.ymax
+        ra = np.deg2rad(np.array([ra_min, ra_max, ra_max, ra_min]))
+        dec = np.deg2rad(np.array([dec_max, dec_max, dec_min, dec_min]))
+        self.pix_list = hp.query_polygon(
+            NSIDE,
+            hp.ang2vec(np.radians(90.) - dec, ra),
+            inclusive=True
+        )
+        if self.logger is not None:
+            msg = str(("HEALPix List: ", self.pix_list))
+            self.logger.info(msg)
+        else:
+            print("HEALPix List: ", self.pix_list)
+
+    def load_norm_filt(self, obj):
+        if obj.type == "star":
+            return self.normF_star
+        elif obj.type == "galaxy" or obj.type == "quasar":
+            # return self.normF_galaxy
+            return None
+        else:
+            return None
+
+    def _load_SED_lib_star(self):
+        self.tempSED_star = h5.File(self.star_SED_path,'r')
+
+    def _load_SED_lib_gals(self):
+        pcs = h5.File(os.path.join(self.galaxy_SED_path, "pcs.h5"), "r")
+        lamb = h5.File(os.path.join(self.galaxy_SED_path, "lamb.h5"), "r")
+        self.lamb_gal = lamb['lamb'][()]
+        self.pcs = pcs['pcs'][()]
+
+    def _load_gals(self, gals, pix_id=None, cat_id=0, agnsed_file=""):
+        ngals = len(gals['ra'])
+
+        # Apply astrometric modeling
+        ra_arr = gals['ra'][:]
+        dec_arr = gals['dec'][:]
+        if self.config["obs_setting"]["enable_astrometric_model"]:
+            ra_list = ra_arr.tolist()
+            dec_list = dec_arr.tolist()
+            pmra_list = np.zeros(ngals).tolist()
+            pmdec_list = np.zeros(ngals).tolist()
+            rv_list = np.zeros(ngals).tolist()
+            parallax_list = [1e-9] * ngals
+            dt = datetime.utcfromtimestamp(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=ngals,
+                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="J2000",
+                input_date_str=date_str,
+                input_time_str=time_str
+            )
+
+        for igals in range(ngals):
+            # # (TEST)
+            # if igals > 100:
+            #     break
+            
+            param = self.initialize_param()
+            param['ra'] = ra_arr[igals]
+            param['dec'] = dec_arr[igals]
+            param['ra_orig'] = gals['ra'][igals]
+            param['dec_orig'] = gals['dec'][igals]
+            # param['mag_use_normal'] = gals['mag_csst_%s'%(self.filt.filter_type)][igals]
+            if self.filt.filter_type == 'NUV':
+                param['mag_use_normal'] = gals['mag_csst_nuv'][igals]
+            else:
+                param['mag_use_normal'] = gals['mag_csst_%s'%(self.filt.filter_type)][igals]
+            if self.filt.is_too_dim(mag=param['mag_use_normal'], margin=self.config["obs_setting"]["mag_lim_margin"]):
+                continue
+
+            param['z'] = gals['redshift'][igals]
+            param['model_tag'] = 'None'
+            param['g1'] = gals['shear'][igals][0]
+            param['g2'] = gals['shear'][igals][1]
+            param['kappa'] = gals['kappa'][igals]
+            param['e1'] = gals['ellipticity_true'][igals][0]
+            param['e2'] = gals['ellipticity_true'][igals][1]
+            
+            # For shape calculation
+            param['e1'], param['e2'], param['ell_total'] = self.rotate_ellipticity(
+                                                                    e1=gals['ellipticity_true'][igals][0],
+                                                                    e2=gals['ellipticity_true'][igals][1],
+                                                                    rotation=self.rotation,
+                                                                    unit='radians')
+            # param['ell_total'] = np.sqrt(param['e1']**2 + param['e2']**2)
+            if param['ell_total'] > 0.9:
+                continue
+            # phi_e = cmath.phase(complex(param['e1'], param['e2']))
+            # param['e1'] = param['ell_total'] * np.cos(phi_e + 2*self.rotation)
+            # param['e2'] = param['ell_total'] * np.sin(phi_e + 2*self.rotation)
+            
+            param['e1_disk'] = param['e1']
+            param['e2_disk'] = param['e2']
+            param['e1_bulge'] = param['e1']
+            param['e2_bulge'] = param['e2']
+
+
+            param['delta_ra'] = 0
+            param['delta_dec'] = 0
+
+            # Masses
+            param['bulgemass'] = gals['bulgemass'][igals]
+            param['diskmass'] = gals['diskmass'][igals]
+
+            param['size'] = gals['size'][igals]
+            if param['size'] > self.max_size:
+                self.max_size = param['size']
+
+            # Sersic index
+            param['disk_sersic_idx'] = 1.
+            param['bulge_sersic_idx'] = 4.
+
+            # Sizes
+            param['bfrac'] = param['bulgemass']/(param['bulgemass'] + param['diskmass'])
+            if param['bfrac'] >= 0.6:
+                param['hlr_bulge'] = param['size']
+                param['hlr_disk'] = param['size'] * (1. - param['bfrac'])
+            else:
+                param['hlr_disk'] = param['size']
+                param['hlr_bulge'] = param['size'] * param['bfrac']
+
+            # SED coefficients
+            param['coeff'] = gals['coeff'][igals]
+            param['detA'] = gals['detA'][igals]
+
+            # Others
+            param['galType'] = gals['type'][igals]
+            param['veldisp'] = gals['veldisp'][igals]
+            
+            # TEST no redening and no extinction
+            param['av'] = 0.0
+            param['redden'] = 0
+
+            # Is this an Quasar?
+            param['qsoindex'] = gals['qsoindex'][igals]
+            if param['qsoindex'] == -1:
+                param['star'] = 0   # Galaxy
+                param['agnsed_file'] = ""
+            else:
+                param['star'] = 2   # Quasar
+                param['agnsed_file'] = agnsed_file
+
+            # NOTE: this cut cannot be put before the SED type has been assigned
+            if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200):
+                continue
+
+            # TEMP
+            self.ids += 1
+            param['id'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(igals)
+            
+            if param['star'] == 0:
+                obj = Galaxy(param, logger=self.logger)
+            elif param['star'] == 2:
+                obj = Quasar(param, logger=self.logger)
+            
+            # Need to deal with additional output columns
+            obj.additional_output_str = self.add_fmt%("n", 0., 0., 0.,
+                                                    param['bulgemass'], param['diskmass'], param['detA'],
+                                                    param['e1'], param['e2'], param['kappa'], param['g1'], param['g2'], param['size'],
+                                                    param['galType'], param['veldisp'])
+            
+            self.objs.append(obj)
+
+    def _load_stars(self, stars, pix_id=None):
+        nstars = len(stars['sourceID'])
+        # Apply astrometric modeling
+        ra_arr = stars["RA"][:]
+        dec_arr = stars["Dec"][:]
+        pmra_arr = stars['pmra'][:]
+        pmdec_arr = stars['pmdec'][:]
+        rv_arr = stars['RV'][:]
+        parallax_arr = stars['parallax'][:]
+        if self.config["obs_setting"]["enable_astrometric_model"]:
+            ra_list = ra_arr.tolist()
+            dec_list = dec_arr.tolist()
+            pmra_list = pmra_arr.tolist()
+            pmdec_list = pmdec_arr.tolist()
+            rv_list = rv_arr.tolist()
+            parallax_list = parallax_arr.tolist()
+            dt = datetime.utcfromtimestamp(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="J2000",
+                input_date_str=date_str,
+                input_time_str=time_str
+            )
+        for istars in range(nstars):
+            # # (TEST)
+            # if istars > 100:
+            #     break
+
+            param = self.initialize_param()
+            param['ra'] = ra_arr[istars]
+            param['dec'] = dec_arr[istars]
+            param['ra_orig'] = stars["RA"][istars]
+            param['dec_orig'] = stars["Dec"][istars]
+            param['pmra'] = pmra_arr[istars]
+            param['pmdec'] = pmdec_arr[istars]
+            param['rv'] = rv_arr[istars]
+            param['parallax'] = parallax_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]
+            self.ids += 1
+            param['id'] = stars['sourceID'][istars]
+            param['sed_type'] = stars['sourceID'][istars]
+            param['model_tag'] = stars['model_tag'][istars]
+            param['teff'] = stars['teff'][istars]
+            param['logg'] = stars['grav'][istars]
+            param['feh'] = stars['feh'][istars]
+            param['z'] = 0.0
+            param['star'] = 1   # Star
+            obj = Star(param, logger=self.logger)
+
+            # Append additional output columns to the .cat file
+            obj.additional_output_str = self.add_fmt%(param["model_tag"], param['teff'], param['logg'], param['feh'],
+                                                    0., 0., 0., 0., 0., 0., 0., 0., 0., -1, 0.)
+
+            self.objs.append(obj)
+
+    def _load(self, **kwargs):
+        self.objs = []
+        self.ids = 0
+        
+        if "star_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["star_cat"] and not self.config["catalog_options"]["galaxy_only"]:
+            star_cat = h5.File(self.star_path, 'r')['catalog']
+            for pix in self.pix_list:
+                try:
+                    stars = star_cat[str(pix)]
+                    self._load_stars(stars, pix_id=pix)
+                    del stars
+                except Exception as e:
+                    self.logger.error(str(e))
+                    print(e)
+        
+        if "galaxy_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["galaxy_cat"] and not self.config["catalog_options"]["star_only"]:
+            for pix in self.pix_list:
+                try:
+                    bundleID  = get_bundleIndex(pix)
+                    bundle_file = "galaxies_C6_bundle{:06}.h5".format(bundleID)
+                    file_path = os.path.join(self.galaxy_path, bundle_file)
+                    gals_cat = h5.File(file_path, 'r')['galaxies']
+                    gals = gals_cat[str(pix)]
+
+                    # Get corresponding AGN SED file
+                    agnsed_file = get_agnsed_file(bundle_file)
+                    agnsed_path = os.path.join(self.AGN_SED_path, agnsed_file)
+                    self.agn_seds[agnsed_file] = fits.open(agnsed_path)[0].data
+
+                    self._load_gals(gals, pix_id=pix, cat_id=bundleID, agnsed_file=agnsed_file)
+
+                    del gals
+                except Exception as e:
+                    traceback.print_exc()
+                    self.logger.error(str(e))
+                    print(e)
+
+        if self.logger is not None:
+            self.logger.info("maximum galaxy size: %.4f"%(self.max_size))
+            self.logger.info("number of objects in catalog: %d"%(len(self.objs)))
+        else:
+            print("number of objects in catalog: ", len(self.objs))
+
+    def load_sed(self, obj, **kwargs):
+        if obj.type == 'star':
+            _, wave, flux = tag_sed(
+                h5file=self.tempSED_star,
+                model_tag=obj.param['model_tag'],
+                teff=obj.param['teff'],
+                logg=obj.param['logg'],
+                feh=obj.param['feh']
+            )
+        elif obj.type == 'galaxy' or obj.type == 'quasar':
+            factor = 10**(-.4 * self.cosmo.distmod(obj.z).value)
+            if obj.type == 'galaxy':
+                flux = np.matmul(self.pcs, obj.coeff) * factor
+                #  if np.any(flux < 0):
+                #     raise ValueError("Glaxy %s: negative SED fluxes"%obj.id)
+                flux[flux < 0] = 0.
+                sedcat = np.vstack((self.lamb_gal, flux)).T
+                sed_data = getObservedSED(
+                    sedCat=sedcat,
+                    redshift=obj.z,
+                    av=obj.param["av"],
+                    redden=obj.param["redden"]
+                )
+                wave, flux = sed_data[0], sed_data[1]
+            elif obj.type == 'quasar':
+                flux = self.agn_seds[obj.agnsed_file][int(obj.qsoindex)] * 1e-17
+                flux[flux < 0] = 0.
+                wave = self.lamb_gal
+        else:
+            raise ValueError("Object type not known")
+        speci = interpolate.interp1d(wave, flux)
+        lamb = np.arange(2000, 11001+0.5, 0.5)
+        y = speci(lamb)
+        # erg/s/cm2/A --> photon/s/m2/A
+        all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13
+        sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX'))
+        
+        if obj.type == 'quasar':
+            # integrate to get the magnitudes
+            sed_photon = np.array([sed['WAVELENGTH'], sed['FLUX']]).T
+            sed_photon = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(sed_photon[:, 1]), interpolant='nearest')
+            sed_photon = galsim.SED(sed_photon, wave_type='A', flux_type='1', fast=False)
+            interFlux = integrate_sed_bandpass(sed=sed_photon, bandpass=self.filt.bandpass_full)
+            obj.param['mag_use_normal'] = getABMAG(interFlux, self.filt.bandpass_full)
+            # mag = getABMAG(interFlux, self.filt.bandpass_full)
+            # print("mag diff = %.3f"%(mag - obj.param['mag_use_normal']))
+        del wave
+        del flux
+        return sed

Catalog/C6_Catalog.py

@@ -83,7 +83,7 @@ class Catalog(CatalogBase):
             self._load_SED_lib_AGN()
 
         if "rotateEll" in config["catalog_options"]:
-            self.rotation = float(int(config["catalog_options"]["rotateEll"]/45.))
+            self.rotation = np.radians(float(config["catalog_options"]["rotateEll"]))
         else:
             self.rotation = 0.
 
@@ -259,7 +259,7 @@ class Catalog(CatalogBase):
             param['id'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(igals)
             
             if param['star'] == 0:
-                obj = Galaxy(param, self.rotation, logger=self.logger)
+                obj = Galaxy(param, logger=self.logger)
             
             # Need to deal with additional output columns
             obj.additional_output_str = self.add_fmt%("n", 0., 0., 0.,

Catalog/C6_fits_Catalog.py

@@ -94,7 +94,7 @@ class Catalog(CatalogBase):
         ###mock_stamp_END
         
         if "rotateEll" in config["catalog_options"]:
-            self.rotation = float(int(config["catalog_options"]["rotateEll"]/45.))
+            self.rotation = np.radians(float(config["catalog_options"]["rotateEll"]))
         else:
             self.rotation = 0.
 
@@ -272,7 +272,7 @@ class Catalog(CatalogBase):
             param['id'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(igals)
             
             if param['star'] == 0:
-                obj = Galaxy(param, self.rotation, logger=self.logger)
+                obj = Galaxy(param, logger=self.logger)
             
             # Need to deal with additional output columns
             obj.additional_output_str = self.add_fmt%("n", 0., 0., 0.,

Catalog/FGS_Catalog.py

@@ -84,7 +84,7 @@ class Catalog(CatalogBase):
             self._load_SED_lib_AGN()
 
         if "rotateEll" in config["catalog_options"]:
-            self.rotation = float(int(config["catalog_options"]["rotateEll"]/45.))
+            self.rotation = np.radians(float(config["catalog_options"]["rotateEll"]))
         else:
             self.rotation = 0.
 
@@ -260,7 +260,7 @@ class Catalog(CatalogBase):
             param['id'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(igals)
             
             if param['star'] == 0:
-                obj = Galaxy(param, self.rotation, logger=self.logger)
+                obj = Galaxy(param, logger=self.logger)
             
             # Need to deal with additional output columns
             obj.additional_output_str = self.add_fmt%("n", 0., 0., 0.,

Catalog/NGPCatalog.py

@@ -52,7 +52,7 @@ class Catalog(CatalogBase):
             self.galaxy_SED_path = os.path.join(config["data_dir"], config["catalog_options"]["SED_templates_path"]["galaxy_SED"])
             self._load_SED_lib_gals()
         if "rotateEll" in config["catalog_options"]:
-            self.rotation = float(int(config["catalog_options"]["rotateEll"]/45.))
+            self.rotation = np.radians(float(config["catalog_options"]["rotateEll"]))
         else:
             self.rotation = 0.
 
@@ -191,7 +191,7 @@ class Catalog(CatalogBase):
             param['id'] = gals['galaxyID'][igals]
             
             if param['star'] == 0:
-                obj = Galaxy(param, self.rotation, logger=self.logger)
+                obj = Galaxy(param, logger=self.logger)
             if param['star'] == 2:
                 obj = Quasar(param, logger=self.logger)
             

Catalog/fd_test_C6.py

+import os
+import galsim
+import random
+import numpy as np
+import h5py as h5
+import healpy as hp
+import astropy.constants as cons
+import traceback
+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 tag_sed, getObservedSED, getABMAG, integrate_sed_bandpass, comoving_dist
+from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position
+
+# (TEST)
+from astropy.cosmology import FlatLambdaCDM
+from astropy import constants
+from astropy import units as U
+
+try:
+    import importlib.resources as pkg_resources
+except ImportError:
+    # Try backported to PY<37 'importlib_resources'
+    import importlib_resources as pkg_resources
+
+NSIDE = 128
+
+def get_bundleIndex(healpixID_ring, bundleOrder=4, healpixOrder=7):
+    assert NSIDE == 2**healpixOrder
+    shift = healpixOrder - bundleOrder
+    shift = 2*shift
+
+    nside_bundle = 2**bundleOrder
+    nside_healpix= 2**healpixOrder
+
+    healpixID_nest= hp.ring2nest(nside_healpix, healpixID_ring)
+    bundleID_nest = (healpixID_nest >> shift)
+    bundleID_ring = hp.nest2ring(nside_bundle, bundleID_nest)
+
+    return bundleID_ring
+
+class Catalog(CatalogBase):
+    def __init__(self, config, chip, pointing, chip_output, filt, **kwargs):
+        super().__init__()
+        self.cat_dir = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["cat_dir"])
+        self.seed_Av = config["catalog_options"]["seed_Av"]
+
+        self.cosmo = FlatLambdaCDM(H0=67.66, Om0=0.3111)
+
+        self.chip_output = chip_output
+        self.filt = filt
+        self.logger = chip_output.logger
+
+        with pkg_resources.path('Catalog.data', 'SLOAN_SDSS.g.fits') as filter_path:
+            self.normF_star = Table.read(str(filter_path))
+        
+        self.config = config
+        self.chip = chip
+        self.pointing = pointing
+
+        self.max_size = 0.
+
+        if "star_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["star_cat"] and not config["catalog_options"]["galaxy_only"]:
+            star_file = config["catalog_options"]["input_path"]["star_cat"]
+            star_SED_file = config["catalog_options"]["SED_templates_path"]["star_SED"]
+            self.star_path = os.path.join(self.cat_dir, star_file)
+            self.star_SED_path = os.path.join(config["data_dir"], star_SED_file)
+            self._load_SED_lib_star()
+        if "galaxy_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["galaxy_cat"] and not config["catalog_options"]["star_only"]:
+            galaxy_dir = config["catalog_options"]["input_path"]["galaxy_cat"]
+            self.galaxy_path = os.path.join(self.cat_dir, galaxy_dir)
+            self.galaxy_SED_path = os.path.join(config["data_dir"], config["catalog_options"]["SED_templates_path"]["galaxy_SED"])
+            self._load_SED_lib_gals()
+
+        if "AGN_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["AGN_cat"] and not config["catalog_options"]["star_only"]:
+            AGN_dir = config["catalog_options"]["input_path"]["AGN_cat"]
+            self.AGN_path = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["AGN_cat"])
+            self.AGN_SED_path = os.path.join(config["data_dir"], config["catalog_options"]["SED_templates_path"]["AGN_SED"])
+            self.AGN_SED_wave_path = os.path.join(config['data_dir'], config["catalog_options"]["SED_templates_path"]["AGN_SED_WAVE"])
+            self._load_SED_lib_AGN()
+
+        if "rotateEll" in config["catalog_options"]:
+            self.rotation = np.radians(float(config["catalog_options"]["rotateEll"]))
+        else:
+            self.rotation = 0.
+
+        # Update output .cat header with catalog specific output columns
+        self._add_output_columns_header()
+
+        self._get_healpix_list()
+        self._load()
+    
+    def _add_output_columns_header(self):
+        self.add_hdr = " model_tag teff logg feh"
+        self.add_hdr += " bulgemass diskmass detA e1 e2 kappa g1 g2 size galType veldisp "
+
+        self.add_fmt = " %10s %8.4f %8.4f %8.4f"
+        self.add_fmt += " %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %4d %8.4f "
+        self.chip_output.update_ouptut_header(additional_column_names=self.add_hdr)
+
+    def _get_healpix_list(self):
+        self.sky_coverage = self.chip.getSkyCoverageEnlarged(self.chip.img.wcs, margin=0.2)
+        ra_min, ra_max, dec_min, dec_max = self.sky_coverage.xmin, self.sky_coverage.xmax, self.sky_coverage.ymin, self.sky_coverage.ymax
+        ra = np.deg2rad(np.array([ra_min, ra_max, ra_max, ra_min]))
+        dec = np.deg2rad(np.array([dec_max, dec_max, dec_min, dec_min]))
+        # vertices = spherical_to_cartesian(1., dec, ra)
+        self.pix_list = hp.query_polygon(
+            NSIDE,
+            hp.ang2vec(np.radians(90.) - dec, ra),
+            inclusive=True
+        )
+        if self.logger is not None:
+            msg = str(("HEALPix List: ", self.pix_list))
+            self.logger.info(msg)
+        else:
+            print("HEALPix List: ", self.pix_list)
+
+    def load_norm_filt(self, obj):
+        if obj.type == "star":
+            return self.normF_star
+        elif obj.type == "galaxy" or obj.type == "quasar":
+            # return self.normF_galaxy
+            return None
+        else:
+            return None
+
+    def _load_SED_lib_star(self):
+        self.tempSED_star = h5.File(self.star_SED_path,'r')
+
+    def _load_SED_lib_gals(self):
+        pcs = h5.File(os.path.join(self.galaxy_SED_path, "pcs.h5"), "r")
+        lamb = h5.File(os.path.join(self.galaxy_SED_path, "lamb.h5"), "r")
+        self.lamb_gal = lamb['lamb'][()]
+        self.pcs = pcs['pcs'][()]
+
+    def _load_SED_lib_AGN(self):
+        from astropy.io import fits
+        self.SED_AGN = fits.open(self.AGN_SED_path)[0].data
+        self.lamb_AGN = np.load(self.AGN_SED_wave_path)
+
+
+    def _load_gals(self, gals, pix_id=None, cat_id=0):
+        ngals = len(gals['ra'])
+
+        # Apply astrometric modeling
+        # in C3 case only aberration
+        ra_arr = gals['ra'][:]
+        dec_arr = gals['dec'][:]
+        if self.config["obs_setting"]["enable_astrometric_model"]:
+            ra_list = ra_arr.tolist()
+            dec_list = dec_arr.tolist()
+            pmra_list = np.zeros(ngals).tolist()
+            pmdec_list = np.zeros(ngals).tolist()
+            rv_list = np.zeros(ngals).tolist()
+            parallax_list = [1e-9] * ngals
+            dt = datetime.utcfromtimestamp(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=ngals,
+                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="J2000",
+                input_date_str=date_str,
+                input_time_str=time_str
+            )
+
+        # for igals in range(ngals):
+        for igals in range(0, ngals, 5):
+            # # (TEST)
+            # if igals > 1000:
+            #     break
+            
+            param = self.initialize_param()
+            param['ra'] = ra_arr[igals]
+            param['dec'] = dec_arr[igals]
+            param['ra_orig'] = gals['ra'][igals]
+            param['dec_orig'] = gals['dec'][igals]
+            # param['mag_use_normal'] = gals['mag_csst_%s'%(self.filt.filter_type)][igals]
+            param['mag_use_normal'] = 20.
+            if self.filt.is_too_dim(mag=param['mag_use_normal'], margin=self.config["obs_setting"]["mag_lim_margin"]):
+                continue
+
+            param['z'] = gals['redshift'][igals]
+            param['model_tag'] = 'None'
+            # param['g1'] = gals['shear'][igals][0]
+            # param['g2'] = gals['shear'][igals][1]
+            param['g1'] = 0.
+            param['g2'] = 0.
+            param['kappa'] = gals['kappa'][igals]
+            param['e1'] = gals['ellipticity_true'][igals][0]
+            param['e2'] = gals['ellipticity_true'][igals][1]
+            
+            # For shape calculation
+            
+            param['ell_total'] = np.sqrt(param['e1']**2 + param['e2']**2)
+            if param['ell_total'] > 0.9:
+                continue
+            # param['e1_disk'] = param['e1']
+            # param['e2_disk'] = param['e2']
+            # param['e1_bulge'] = param['e1']
+            # param['e2_bulge'] = param['e2']
+            param['e1_disk'] = 0.
+            param['e2_disk'] = 0.
+            param['e1_bulge'] = 0.
+            param['e2_bulge'] = 0.
+
+
+            param['delta_ra'] = 0
+            param['delta_dec'] = 0
+
+            # Masses
+            param['bulgemass'] = gals['bulgemass'][igals]
+            param['diskmass'] = gals['diskmass'][igals]
+
+            # param['size'] = gals['size'][igals]
+            param['size'] = 1.
+            if param['size'] > self.max_size:
+                self.max_size = param['size']
+
+            # Sersic index
+            param['disk_sersic_idx'] = 1.
+            param['bulge_sersic_idx'] = 4.
+
+            # Sizes
+            # param['bfrac'] = param['bulgemass']/(param['bulgemass'] + param['diskmass'])
+            param['bfrac'] = 0.
+            if param['bfrac'] >= 0.6:
+                param['hlr_bulge'] = param['size']
+                param['hlr_disk'] = param['size'] * (1. - param['bfrac'])
+            else:
+                param['hlr_disk'] = param['size']
+                param['hlr_bulge'] = param['size'] * param['bfrac']
+
+            # SED coefficients
+            param['coeff'] = gals['coeff'][igals]
+            param['detA'] = gals['detA'][igals]
+
+            # Others
+            param['galType'] = gals['type'][igals]
+            param['veldisp'] = gals['veldisp'][igals]
+            
+            # TEST no redening and no extinction
+            param['av'] = 0.0
+            param['redden'] = 0
+
+            param['star'] = 0   # Galaxy
+
+            # NOTE: this cut cannot be put before the SED type has been assigned
+            if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200):
+                continue
+
+            # TEMP
+            self.ids += 1
+            # param['id'] = self.ids
+            param['id'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(igals)
+            
+            if param['star'] == 0:
+                obj = Galaxy(param, logger=self.logger)
+            
+            # Need to deal with additional output columns
+            obj.additional_output_str = self.add_fmt%("n", 0., 0., 0.,
+                                                    param['bulgemass'], param['diskmass'], param['detA'],
+                                                    param['e1'], param['e2'], param['kappa'], param['g1'], param['g2'], param['size'],
+                                                    param['galType'], param['veldisp'])
+            
+            self.objs.append(obj)
+
+    def _load_stars(self, stars, pix_id=None):
+        nstars = len(stars['sourceID'])
+        # Apply astrometric modeling
+        ra_arr = stars["RA"][:]
+        dec_arr = stars["Dec"][:]
+        pmra_arr = stars['pmra'][:]
+        pmdec_arr = stars['pmdec'][:]
+        rv_arr = stars['RV'][:]
+        parallax_arr = stars['parallax'][:]
+        if self.config["obs_setting"]["enable_astrometric_model"]:
+            ra_list = ra_arr.tolist()
+            dec_list = dec_arr.tolist()
+            pmra_list = pmra_arr.tolist()
+            pmdec_list = pmdec_arr.tolist()
+            rv_list = rv_arr.tolist()
+            parallax_list = parallax_arr.tolist()
+            dt = datetime.utcfromtimestamp(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="J2000",
+                input_date_str=date_str,
+                input_time_str=time_str
+            )
+        for istars in range(nstars):
+            # # (TEST)
+            # if istars > 100:
+            #     break
+
+            param = self.initialize_param()
+            param['ra'] = ra_arr[istars]
+            param['dec'] = dec_arr[istars]
+            param['ra_orig'] = stars["RA"][istars]
+            param['dec_orig'] = stars["Dec"][istars]
+            param['pmra'] = pmra_arr[istars]
+            param['pmdec'] = pmdec_arr[istars]
+            param['rv'] = rv_arr[istars]
+            param['parallax'] = parallax_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]
+            # if param['mag_use_normal'] >= 26.5:
+            #     continue
+            self.ids += 1
+            param['id'] = stars['sourceID'][istars]
+            param['sed_type'] = stars['sourceID'][istars]
+            param['model_tag'] = stars['model_tag'][istars]
+            param['teff'] = stars['teff'][istars]
+            param['logg'] = stars['grav'][istars]
+            param['feh'] = stars['feh'][istars]
+            param['z'] = 0.0
+            param['star'] = 1   # Star
+            obj = Star(param, logger=self.logger)
+
+            # Append additional output columns to the .cat file
+            obj.additional_output_str = self.add_fmt%(param["model_tag"], param['teff'], param['logg'], param['feh'],
+                                                    0., 0., 0., 0., 0., 0., 0., 0., 0., -1, 0.)
+
+            self.objs.append(obj)
+
+    def _load_AGNs(self):
+        data = Table.read(self.AGN_path)
+        ra_arr = data['ra']
+        dec_arr = data['dec']
+        nAGNs = len(data)
+        if self.config["obs_setting"]["enable_astrometric_model"]:
+            ra_list = ra_arr.tolist()
+            dec_list = dec_arr.tolist()
+            pmra_list = np.zeros(nAGNs).tolist()
+            pmdec_list = np.zeros(nAGNs).tolist()
+            rv_list = np.zeros(nAGNs).tolist()
+            parallax_list = [1e-9] * nAGNs
+            dt = datetime.utcfromtimestamp(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=nAGNs,
+                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="J2000",
+                input_date_str=date_str,
+                input_time_str=time_str
+            )
+        for iAGNs in range(nAGNs):
+            param = self.initialize_param()
+            param['ra'] = ra_arr[iAGNs]
+            param['dec'] = dec_arr[iAGNs]
+            param['ra_orig'] = data['ra'][iAGNs]
+            param['dec_orig'] = data['dec'][iAGNs]
+            param['z'] = data['z'][iAGNs]
+            param['appMag'] = data['appMag'][iAGNs]
+            param['absMag'] = data['absMag'][iAGNs]
+            
+            # NOTE: this cut cannot be put before the SED type has been assigned
+            if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200):
+                continue
+
+            # TEST no redening and no extinction
+            param['av'] = 0.0
+            param['redden'] = 0
+
+            param['star'] = 2   # Quasar
+            param['id'] = data['igmlos'][iAGNs]
+
+            if param['star'] == 2:
+                obj = Quasar(param, logger=self.logger)
+
+            # Append additional output columns to the .cat file
+            obj.additional_output_str = self.add_fmt%("n", 0., 0., 0.,
+                                                    0., 0., 0., 0., 0., 0., 0., 0., 0., -1, 0.)
+            self.objs.append(obj)
+
+    def _load(self, **kwargs):
+        self.objs = []
+        self.ids = 0
+        
+        # if "star_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["star_cat"] and not self.config["catalog_options"]["galaxy_only"]:
+        #     star_cat = h5.File(self.star_path, 'r')['catalog']
+        #     for pix in self.pix_list:
+        #         try:
+        #             stars = star_cat[str(pix)]
+        #             self._load_stars(stars, pix_id=pix)
+        #             del stars
+        #         except Exception as e:
+        #             self.logger.error(str(e))
+        #             print(e)
+        
+        if "galaxy_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["galaxy_cat"] and not self.config["catalog_options"]["star_only"]:
+            for pix in self.pix_list:
+                try:
+                    bundleID  = get_bundleIndex(pix)
+                    file_path = os.path.join(self.galaxy_path, "galaxies_C6_bundle{:06}.h5".format(bundleID))
+                    gals_cat = h5.File(file_path, 'r')['galaxies']
+                    gals = gals_cat[str(pix)]
+                    self._load_gals(gals, pix_id=pix, cat_id=bundleID)
+                    del gals
+                except Exception as e:
+                    traceback.print_exc()
+                    self.logger.error(str(e))
+                    print(e)
+
+        # if "AGN_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["AGN_cat"] and not self.config["catalog_options"]["star_only"]:
+        #     try:
+        #         self._load_AGNs()
+        #     except Exception as e:
+        #         traceback.print_exc()
+        #         self.logger.error(str(e))
+        #         print(e)
+
+        if self.logger is not None:
+            self.logger.info("maximum galaxy size: %.4f"%(self.max_size))
+            self.logger.info("number of objects in catalog: %d"%(len(self.objs)))
+        else:
+            print("number of objects in catalog: ", len(self.objs))
+
+    def load_sed(self, obj, **kwargs):
+        if obj.type == 'star':
+            _, wave, flux = tag_sed(
+                h5file=self.tempSED_star,
+                model_tag=obj.param['model_tag'],
+                teff=obj.param['teff'],
+                logg=obj.param['logg'],
+                feh=obj.param['feh']
+            )
+        elif obj.type == 'galaxy' or obj.type == 'quasar':
+            factor = 10**(-.4 * self.cosmo.distmod(obj.z).value)
+            if obj.type == 'galaxy':
+                flux = np.matmul(self.pcs, obj.coeff) * factor
+                #  if np.any(flux < 0):
+                #     raise ValueError("Glaxy %s: negative SED fluxes"%obj.id)
+                flux[flux < 0] = 0.
+                sedcat = np.vstack((self.lamb_gal, flux)).T
+                sed_data = getObservedSED(
+                    sedCat=sedcat,
+                    redshift=obj.z,
+                    av=obj.param["av"],
+                    redden=obj.param["redden"]
+                )
+                wave, flux = sed_data[0], sed_data[1]
+            elif obj.type == 'quasar':
+                flux = self.SED_AGN[int(obj.id)] * 1e-17
+                # if np.any(flux < 0):
+                #     raise ValueError("Glaxy %s: negative SED fluxes"%obj.id)
+                flux[flux < 0] = 0.
+                # sedcat = np.vstack((self.lamb_AGN, flux)).T
+                wave = self.lamb_AGN
+            # print("sed (erg/s/cm2/A) = ", sed_data)
+            # np.savetxt(os.path.join(self.config["work_dir"], "%s_sed.txt"%(obj.id)), sedcat)
+        else:
+            raise ValueError("Object type not known")
+        speci = interpolate.interp1d(wave, flux)
+        lamb = np.arange(2000, 11001+0.5, 0.5)
+        y = speci(lamb)
+        # erg/s/cm2/A --> photon/s/m2/A
+        all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13
+        sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX'))
+        
+        if obj.type == 'quasar':
+            # integrate to get the magnitudes
+            sed_photon = np.array([sed['WAVELENGTH'], sed['FLUX']]).T
+            sed_photon = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(sed_photon[:, 1]), interpolant='nearest')
+            sed_photon = galsim.SED(sed_photon, wave_type='A', flux_type='1', fast=False)
+            interFlux = integrate_sed_bandpass(sed=sed_photon, bandpass=self.filt.bandpass_full)
+            obj.param['mag_use_normal'] = getABMAG(interFlux, self.filt.bandpass_full)
+            # if obj.param['mag_use_normal'] >= 30:
+            #     print("obj ID = %d"%obj.id)
+            #     print("mag_use_normal = %.3f"%obj.param['mag_use_normal'])
+            #     print("integrated flux = %.7f"%(interFlux))
+            #     print("app mag = %.3f"%obj.param['appMag'])
+            #     np.savetxt('./AGN_SED_test/sed_objID_%d.txt'%obj.id, np.transpose([self.lamb_AGN, self.SED_AGN[int(obj.id)]]))
+            # print("obj ID = %d"%obj.id)
+            # print("mag_use_normal = %.3f"%obj.param['mag_use_normal'])
+            # print("integrated flux = %.7f"%(interFlux))
+            # print("app mag = %.3f"%obj.param['appMag'])
+            # print("abs mag = %.3f"%obj.param['absMag'])
+            # mag = getABMAG(interFlux, self.filt.bandpass_full)
+            # print("mag diff = %.3f"%(mag - obj.param['mag_use_normal']))
+        del wave
+        del flux
+        return sed

Catalog/testCat_galaxy.py

@@ -81,7 +81,7 @@ class Catalog(CatalogBase):
             self._load_SED_lib_AGN()
 
         if "rotateEll" in config["catalog_options"]:
-            self.rotation = float(int(config["catalog_options"]["rotateEll"]/45.))
+            self.rotation = np.radians(float(config["catalog_options"]["rotateEll"]))
         else:
             self.rotation = 0.
 
@@ -253,7 +253,7 @@ class Catalog(CatalogBase):
             param['id'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(igals)
             
             if param['star'] == 0:
-                obj = Galaxy(param, self.rotation, logger=self.logger)
+                obj = Galaxy(param, logger=self.logger)
             
             # Need to deal with additional output columns
             obj.additional_output_str = self.add_fmt%("n", 0., 0., 0.,

Catalog/wcs_test_C6.py

+import os
+import galsim
+import random
+import numpy as np
+import h5py as h5
+import healpy as hp
+import astropy.constants as cons
+import traceback
+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 tag_sed, getObservedSED, getABMAG, integrate_sed_bandpass, comoving_dist
+from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position
+
+# (TEST)
+from astropy.cosmology import FlatLambdaCDM
+from astropy import constants
+from astropy import units as U
+
+try:
+    import importlib.resources as pkg_resources
+except ImportError:
+    # Try backported to PY<37 'importlib_resources'
+    import importlib_resources as pkg_resources
+
+NSIDE = 128
+
+def get_bundleIndex(healpixID_ring, bundleOrder=4, healpixOrder=7):
+    assert NSIDE == 2**healpixOrder
+    shift = healpixOrder - bundleOrder
+    shift = 2*shift
+
+    nside_bundle = 2**bundleOrder
+    nside_healpix= 2**healpixOrder
+
+    healpixID_nest= hp.ring2nest(nside_healpix, healpixID_ring)
+    bundleID_nest = (healpixID_nest >> shift)
+    bundleID_ring = hp.nest2ring(nside_bundle, bundleID_nest)
+
+    return bundleID_ring
+
+class Catalog(CatalogBase):
+    def __init__(self, config, chip, pointing, chip_output, filt, **kwargs):
+        super().__init__()
+        self.cat_dir = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["cat_dir"])
+        self.seed_Av = config["catalog_options"]["seed_Av"]
+
+        self.cosmo = FlatLambdaCDM(H0=67.66, Om0=0.3111)
+
+        self.chip_output = chip_output
+        self.filt = filt
+        self.logger = chip_output.logger
+
+        with pkg_resources.path('Catalog.data', 'SLOAN_SDSS.g.fits') as filter_path:
+            self.normF_star = Table.read(str(filter_path))
+        
+        self.config = config
+        self.chip = chip
+        self.pointing = pointing
+
+        self.max_size = 0.
+
+        if "star_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["star_cat"] and not config["catalog_options"]["galaxy_only"]:
+            star_file = config["catalog_options"]["input_path"]["star_cat"]
+            star_SED_file = config["catalog_options"]["SED_templates_path"]["star_SED"]
+            self.star_path = os.path.join(self.cat_dir, star_file)
+            self.star_SED_path = os.path.join(config["data_dir"], star_SED_file)
+            self._load_SED_lib_star()
+        if "galaxy_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["galaxy_cat"] and not config["catalog_options"]["star_only"]:
+            galaxy_dir = config["catalog_options"]["input_path"]["galaxy_cat"]
+            self.galaxy_path = os.path.join(self.cat_dir, galaxy_dir)
+            self.galaxy_SED_path = os.path.join(config["data_dir"], config["catalog_options"]["SED_templates_path"]["galaxy_SED"])
+            self._load_SED_lib_gals()
+
+        if "AGN_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["AGN_cat"] and not config["catalog_options"]["star_only"]:
+            AGN_dir = config["catalog_options"]["input_path"]["AGN_cat"]
+            self.AGN_path = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["AGN_cat"])
+            self.AGN_SED_path = os.path.join(config["data_dir"], config["catalog_options"]["SED_templates_path"]["AGN_SED"])
+            self.AGN_SED_wave_path = os.path.join(config['data_dir'], config["catalog_options"]["SED_templates_path"]["AGN_SED_WAVE"])
+            self._load_SED_lib_AGN()
+
+        if "rotateEll" in config["catalog_options"]:
+            self.rotation = np.radians(float(config["catalog_options"]["rotateEll"]))
+        else:
+            self.rotation = 0.
+
+        # Update output .cat header with catalog specific output columns
+        self._add_output_columns_header()
+
+        self._get_healpix_list()
+        self._load()
+    
+    def _add_output_columns_header(self):
+        self.add_hdr = " model_tag teff logg feh"
+        self.add_hdr += " bulgemass diskmass detA e1 e2 kappa g1 g2 size galType veldisp "
+
+        self.add_fmt = " %10s %8.4f %8.4f %8.4f"
+        self.add_fmt += " %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %4d %8.4f "
+        self.chip_output.update_ouptut_header(additional_column_names=self.add_hdr)
+
+    def _get_healpix_list(self):
+        self.sky_coverage = self.chip.getSkyCoverageEnlarged(self.chip.img.wcs, margin=0.2)
+        ra_min, ra_max, dec_min, dec_max = self.sky_coverage.xmin, self.sky_coverage.xmax, self.sky_coverage.ymin, self.sky_coverage.ymax
+        ra = np.deg2rad(np.array([ra_min, ra_max, ra_max, ra_min]))
+        dec = np.deg2rad(np.array([dec_max, dec_max, dec_min, dec_min]))
+        # vertices = spherical_to_cartesian(1., dec, ra)
+        self.pix_list = hp.query_polygon(
+            NSIDE,
+            hp.ang2vec(np.radians(90.) - dec, ra),
+            inclusive=True
+        )
+        if self.logger is not None:
+            msg = str(("HEALPix List: ", self.pix_list))
+            self.logger.info(msg)
+        else:
+            print("HEALPix List: ", self.pix_list)
+
+    def load_norm_filt(self, obj):
+        if obj.type == "star":
+            return self.normF_star
+        elif obj.type == "galaxy" or obj.type == "quasar":
+            # return self.normF_galaxy
+            return None
+        else:
+            return None
+
+    def _load_SED_lib_star(self):
+        self.tempSED_star = h5.File(self.star_SED_path,'r')
+
+    def _load_SED_lib_gals(self):
+        pcs = h5.File(os.path.join(self.galaxy_SED_path, "pcs.h5"), "r")
+        lamb = h5.File(os.path.join(self.galaxy_SED_path, "lamb.h5"), "r")
+        self.lamb_gal = lamb['lamb'][()]
+        self.pcs = pcs['pcs'][()]
+
+    def _load_SED_lib_AGN(self):
+        from astropy.io import fits
+        self.SED_AGN = fits.open(self.AGN_SED_path)[0].data
+        self.lamb_AGN = np.load(self.AGN_SED_wave_path)
+
+
+    def _load_gals(self, gals, pix_id=None, cat_id=0):
+        ngals = len(gals['ra'])
+
+        # Apply astrometric modeling
+        # in C3 case only aberration
+        ra_arr = gals['ra'][:]
+        dec_arr = gals['dec'][:]
+        if self.config["obs_setting"]["enable_astrometric_model"]:
+            ra_list = ra_arr.tolist()
+            dec_list = dec_arr.tolist()
+            pmra_list = np.zeros(ngals).tolist()
+            pmdec_list = np.zeros(ngals).tolist()
+            rv_list = np.zeros(ngals).tolist()
+            parallax_list = [1e-9] * ngals
+            dt = datetime.utcfromtimestamp(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=ngals,
+                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="J2000",
+                input_date_str=date_str,
+                input_time_str=time_str
+            )
+
+        for igals in range(ngals):
+            # # (TEST)
+            # if igals > 100:
+            #     break
+            
+            param = self.initialize_param()
+            param['ra'] = ra_arr[igals]
+            param['dec'] = dec_arr[igals]
+            param['ra_orig'] = gals['ra'][igals]
+            param['dec_orig'] = gals['dec'][igals]
+            param['mag_use_normal'] = gals['mag_csst_%s'%(self.filt.filter_type)][igals]
+            # param['mag_use_normal'] = 20.
+            if self.filt.is_too_dim(mag=param['mag_use_normal'], margin=self.config["obs_setting"]["mag_lim_margin"]):
+                continue
+
+            param['z'] = gals['redshift'][igals]
+            param['model_tag'] = 'None'
+            # param['g1'] = gals['shear'][igals][0]
+            # param['g2'] = gals['shear'][igals][1]
+            param['g1'] = 0.
+            param['g2'] = 0.
+            param['kappa'] = gals['kappa'][igals]
+            # param['e1'] = gals['ellipticity_true'][igals][0]
+            # param['e2'] = gals['ellipticity_true'][igals][1]
+            param['e1'] = 0.
+            param['e2'] = 0.
+            
+            # For shape calculation
+            
+            param['ell_total'] = np.sqrt(param['e1']**2 + param['e2']**2)
+            if param['ell_total'] > 0.9:
+                continue
+            param['e1_disk'] = param['e1']
+            param['e2_disk'] = param['e2']
+            param['e1_bulge'] = param['e1']
+            param['e2_bulge'] = param['e2']
+
+
+            param['delta_ra'] = 0
+            param['delta_dec'] = 0
+
+            # Masses
+            param['bulgemass'] = gals['bulgemass'][igals]
+            param['diskmass'] = gals['diskmass'][igals]
+
+            param['size'] = gals['size'][igals]
+            if param['size'] > self.max_size:
+                self.max_size = param['size']
+
+            # Sersic index
+            param['disk_sersic_idx'] = 1.
+            param['bulge_sersic_idx'] = 4.
+
+            # Sizes
+            param['bfrac'] = param['bulgemass']/(param['bulgemass'] + param['diskmass'])
+            if param['bfrac'] >= 0.6:
+                param['hlr_bulge'] = param['size']
+                param['hlr_disk'] = param['size'] * (1. - param['bfrac'])
+            else:
+                param['hlr_disk'] = param['size']
+                param['hlr_bulge'] = param['size'] * param['bfrac']
+
+            # SED coefficients
+            param['coeff'] = gals['coeff'][igals]
+            param['detA'] = gals['detA'][igals]
+
+            # Others
+            param['galType'] = gals['type'][igals]
+            param['veldisp'] = gals['veldisp'][igals]
+            
+            # TEST no redening and no extinction
+            param['av'] = 0.0
+            param['redden'] = 0
+
+            param['star'] = 0   # Galaxy
+
+            # NOTE: this cut cannot be put before the SED type has been assigned
+            if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200):
+                continue
+
+            # TEMP
+            self.ids += 1
+            # param['id'] = self.ids
+            param['id'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(igals)
+            
+            if param['star'] == 0:
+                obj = Galaxy(param, logger=self.logger)
+            
+            # Need to deal with additional output columns
+            obj.additional_output_str = self.add_fmt%("n", 0., 0., 0.,
+                                                    param['bulgemass'], param['diskmass'], param['detA'],
+                                                    param['e1'], param['e2'], param['kappa'], param['g1'], param['g2'], param['size'],
+                                                    param['galType'], param['veldisp'])
+            
+            self.objs.append(obj)
+
+    def _load_stars(self, stars, pix_id=None):
+        nstars = len(stars['sourceID'])
+        # Apply astrometric modeling
+        ra_arr = stars["RA"][:]
+        dec_arr = stars["Dec"][:]
+        pmra_arr = stars['pmra'][:]
+        pmdec_arr = stars['pmdec'][:]
+        rv_arr = stars['RV'][:]
+        parallax_arr = stars['parallax'][:]
+        if self.config["obs_setting"]["enable_astrometric_model"]:
+            ra_list = ra_arr.tolist()
+            dec_list = dec_arr.tolist()
+            pmra_list = pmra_arr.tolist()
+            pmdec_list = pmdec_arr.tolist()
+            rv_list = rv_arr.tolist()
+            parallax_list = parallax_arr.tolist()
+            dt = datetime.utcfromtimestamp(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="J2000",
+                input_date_str=date_str,
+                input_time_str=time_str
+            )
+        for istars in range(nstars):
+            # # (TEST)
+            # if istars > 100:
+            #     break
+
+            param = self.initialize_param()
+            param['ra'] = ra_arr[istars]
+            param['dec'] = dec_arr[istars]
+            param['ra_orig'] = stars["RA"][istars]
+            param['dec_orig'] = stars["Dec"][istars]
+            param['pmra'] = pmra_arr[istars]
+            param['pmdec'] = pmdec_arr[istars]
+            param['rv'] = rv_arr[istars]
+            param['parallax'] = parallax_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]
+            param['mag_use_normal'] = 20.
+            # if param['mag_use_normal'] >= 26.5:
+            #     continue
+            self.ids += 1
+            param['id'] = stars['sourceID'][istars]
+            param['sed_type'] = stars['sourceID'][istars]
+            param['model_tag'] = stars['model_tag'][istars]
+            param['teff'] = stars['teff'][istars]
+            param['logg'] = stars['grav'][istars]
+            param['feh'] = stars['feh'][istars]
+            param['z'] = 0.0
+            param['star'] = 1   # Star
+            obj = Star(param, logger=self.logger)
+
+            # Append additional output columns to the .cat file
+            obj.additional_output_str = self.add_fmt%(param["model_tag"], param['teff'], param['logg'], param['feh'],
+                                                    0., 0., 0., 0., 0., 0., 0., 0., 0., -1, 0.)
+
+            self.objs.append(obj)
+
+    def _load_AGNs(self):
+        data = Table.read(self.AGN_path)
+        ra_arr = data['ra']
+        dec_arr = data['dec']
+        nAGNs = len(data)
+        if self.config["obs_setting"]["enable_astrometric_model"]:
+            ra_list = ra_arr.tolist()
+            dec_list = dec_arr.tolist()
+            pmra_list = np.zeros(nAGNs).tolist()
+            pmdec_list = np.zeros(nAGNs).tolist()
+            rv_list = np.zeros(nAGNs).tolist()
+            parallax_list = [1e-9] * nAGNs
+            dt = datetime.utcfromtimestamp(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=nAGNs,
+                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="J2000",
+                input_date_str=date_str,
+                input_time_str=time_str
+            )
+        for iAGNs in range(nAGNs):
+            param = self.initialize_param()
+            param['ra'] = ra_arr[iAGNs]
+            param['dec'] = dec_arr[iAGNs]
+            param['ra_orig'] = data['ra'][iAGNs]
+            param['dec_orig'] = data['dec'][iAGNs]
+            param['z'] = data['z'][iAGNs]
+            param['appMag'] = data['appMag'][iAGNs]
+            param['absMag'] = data['absMag'][iAGNs]
+            
+            # NOTE: this cut cannot be put before the SED type has been assigned
+            if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200):
+                continue
+
+            # TEST no redening and no extinction
+            param['av'] = 0.0
+            param['redden'] = 0
+
+            param['star'] = 2   # Quasar
+            param['id'] = data['igmlos'][iAGNs]
+
+            if param['star'] == 2:
+                obj = Quasar(param, logger=self.logger)
+
+            # Append additional output columns to the .cat file
+            obj.additional_output_str = self.add_fmt%("n", 0., 0., 0.,
+                                                    0., 0., 0., 0., 0., 0., 0., 0., 0., -1, 0.)
+            self.objs.append(obj)
+
+    def _load(self, **kwargs):
+        self.objs = []
+        self.ids = 0
+        
+        if "star_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["star_cat"] and not self.config["catalog_options"]["galaxy_only"]:
+            star_cat = h5.File(self.star_path, 'r')['catalog']
+            for pix in self.pix_list:
+                try:
+                    stars = star_cat[str(pix)]
+                    self._load_stars(stars, pix_id=pix)
+                    del stars
+                except Exception as e:
+                    self.logger.error(str(e))
+                    print(e)
+        
+        if "galaxy_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["galaxy_cat"] and not self.config["catalog_options"]["star_only"]:
+            for pix in self.pix_list:
+                try:
+                    bundleID  = get_bundleIndex(pix)
+                    file_path = os.path.join(self.galaxy_path, "galaxies_C6_bundle{:06}.h5".format(bundleID))
+                    gals_cat = h5.File(file_path, 'r')['galaxies']
+                    gals = gals_cat[str(pix)]
+                    self._load_gals(gals, pix_id=pix, cat_id=bundleID)
+                    del gals
+                except Exception as e:
+                    traceback.print_exc()
+                    self.logger.error(str(e))
+                    print(e)
+
+        if "AGN_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["AGN_cat"] and not self.config["catalog_options"]["star_only"]:
+            try:
+                self._load_AGNs()
+            except Exception as e:
+                traceback.print_exc()
+                self.logger.error(str(e))
+                print(e)
+
+        if self.logger is not None:
+            self.logger.info("maximum galaxy size: %.4f"%(self.max_size))
+            self.logger.info("number of objects in catalog: %d"%(len(self.objs)))
+        else:
+            print("number of objects in catalog: ", len(self.objs))
+
+    def load_sed(self, obj, **kwargs):
+        if obj.type == 'star':
+            _, wave, flux = tag_sed(
+                h5file=self.tempSED_star,
+                model_tag=obj.param['model_tag'],
+                teff=obj.param['teff'],
+                logg=obj.param['logg'],
+                feh=obj.param['feh']
+            )
+        elif obj.type == 'galaxy' or obj.type == 'quasar':
+            factor = 10**(-.4 * self.cosmo.distmod(obj.z).value)
+            if obj.type == 'galaxy':
+                flux = np.matmul(self.pcs, obj.coeff) * factor
+                #  if np.any(flux < 0):
+                #     raise ValueError("Glaxy %s: negative SED fluxes"%obj.id)
+                flux[flux < 0] = 0.
+                sedcat = np.vstack((self.lamb_gal, flux)).T
+                sed_data = getObservedSED(
+                    sedCat=sedcat,
+                    redshift=obj.z,
+                    av=obj.param["av"],
+                    redden=obj.param["redden"]
+                )
+                wave, flux = sed_data[0], sed_data[1]
+            elif obj.type == 'quasar':
+                flux = self.SED_AGN[int(obj.id)] * 1e-17
+                # if np.any(flux < 0):
+                #     raise ValueError("Glaxy %s: negative SED fluxes"%obj.id)
+                flux[flux < 0] = 0.
+                # sedcat = np.vstack((self.lamb_AGN, flux)).T
+                wave = self.lamb_AGN
+            # print("sed (erg/s/cm2/A) = ", sed_data)
+            # np.savetxt(os.path.join(self.config["work_dir"], "%s_sed.txt"%(obj.id)), sedcat)
+        else:
+            raise ValueError("Object type not known")
+        speci = interpolate.interp1d(wave, flux)
+        lamb = np.arange(2000, 11001+0.5, 0.5)
+        y = speci(lamb)
+        # erg/s/cm2/A --> photon/s/m2/A
+        all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13
+        sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX'))
+        
+        if obj.type == 'quasar':
+            # integrate to get the magnitudes
+            sed_photon = np.array([sed['WAVELENGTH'], sed['FLUX']]).T
+            sed_photon = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(sed_photon[:, 1]), interpolant='nearest')
+            sed_photon = galsim.SED(sed_photon, wave_type='A', flux_type='1', fast=False)
+            interFlux = integrate_sed_bandpass(sed=sed_photon, bandpass=self.filt.bandpass_full)
+            obj.param['mag_use_normal'] = getABMAG(interFlux, self.filt.bandpass_full)
+            # if obj.param['mag_use_normal'] >= 30:
+            #     print("obj ID = %d"%obj.id)
+            #     print("mag_use_normal = %.3f"%obj.param['mag_use_normal'])
+            #     print("integrated flux = %.7f"%(interFlux))
+            #     print("app mag = %.3f"%obj.param['appMag'])
+            #     np.savetxt('./AGN_SED_test/sed_objID_%d.txt'%obj.id, np.transpose([self.lamb_AGN, self.SED_AGN[int(obj.id)]]))
+            # print("obj ID = %d"%obj.id)
+            # print("mag_use_normal = %.3f"%obj.param['mag_use_normal'])
+            # print("integrated flux = %.7f"%(interFlux))
+            # print("app mag = %.3f"%obj.param['appMag'])
+            # print("abs mag = %.3f"%obj.param['absMag'])
+            # mag = getABMAG(interFlux, self.filt.bandpass_full)
+            # print("mag diff = %.3f"%(mag - obj.param['mag_use_normal']))
+        del wave
+        del flux
+        return sed

ObservationSim/Config/Header/ImageHeader.py

@@ -18,6 +18,10 @@ from astropy.wcs.utils import fit_wcs_from_points
 from astropy.time import Time
 from astropy import wcs
 
+from datetime import datetime
+# import socket
+import platform
+
 def chara2digit(char):
     """ Function to judge and convert characters to digitals
 
@@ -287,7 +291,8 @@ def WCS_def(xlen = 9216, ylen = 9232, gapy = 898.0, gapx1 = 534, gapx2 = 1309, r
         r_dat['CRPIX1'] = -x1
         r_dat['CRPIX2'] = -y1
 
-        cd = np.array([[ pixel_scale,  0], [0, pixel_scale]])/3600.*flag_x[1]
+        # cd = np.array([[ pixel_scale,  0], [0, pixel_scale]])/3600.*flag_x[1]
+        cd = np.array([[ pixel_scale,  0], [0, -pixel_scale]])/3600.
         cd_rot = rotate_CD_matrix(cd, pa_aper)
         r_dat['CD1_1'] = cd_rot[0,0]
         r_dat['CD1_2'] = cd_rot[0,1]
@@ -381,7 +386,7 @@ def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec
     h_prim['EXPTIME'] = exptime
     
     # Define file types
-    file_type = {'SCI':'SCIE', 'BIAS':'BIAS', 'DARK':'DARK', 'FLAT':'FLAT', 'CRS':'cosmic_ray', 'CRD':'cosmic_ray','CALS':'CALS','CALF':'CALF'}
+    file_type = {'SCI':'SCIE', 'BIAS':'BIAS', 'DARK':'DARK', 'FLAT':'FLAT', 'CRS':'CRS', 'CRD':'CRD','CALS':'CALS','CALF':'CALF'}
     h_prim['FILETYPE'] = file_type[im_type]
 
     co = coord.SkyCoord(ra, dec, unit='deg')
@@ -417,7 +422,7 @@ def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec
     file_end_time = end_time_str[0:4] + end_time_str[5:7]+end_time_str[8:10] + end_time_str[11:13] + end_time_str[14:16] + end_time_str[17:19]
     # h_prim['FILENAME'] = 'CSST_MSC_MS_' + im_type + '_' + file_start_time + '_' + file_end_time + '_' + OBS_id + '_' + CCDID[
     #     k - 1].rjust(2, '0') + '_L0_V01'
-    h_prim['FILENAME'] = 'CSST_MSC_MS_' + im_type + '_' + file_start_time + '_' + file_end_time + '_' + OBS_id + '_' + chip_name + '_L0_V01'
+    h_prim['FILENAME'] = 'CSST_MSC_MS_' + file_type[im_type] + '_' + file_start_time + '_' + file_end_time + '_' + OBS_id + '_' + chip_name + '_L0_V01'
 
 
     h_prim['POSI0_X'] = sat_pos[0]
@@ -433,12 +438,15 @@ def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec
     # Get version of CSSTSim Package
     from pkg_resources import get_distribution
     # h_prim['SIM_VER'] = (get_distribution("CSSTSim").version, "Version of CSST MSC simulation software")
-    h_prim['FITSCREA'] = get_distribution("CSSTSim").version
+    currentDateAndTime = datetime.now()
+    compute_name = platform.node()
+    h_prim['FITSCREA'] = get_distribution("CSSTSim").version +'_' + currentDateAndTime.strftime("%Y%m%d") + '_' +compute_name
+    h_prim['EPOCH'] = round((Time(h_prim['EXPSTART'], format='mjd', scale='tcb')).jyear, 1)
 
     return h_prim
 
 def generateExtensionHeader(chip, xlen = 9216, ylen = 9232,ra = 60, dec = -40, pa = -23.433, gain = 1.0, readout = 5.0, dark = 0.02, saturation=90000, pixel_scale = 0.074, pixel_size=1e-2, 
-                            extName='SCI', row_num = None, col_num = None, xcen=None, ycen=None):
+                            extName='SCI', row_num = None, col_num = None, xcen=None, ycen=None, timestamp = 1621915200,exptime = 150., readoutTime = 40.):
 
     e_header_fn = os.path.split(os.path.realpath(__file__))[0] + '/extension_header.header'
     f = open(os.path.split(os.path.realpath(__file__))[0] + '/filter.lst')
@@ -460,47 +468,60 @@ def generateExtensionHeader(chip, xlen = 9216, ylen = 9232,ra = 60, dec = -40, p
     # h_ext['CCDCHIP'] = CCDID[k - 1].rjust(2, '0')
     # h_ext['CCDLABEL'] = filters[k-1] + '-' + filterID[k-1]
     # h_ext['FILTER'] = filters[k-1]
-    h_ext['CCDCHIP'] = chip.chip_name
-    h_ext['CCDLABEL'] = chip.filter_type + '-' + str(chip.chipID).rjust(2, '0')
+    h_ext['CHIPID'] = str(chip.chipID).rjust(2, '0')
+    h_ext['CHIPLAB'] = chip.chip_name
     h_ext['FILTER'] = chip.filter_type
     h_ext['NAXIS1'] = xlen
     h_ext['NAXIS2'] = ylen
     h_ext['EXTNAME'] = extName
-    h_ext['GAIN1'] = gain
-    h_ext['GAIN2'] = gain
-    h_ext['GAIN3'] = gain
-    h_ext['GAIN4'] = gain
-    h_ext['GAIN5'] = gain
-    h_ext['GAIN6'] = gain
-    h_ext['GAIN7'] = gain
-    h_ext['GAIN8'] = gain
-    h_ext['GAIN9'] = gain
-    h_ext['GAIN10'] = gain
-    h_ext['GAIN11'] = gain
-    h_ext['GAIN12'] = gain
-    h_ext['GAIN13'] = gain
-    h_ext['GAIN14'] = gain
-    h_ext['GAIN15'] = gain
-    h_ext['GAIN16'] = gain
-    h_ext['RDNOIS1'] = readout
-    h_ext['RDNOIS2'] = readout
-    h_ext['RDNOIS3'] = readout
-    h_ext['RDNOIS4'] = readout
-    h_ext['RDNOIS5'] = readout
-    h_ext['RDNOIS6'] = readout
-    h_ext['RDNOIS7'] = readout
-    h_ext['RDNOIS8'] = readout
-    h_ext['RDNOIS9'] = readout
-    h_ext['RDNOIS10'] = readout
-    h_ext['RDNOIS11'] = readout
-    h_ext['RDNOIS12'] = readout
-    h_ext['RDNOIS13'] = readout
-    h_ext['RDNOIS14'] = readout
-    h_ext['RDNOIS15'] = readout
-    h_ext['RDNOIS16'] = readout
+    h_ext['GAIN01'] = chip.gain_channel[0]
+    h_ext['GAIN02'] = chip.gain_channel[1]
+    h_ext['GAIN03'] = chip.gain_channel[2]
+    h_ext['GAIN04'] = chip.gain_channel[3]
+    h_ext['GAIN05'] = chip.gain_channel[4]
+    h_ext['GAIN06'] = chip.gain_channel[5]
+    h_ext['GAIN07'] = chip.gain_channel[6]
+    h_ext['GAIN08'] = chip.gain_channel[7]
+    h_ext['GAIN09'] = chip.gain_channel[8]
+    h_ext['GAIN10'] = chip.gain_channel[9]
+    h_ext['GAIN11'] = chip.gain_channel[10]
+    h_ext['GAIN12'] = chip.gain_channel[11]
+    h_ext['GAIN13'] = chip.gain_channel[12]
+    h_ext['GAIN14'] = chip.gain_channel[13]
+    h_ext['GAIN15'] = chip.gain_channel[14]
+    h_ext['GAIN16'] = chip.gain_channel[15]
+    h_ext['RON01'] = readout
+    h_ext['RON02'] = readout
+    h_ext['RON03'] = readout
+    h_ext['RON04'] = readout
+    h_ext['RON05'] = readout
+    h_ext['RON06'] = readout
+    h_ext['RON07'] = readout
+    h_ext['RON08'] = readout
+    h_ext['RON09'] = readout
+    h_ext['RON10'] = readout
+    h_ext['RON11'] = readout
+    h_ext['RON12'] = readout
+    h_ext['RON13'] = readout
+    h_ext['RON14'] = readout
+    h_ext['RON15'] = readout
+    h_ext['RON16'] = readout
 
     h_ext['PIXSCAL1'] = pixel_scale
     h_ext['PIXSCAL2'] = pixel_scale
+    h_ext['EXPTIME'] = exptime
+    h_ext['DARKTIME'] = exptime + 2
+
+    datetime_obs = datetime.utcfromtimestamp(timestamp)
+    tstart = Time(datetime_obs)
+    tstart_read = Time(tstart.mjd + exptime / 86400., format="mjd")
+    tend_read = Time(tstart.mjd + (exptime + readoutTime) / 86400., format="mjd")
+    t_s1 = str(tstart_read.datetime).split()
+    h_ext['READT0'] = t_s1[0]+'T'+t_s1[1]
+    t_s2 = str(tend_read.datetime).split()
+    h_ext['READT1'] = t_s2[0] + 'T' + t_s2[1]
+
+
 
     # h_ext['POS_ANG'] = pa
     header_wcs = WCS_def(xlen=xlen, ylen=ylen, gapy=898.0, gapx1=534, gapx2=1309, ra_ref=ra, dec_ref=dec, pa=pa, pixel_scale=pixel_scale, pixel_size=pixel_size, 
@@ -514,8 +535,8 @@ def generateExtensionHeader(chip, xlen = 9216, ylen = 9232,ra = 60, dec = -40, p
     h_ext['CD1_2'] = header_wcs['CD1_2']
     h_ext['CD2_1'] = header_wcs['CD2_1']
     h_ext['CD2_2'] = header_wcs['CD2_2']
-    h_ext['EQUINOX'] = header_wcs['EQUINOX']
-    h_ext['WCSDIM'] = header_wcs['WCSDIM']
+    # h_ext['EQUINOX'] = header_wcs['EQUINOX']
+    # h_ext['WCSDIM'] = header_wcs['WCSDIM']
     h_ext['CTYPE1'] = header_wcs['CTYPE1']
     h_ext['CTYPE2'] = header_wcs['CTYPE2']
 

ObservationSim/Config/Header/extension_header.header

-XTENSION= 'IMAGE   '           / extension type                                 BITPIX  =                   16 / bits per data value                            NAXIS   =                    2 / number of data axes                            NAXIS1  =                 9216 / length of first array axis                     NAXIS2  =                 9232 / length of second array axis                    PCOUNT  =                    0                                                  GCOUNT  =                    1                                                  EXTNAME = 'SCI     '                                                            EXTVER  =                    1                                                  BSCALE  =                    1                                                  BZERO   =                32768                                                  BUNIT   = 'ADU     '            / physical unit of array values                 COMMENT ==================================================================      COMMENT Detector information                                                    COMMENT ==================================================================      DETECTOR= 'CCD'                 / detector name                                 DETTEMP1=                   173 / detector temperature at EXPSTART              DETTEMP2=                   173 / detector temperature at EXPEND                DETTEMP3=                   173 / detector temperature when readout is finished DETSIZE = '9560x9264'           / detector size                                 DATASEC = '9216x9232'           / data size                                     PIXSCAL1=                 0.074 / pixel scale for axis 1                        PIXSCAL2=                 0.074 / pixel scale for axis 2                        PIXSIZE1=                    10 / pixel size  for axis 1                        PIXSIZE2=                    10 / pixel size  for axis 1                        COMMENT ==================================================================      COMMENT CCD chip information                                                    COMMENT ==================================================================      CCDCHIP =                     8 / CCD chip ID                                   CCDLABEL= 'y-1'                 / CCD chip label                                FILTER  = 'y'                   / filter name                                   NCHAN   =                    16 / number of readout channels                    NCHAN1  =                     8 / number of horizintal channels                 NCHAN2  =                     2 / number of vertical channels                   PSCAN1  =                    27 / horizontal prescan width, per readout channel PSCAN2  =                     8 / horizontal prescan width, per readout channel OSCAN1  =                    16 / horizontal overscan width,per readout channel OSCAN2  =                    16 / horizontal overscan width,per readout channel COMMENT ==================================================================      COMMENT WORLD COORDINATE SYSTEM AND RELATED PARAMETERS                          COMMENT ==================================================================      WCSAXES =                    2 / number of World Coordinate System axes         CRPIX1  =             -10017.0 / x-coordinate of reference pixel                CRPIX2  =              24876.0 / y-coordinate of reference pixel                CRVAL1  =            62.228226 / first axis value at reference pixel            CRVAL2  =           -42.316932 / second axis value at reference pixel           CTYPE1  = 'RA---TAN'           / the coordinate type for the first axis         CTYPE2  = 'DEC--TAN'           / the coordinate type for the second axis        CD1_1   = 1.88602083707394E-05 / partial of first axis coordinate w.r.t. x      CD2_1   = -8.1745583617600E-06 / partial of first axis coordinate w.r.t. x      CD1_2   = 8.17455836176000E-06 / partial of second axis coordinate w.r.t. x     CD2_2   = 1.88602083707394E-05 / partial of second axis coordinate w.r.t. x     OTHERS  =                 ''   /                                                COMMENT ==================================================================      COMMENT Readout information                                                     COMMENT ==================================================================      GAIN1   =                  1.1 / CCD gain (channel 1)                           GAIN2   =                  1.1 / CCD gain (channel 2)                           GAIN3   =                  1.1 / CCD gain (channel 3)                           GAIN4   =                  1.1 / CCD gain (channel 4)                           GAIN5   =                  1.1 / CCD gain (channel 5)                           GAIN6   =                  1.1 / CCD gain (channel 6)                           GAIN7   =                  1.1 / CCD gain (channel 7)                           GAIN8   =                  1.1 / CCD gain (channel 8)                           GAIN9   =                  1.1 / CCD gain (channel 9)                           GAIN10  =                  1.1 / CCD gain (channel 10)                          GAIN11  =                  1.1 / CCD gain (channel 11)                          GAIN12  =                  1.1 / CCD gain (channel 12)                          GAIN13  =                  1.1 / CCD gain (channel 13)                          GAIN14  =                  1.1 / CCD gain (channel 14)                          GAIN15  =                  1.1 / CCD gain (channel 15)                          GAIN16  =                  1.1 / CCD gain (channel 16)                          RDNOIS1 =                  5.0 / read noise  (channel 1)                        RDNOIS2 =                  5.0 / read noise  (channel 2)                        RDNOIS3 =                  5.0 / read noise  (channel 3)                        RDNOIS4 =                  5.0 / read noise  (channel 4)                        RDNOIS5 =                  5.0 / read noise  (channel 5)                        RDNOIS6 =                  5.0 / read noise  (channel 6)                        RDNOIS7 =                  5.0 / read noise  (channel 7)                        RDNOIS8 =                  5.0 / read noise  (channel 8)                        RDNOIS9 =                  5.0 / read noise  (channel 9)                        RDNOIS10=                  5.0 / read noise  (channel 10)                       RDNOIS11=                  5.0 / read noise  (channel 11)                       RDNOIS12=                  5.0 / read noise  (channel 12)                       RDNOIS13=                  5.0 / read noise  (channel 13)                       RDNOIS14=                  5.0 / read noise  (channel 14)                       RDNOIS15=                  5.0 / read noise  (channel 15)                       RDNOIS16=                  5.0 / read noise  (channel 16)                       RDSPEED =                 10.0 / read speed (in MHz)                            COMMENT ==================================================================      COMMENT Shutter information                                                     COMMENT ==================================================================      SHUTSTAT=                    T / shutter status                                 SHTOPEN0=                  0.0 / shutter open time (begin)                      SHTOPEN1=                  0.0 / shutter open time (end)                        SHTCLOS0=                  0.0 / shutter close time (begin)                     SHTCLOS1=                  0.0 / shutter close time (end)                       COMMENT ==================================================================      COMMENT LED information                                                         COMMENT ==================================================================      LEDSTAT1= '0000000000000000000000000000' / main LED status                      LEDEXP1 =                  0.0 / main LED flash time (s)                        LEDSTAT2= '0000000000000000000000000000' / backup LED status                    LEDEXP2 =                  0.0 / backup LED flash time (s)                      LEDTEMPA=                173.0 / LED temperature (main LED in K)                LEDTEMPB=                173.0 / LED temperature (backup LED in K)              COMMENT ==================================================================      COMMENT Other information                                                       COMMENT ==================================================================      CHECKSUM= '''abcde'''          / SHA256 checksum of global headers              END                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          
+XTENSION= 'IMAGE   '           / extension type                                 BITPIX  =                   16 / bits per data value                            NAXIS   =                    2 / number of data axes                            NAXIS1  =                 9216 / length of first array axis                     NAXIS2  =                 9232 / length of second array axis                    PCOUNT  =                    0                                                  GCOUNT  =                    1                                                  EXTNAME = 'SCI     '                                                            EXTVER  =                    1                                                  BSCALE  =                    1                                                  BZERO   =                32768                                                  BUNIT   = 'ADU     '            / physical unit of array values                 COMMENT ==================================================================      COMMENT Detector information                                                    COMMENT ==================================================================      CAMERA  = 'MS'                  / camera of main survey                         DETSN   = '12345678'            / detector serial number                        DETNAME = 'CCD'                 / detector type                                 DETTEMP1=                 173.0 / detector temperature at EXPSTART(in Kelvin)   DETTEMP2=                 173.0 / detector temperature at EXPEND(in Kelvin)     DETTEMP3=                 173.0 / detector temperature at READT1(in Kelvin)     DETSIZE = '9560x9264'           / detector size                                 DATASECT= '9216x9232'           / data section                                  PIXSCAL1=                 0.074 / pixel scale for axis 1                        PIXSCAL2=                 0.074 / pixel scale for axis 2                        PIXSIZE1=                    10 / pixel size for axis 1 (in um)                 PIXSIZE2=                    10 / pixel size for axis 2 (in um)                 COMMENT ==================================================================      COMMENT CCD chip information                                                    COMMENT ==================================================================      CHIPID  = '08'                  / chip ID                                       CHIPLAB = 'y-1'                 / chip label                                    FILTER  = 'y'                   / filter name                                   NCHAN   =                    16 / number of readout channels                    PSCAN1  =                    27 / horizontal prescan width, per readout channel PSCAN2  =                     8 / vertical prescan width, per readout channel   OSCAN1  =                    16 / horizontal overscan width,per readout channel OSCAN2  =                    16 / vertical overscan width,per readout channel   COMMENT ==================================================================      COMMENT WORLD COORDINATE SYSTEM AND RELATED PARAMETERS                          COMMENT ==================================================================      WCSAXES =                    2 / number of World Coordinate System axes         CRPIX1  =             -10017.0 / x-coordinate of reference pixel                CRPIX2  =              24876.0 / y-coordinate of reference pixel                CRVAL1  =            62.228226 / first axis value at reference pixel            CRVAL2  =           -42.316932 / second axis value at reference pixel           CTYPE1  = 'RA---TAN'           / the coordinate type for the first axis         CTYPE2  = 'DEC--TAN'           / the coordinate type for the second axis        CD1_1   = 1.88602083707394E-05 / partial of first axis coordinate w.r.t.x       CD1_2   = 8.17455836176000E-06 / partial of first axis coordinate w.r.t.y       CD2_1   = -8.1745583617600E-06 / partial of second axis coordinate w.r.t.x      CD2_2   = 1.88602083707394E-05 / partial of second axis coordinate w.r.t.y      OTHERS  =                 ''   /                                                COMMENT ==================================================================      COMMENT Readout information                                                     COMMENT ==================================================================      GAINLVL = '01'                 / gain level                                     GAIN01  =                  1.1 / gain (channel 01)                              GAIN02  =                  1.1 / gain (channel 02)                              GAIN03  =                  1.1 / gain (channel 03)                              GAIN04  =                  1.1 / gain (channel 04)                              GAIN05  =                  1.1 / gain (channel 05)                              GAIN06  =                  1.1 / gain (channel 06)                              GAIN07  =                  1.1 / gain (channel 07)                              GAIN08  =                  1.1 / gain (channel 08)                              GAIN09  =                  1.1 / gain (channel 09)                              GAIN10  =                  1.1 / gain (channel 10)                              GAIN11  =                  1.1 / gain (channel 11)                              GAIN12  =                  1.1 / gain (channel 12)                              GAIN13  =                  1.1 / gain (channel 13)                              GAIN14  =                  1.1 / gain (channel 14)                              GAIN15  =                  1.1 / gain (channel 15)                              GAIN16  =                  1.1 / gain (channel 16)                              RON01   =                  5.0 / read noise (channel 01)                        RON02   =                  5.0 / read noise (channel 02)                        RON03   =                  5.0 / read noise (channel 03)                        RON04   =                  5.0 / read noise (channel 04)                        RON05   =                  5.0 / read noise (channel 05)                        RON06   =                  5.0 / read noise (channel 06)                        RON07   =                  5.0 / read noise (channel 07)                        RON08   =                  5.0 / read noise (channel 08)                        RON09   =                  5.0 / read noise (channel 09)                        RON10   =                  5.0 / read noise (channel 10)                        RON11   =                  5.0 / read noise (channel 11)                        RON12   =                  5.0 / read noise (channel 12)                        RON13   =                  5.0 / read noise (channel 13)                        RON14   =                  5.0 / read noise (channel 14)                        RON15   =                  5.0 / read noise (channel 15)                        RON16   =                  5.0 / read noise (channel 16)                        READT0  = '2024-00-00T00:00:00'/ readout start time(UTC)                        READT1  = '2024-00-00T00:00:00'/ readout end time(UTC)                          ROSPEED =                 10.0 / readout speed (in MHz)                         EXPTIME =                150.0 / exposure duration                              DARKTIME=                150.0 / dark current time                              COMMENT ==================================================================      COMMENT Shutter information                                                     COMMENT ==================================================================      SHTSTAT =                    T / shutter status                                 SHTOPEN0=                  0.0 / shutter open time (begin)                      SHTOPEN1=                  0.0 / shutter open time (end)                        SHTCLOS0=                  0.0 / shutter close time (begin)                     SHTCLOS1=                  0.0 / shutter close time (end)                       COMMENT ==================================================================      COMMENT LED information                                                         COMMENT ==================================================================      LEDFLAG =                    0 / main/backup LED                                LEDSTAT = '00000000000000'     / LED status                                     LEDEXPT =                  0.0 / LED flash time (s)                             LEDTEMP =                173.0 / LED temperature (in K)                         COMMENT ==================================================================      COMMENT Other information                                                       COMMENT ==================================================================      CHECKSUM= '''abcde'''          / HDU checksum updated yyyy-mm-ddTHH:MM:SS       DATASUM = '''abcde'''          / data unit checksum updated yyyy-mm-ddTHH:MM:SS END                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          

ObservationSim/Config/Header/global_header.header

-SIMPLE  =                    T / conforms to FITS standard                      BITPIX  =                    8 / array data type                                NAXIS   =                    0 / number of array dimensions                     EXTEND  =                    T                                                  NEXTEND =                    1                                                  GROUPS  =                    F                                                  DATE    = '2021-03-04T09:30:00'/ date this file was written                     FILENAME= 'MSC_MS_210304093000_100000000_06_raw' / file name                    FILETYPE= 'SCI     '           / observation type                               TELESCOP= 'CSST    '           / telescope used to acquire data                 INSTRUME= 'MSC     '           / instrument used to acquire data                RADECSYS= 'ICRS    '           / frame of reference of coordinates              EQUINOX =               2000.0 /                                                FITSCREA= 'C6'                 / FITS create software version                   COMMENT ==================================================================      COMMENT Object information                                                      COMMENT ==================================================================      OBJECT  = '00000000'           / object name                                    TARGET  = '+000000000000'      / target name (hhmmss+ddmmss)                    OBSID   = '00000000'           / observation ID                                 OBJ_RA  =            62.228226 / R.A. of the object (degrees)                   OBJ_DEC =           -42.316932 / declination of the object (degrees)            COMMENT ==================================================================      COMMENT Telescope information                                                   COMMENT ==================================================================      REFFRAME= 'CSSTGSC-1.0'        / guide star catalog version                     DATE-OBS= '2021-03-04T09:30:00'/ date of the observation (yyyy-mm-dd hh:mm:ss)  SATESWV = '0001'               / software version in the satellite              EXPSTART=              59130.5 / exposure start time (MJD)                      CABSTART=              59130.5 / (MJD)                                          SUNANGL0=                 50.0 / angle between sun and optical axis at CABSTART MOONANG0=                 30.0 / angle moon and opt axis at CABSTART at CABST   TEL_ALT0=                 20.0 / angle opt axis and the ground-piston at CABST  POS_ANG0=                 20.0 / angle bwt y axis and the North Pole at CABST   POSI0_X =                  0.0 / the orbital position in X at CABSTART          POSI0_Y =                  0.0 / the orbital position in Y at CABSTART          POSI0_Z =                  0.0 / the orbital position in Z at CABSTART          VELO0_X =                  0.0 / the orbital velocity in X at CABSTART          VELO0_Y =                  0.0 / the orbital velocity in Y at CABSTART          VELO0_Z =                  0.0 / the orbital velocity in Z at CABSTART          EULER0_1=                  0.0 / Euler angle 1 at  CABSTART                     EULER0_2=                  0.0 / Euler angle 2 at  CABSTART                     EULER0_3=                  0.0 / Euler angle 3 at  CABSTART                     RA_PNT0 =                  0.0 / RA of the pointing (degrees) at CABSTART       DEC_PNT0=                  0.0 / DEC of the pointing (degrees) at CABSTART      EXPEND  =                  0.0 / exposure end time (MJD)                        CABEND  =                  0.0 / (MJD)                                          SUNANGL1=                 50.0 / angle between sun and optical axis at CABEND   MOONANG1=                 30.0 / angle btw moon and optical axis at CABEND      TEL_ALT1=                 20.0 / angle opt axis and the ground-piston at CAEND  POS_ANG1=                 20.0 / angle bwt y axis and the North Pole at CABEND  POSI1_X =                  0.0 / the orbital position in X at CABEND            POSI1_Y =                  0.0 / the orbital position in Y at CABEND            POSI1_Z =                  0.0 / the orbital position in Z at CABEND            VELO1_X =                  0.0 / The orbital velocity in X at CABEND            VELO1_Y =                  0.0 / The orbital velocity in Y at CABEND            VELO1_Z =                  0.0 / The orbital velocity in Z at CABEND            EULER1_1=                  0.0 / Euler angle 1 at  CABEND                       EULER1_2=                  0.0 / Euler angle 2 at  CABEND                       EULER1_3=                  0.0 / Euler angle 3 at  CABEND                       RA_PNT1 =                  0.0 / RA of the pointing (degrees) at CABEND         DEC_PNT1=                  0.0 / DEC of the pointing (degrees) at CABEND        EXPTIME =                150.0 / exposure duration                              EPOCH   =                 2000 / coordinate epoch                               COMMENT Other information                                                       COMMENT ==================================================================      CHECKSUM= 'abcdefg '           / SHA256 checksum of global headers              END                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             
+SIMPLE  =                    T / conforms to FITS standard                      BITPIX  =                    8 / array data type                                NAXIS   =                    0 / number of array dimensions                     EXTEND  =                    T                                                  NEXTEND =                    1 / number of array dimensions                     GROUPS  =                    F / '   '                                          DATE    = '2021-03-04T09:30:00'/ the date on which this file was written        FILENAME= 'MSC_MS_210304093000_100000000_06_raw' / file name                    FILETYPE= 'SCIE    '           / observation type                               TELESCOP= 'CSST    '           / telescope used to acquire data                 INSTRUME= 'MSC     '           / instrument used to acquire data                RADECSYS= 'ICRS    '           / reference coordinates system                   EQUINOX =               2000.0 /                                                FITSCREA= 'C6'                 / FITS create software version                   COMMENT ==================================================================      COMMENT Object information                                                      COMMENT ==================================================================      OBJECT  = '00000000'           / object name                                    TARGET  = '+000000000000'      / target name (hhmmss+ddmmss)                    OBSID   = '00000000'           / observation ID                                 OBJ_RA  =            62.228226 / R.A. of the object (degrees)                   OBJ_DEC =           -42.316932 / declination of the object (degrees)            COMMENT ==================================================================      COMMENT Telescope information                                                   COMMENT ==================================================================      REFFRAME= 'CSSTGSC-1.0'        / guide star catalog version                     DATE-OBS= '2021-03-04T09:30:00'/ date of the observation (yyyy-mm-dd hh:mm:ss)  SATESWV = '0001'               / software version in the satellite              EXPSTART=              59130.5 / exposure start time (MJD)                      CABSTART=              59130.5 / (MJD)                                          SUNANGL0=                 50.0 / angle between sun and opt axis at CABSTART     MOONANG0=                 30.0 / angle between moon and opt axis at CABSTART    TEL_ALT0=                 20.0 / angle between opt axis and Elimb at CABSTART   POS_ANG0=                 20.0 / angle between y axis and NP at CABSTART        POSI0_X =                  0.0 / the orbital position in X at CABSTART          POSI0_Y =                  0.0 / the orbital position in Y at CABSTART          POSI0_Z =                  0.0 / the orbital position in Z at CABSTART          VELO0_X =                  0.0 / the orbital velocity in X at CABSTART          VELO0_Y =                  0.0 / the orbital velocity in Y at CABSTART          VELO0_Z =                  0.0 / the orbital velocity in Z at CABSTART          EULER0_1=                  0.0 / euler angle 1 at CABSTART                      EULER0_2=                  0.0 / euler angle 2 at CABSTART                      EULER0_3=                  0.0 / euler angle 3 at CABSTART                      RA_PNT0 =                  0.0 / RA of the pointing (degrees) at CABSTART       DEC_PNT0=                  0.0 / DEC of the pointing (degrees) at CABSTART      EXPEND  =                  0.0 / exposure end time (MJD)                        CABEND  =                  0.0 / (MJD)                                          SUNANGL1=                 50.0 / angle between sun and opt axis at CABEND       MOONANG1=                 30.0 / angle between moon and opt axis at CABEND      TEL_ALT1=                 20.0 / angle between opt axis and Elimb at CABEND     POS_ANG1=                 20.0 / angle between y axis and NP at CABEND          POSI1_X =                  0.0 / the orbital position in X at CABEND            POSI1_Y =                  0.0 / the orbital position in Y at CABEND            POSI1_Z =                  0.0 / the orbital position in Z at CABEND            VELO1_X =                  0.0 / the orbital velocity in X at CABEND            VELO1_Y =                  0.0 / the orbital velocity in Y at CABEND            VELO1_Z =                  0.0 / the orbital velocity in Z at CABEND            EULER1_1=                  0.0 / euler angle 1 at CABEND                        EULER1_2=                  0.0 / euler angle 2 at CABEND                        EULER1_3=                  0.0 / euler angle 3 at CABEND                        RA_PNT1 =                  0.0 / RA of the pointing (degrees) at CABEND         DEC_PNT1=                  0.0 / DEC of the pointing (degrees) at CABEND        EXPTIME =                150.0 / exposure duration                              EPOCH   =               2000.0 / coordinate epoch                               COMMENT Other information                                                       COMMENT ==================================================================      CHECKSUM= 'abcdefg '           / HDU checksum updated yyyy-mm-ddTHH:MM:SS       DATASUM = 'abcdefg '           / data unit checksum updated yyyy-mm-ddTHH:MM:SS END                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             

ObservationSim/Config/Pointing.py

@@ -3,16 +3,18 @@ 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'):
+    def __init__(self, id=0, ra=0., dec=0., img_pa=0., timestamp=1621915200, sat_x=0., sat_y=0., sat_z=0., sun_x=0., sun_y=0., sun_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.sun_x, self.sun_y, self.sun_z = sun_x, sun_y, sun_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
+        self.jdt = 0.
 
     def read_pointing_columns(self, columns, id=0, t=1621915200, pointing_type='MS'):
         self.id = id
@@ -24,10 +26,14 @@ class Pointing(object):
         if col_len > 5:
             jdt = np.double(columns[5])
             t_temp = Time(jdt, format='jd')
+            self.jdt = jdt
             self.timestamp = t_temp.unix
             self.sat_x = float(columns[6])
             self.sat_y = float(columns[7])
             self.sat_z = float(columns[8])
+            self.sun_x = float(columns[9])
+            self.sun_y = float(columns[10])
+            self.sun_z = float(columns[1])
             self.sat_vx = float(columns[15])
             self.sat_vy = float(columns[16])
             self.sat_vz = float(columns[17])

ObservationSim/Instrument/Chip/Chip.py

@@ -27,7 +27,9 @@ class Chip(FocalPlane):
         # self.npix_x = 9216
         # self.npix_y = 9232
         # self.pix_scale  = 0.074 # pixel scale
-
+        self.nsecy = 2
+        self.nsecx = 8
+        self.gain_channel = np.ones(self.nsecy* self.nsecx)
         self._set_attributes_from_config(config)
 
         self.logger = logger
@@ -75,7 +77,8 @@ class Chip(FocalPlane):
                 self.fdModel = None
             else:
                 try:
-                    with pkg_resources.files('ObservationSim.Instrument.data.field_distortion').joinpath("FieldDistModelGlobal_mainFP_v1.0.pickle") as field_distortion:
+                    # with pkg_resources.files('ObservationSim.Instrument.data.field_distortion').joinpath("FieldDistModelGlobal_mainFP_v1.0.pickle") as field_distortion:
+                    with pkg_resources.files('ObservationSim.Instrument.data.field_distortion').joinpath("FieldDistModel_v2.0.pickle") as field_distortion:
                         with open(field_distortion, "rb") as f:
                             self.fdModel = pickle.load(f)
                 except AttributeError:
@@ -115,7 +118,7 @@ class Chip(FocalPlane):
         self._getCRdata()
 
         # Define the sensor model
-        if config["ins_effects"]["bright_fatter"] == True and self.survey_type == "photometric":
+        if "bright_fatter" in config["ins_effects"] and config["ins_effects"]["bright_fatter"] == True and self.survey_type == "photometric":
             self.sensor = galsim.SiliconSensor(strength=self.df_strength, treering_func=treering_func)
         else:
             self.sensor = galsim.Sensor()
@@ -144,18 +147,18 @@ class Chip(FocalPlane):
     def _getSurveyType(self):
         if self.filter_type in ["GI", "GV", "GU"]:
             return "spectroscopic"
-        elif self.filter_type in ["nuv", "u", "g", 'r', 'i', 'z', 'y', 'FGS']:
+        elif self.filter_type in ["NUV", "u", "g", 'r', 'i', 'z', 'y', 'FGS']:
             return "photometric"
         # elif self.filter_type in ["FGS"]:
         #     return "FGS"
 
     def _getChipEffCurve(self, filter_type):
         # CCD efficiency curves
-        if filter_type in ['nuv', 'u', 'GU']: filename = 'UV0.txt'
+        if filter_type in ['NUV', 'u', 'GU']: filename = 'UV0.txt'
         if filter_type in ['g', 'r', 'GV', 'FGS']: filename = 'Astro_MB.txt' # TODO, need to switch to the right efficiency curvey for FGS CMOS
         if filter_type in ['i', 'z', 'y', 'GI']: filename = 'Basic_NIR.txt'
         # Mirror efficiency:
-        # if filter_type == 'nuv': mirror_eff = 0.54
+        # if filter_type == 'NUV': mirror_eff = 0.54
         # if filter_type == 'u': mirror_eff = 0.68
         # if filter_type in ['g', 'r', 'i', 'z', 'y']: mirror_eff = 0.8
         # if filter_type in ['GU', 'GV', 'GI']: mirror_eff = 1. # Not sure if this is right
@@ -182,7 +185,7 @@ class Chip(FocalPlane):
     def getChipFilter(self, chipID=None, filter_layout=None):
         """Return the filter index and type for a given chip #(chipID)
         """
-        filter_type_list = ["nuv","u", "g", "r", "i","z","y","GU", "GV", "GI", "FGS"]
+        filter_type_list = ["NUV","u", "g", "r", "i","z","y","GU", "GV", "GI", "FGS"]
         if filter_layout is not None:
             return filter_layout[chipID][0], filter_layout[chipID][1]
         if chipID == None:
@@ -195,7 +198,7 @@ class Chip(FocalPlane):
         if chipID in [7, 24]:          filter_type = "i"
         if chipID in [14, 17]:         filter_type = "u"
         if chipID in [9, 22]:          filter_type = "r"
-        if chipID in [12, 13, 18, 19]: filter_type = "nuv"
+        if chipID in [12, 13, 18, 19]: filter_type = "NUV"
         if chipID in [8, 23]:          filter_type = "g"
         if chipID in [1, 10, 21, 30]:  filter_type = "GI"
         if chipID in [2, 5, 26, 29]:   filter_type = "GV"
@@ -304,7 +307,10 @@ class Chip(FocalPlane):
         if self.chipID == 30: confFile = ['CSST_GI8.conf', 'CSST_GI7.conf']
         return confFile
 
-    def generateHeader(self, ra_cen, dec_cen, img_rot, im_type, pointing_ID, date_obs, time_obs, exptime=150.):
+    def generateHeader(self, ra_cen, dec_cen, img_rot, im_type, pointing_ID, exptime=150., timestamp = 1621915200):
+        datetime_obs = datetime.utcfromtimestamp(timestamp)
+        date_obs = datetime_obs.strftime("%y%m%d")
+        time_obs = datetime_obs.strftime("%H%M%S")
         h_prim = generatePrimaryHeader(
             xlen=self.npix_x, 
             ylen=self.npix_y, 
@@ -333,23 +339,30 @@ class Chip(FocalPlane):
             pixel_size=self.pix_size,
             xcen=self.x_cen,
             ycen=self.y_cen,
-            extName='raw')
+            extName='SCI',
+            timestamp = timestamp,
+            exptime = exptime,
+            readoutTime = 40.)
         return h_prim, h_ext
 
-    def outputCal(self, img, ra_cen, dec_cen, img_rot, im_type, pointing_ID, date_obs, time_obs, output_dir, exptime=150.):
+    def outputCal(self, img, ra_cen, dec_cen, img_rot, im_type, pointing_ID, output_dir, exptime=150., timestamp = 1621915200):
         h_prim, h_ext = self.generateHeader(
             ra_cen=ra_cen,
             dec_cen=dec_cen,
             img_rot=img_rot,
             im_type=im_type,
             pointing_ID=pointing_ID,
-            date_obs=date_obs,
-            time_obs=time_obs,
-            exptime=exptime)
+            exptime=exptime,
+            timestamp = timestamp)
         hdu1 = fits.PrimaryHDU(header=h_prim)
         hdu1.add_checksum()
+        hdu1.header.comments['CHECKSUM'] = 'HDU checksum'
+        hdu1.header.comments['DATASUM'] = 'data unit checksum'
         hdu2 = fits.ImageHDU(img.array, header=h_ext)
         hdu2.add_checksum()
+        hdu2.header.comments['XTENSION'] = 'extension type'
+        hdu2.header.comments['CHECKSUM'] = 'HDU checksum'
+        hdu2.header.comments['DATASUM'] = 'data unit checksum'
         hdu1 = fits.HDUList([hdu1, hdu2])
         fname = os.path.join(output_dir, h_prim['FILENAME']+'.fits')
         hdu1.writeto(fname, output_verify='ignore', overwrite=True)
@@ -456,9 +469,9 @@ class Chip(FocalPlane):
             del cr_map
             # crmap_gsimg.write("%s/CosmicRay_%s_1.fits" % (chip_output.subdir, self.chipID))
             # crmap_gsimg.write("%s/CosmicRay_%s.fits" % (chip_output.subdir, self.chipID))
-            datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
-            date_obs = datetime_obs.strftime("%y%m%d")
-            time_obs = datetime_obs.strftime("%H%M%S")
+            # datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
+            # date_obs = datetime_obs.strftime("%y%m%d")
+            # time_obs = datetime_obs.strftime("%H%M%S")
             self.outputCal(
                 img=crmap_gsimg,
                 ra_cen=ra_cen,
@@ -466,10 +479,9 @@ class Chip(FocalPlane):
                 img_rot=img_rot,
                 im_type='CRS',
                 pointing_ID=pointing_ID,
-                date_obs=date_obs,
-                time_obs=time_obs,
                 output_dir=chip_output.subdir,
-                exptime=exptime)
+                exptime=exptime,
+                timestamp=timestamp_obs)
             del crmap_gsimg
 
         # Apply PRNU effect and output PRNU flat file:
@@ -503,21 +515,6 @@ class Chip(FocalPlane):
         if config["ins_effects"]["add_badcolumns"] == True:
             img = effects.BadColumns(img, seed=SeedBadColumns, chipid=self.chipID, logger=self.logger)
 
-        # Add Bias level
-        if config["ins_effects"]["add_bias"] == True:
-            if self.logger is not None:
-                self.logger.info("  Adding Bias level and 16-channel non-uniformity")
-            else:
-                print("  Adding Bias level and 16-channel non-uniformity")
-            if config["ins_effects"]["bias_16channel"] == True:
-                img = effects.AddBiasNonUniform16(img, 
-                    bias_level=float(self.bias_level), 
-                    nsecy = 2, nsecx=8, 
-                    seed=SeedBiasNonuni+self.chipID,
-                    logger=self.logger)
-            elif config["ins_effects"]["bias_16channel"] == False:
-                img += self.bias_level
-
         # Apply Nonlinearity on the chip image
         if config["ins_effects"]["non_linear"] == True:
             if self.logger is not None:
@@ -541,6 +538,21 @@ class Chip(FocalPlane):
             else:
                 print("  Apply CTE Effect")
             img = effects.CTE_Effect(GSImage=img, threshold=27)
+        
+        # Add Bias level
+        if config["ins_effects"]["add_bias"] == True:
+            if self.logger is not None:
+                self.logger.info("  Adding Bias level and 16-channel non-uniformity")
+            else:
+                print("  Adding Bias level and 16-channel non-uniformity")
+            if config["ins_effects"]["bias_16channel"] == True:
+                img = effects.AddBiasNonUniform16(img, 
+                    bias_level=float(self.bias_level), 
+                    nsecy = 2, nsecx=8, 
+                    seed=SeedBiasNonuni+self.chipID,
+                    logger=self.logger)
+            elif config["ins_effects"]["bias_16channel"] == False:
+                img += self.bias_level
 
         # Add Read-out Noise
         if config["ins_effects"]["add_readout"] == True:
@@ -555,7 +567,7 @@ class Chip(FocalPlane):
         else:
             print("  Applying Gain (and 16 channel non-uniformity) & Quantization", flush=True)
         if config["ins_effects"]["gain_16channel"] == True:
-            img = effects.ApplyGainNonUniform16(
+            img, self.gain_channel = effects.ApplyGainNonUniform16(
                 img, gain=self.gain, 
                 nsecy = 2, nsecx=8, 
                 seed=SeedGainNonuni+self.chipID,
@@ -610,7 +622,7 @@ class Chip(FocalPlane):
                 if config["ins_effects"]["add_badcolumns"] == True:
                     BiasCombImg = effects.BadColumns(BiasCombImg-float(self.bias_level)+5, seed=SeedBadColumns, chipid=self.chipID, logger=self.logger) + float(self.bias_level)-5
 
-                BiasCombImg = effects.ApplyGainNonUniform16(BiasCombImg, gain=self.gain, 
+                BiasCombImg, self.gain_channel = effects.ApplyGainNonUniform16(BiasCombImg, gain=self.gain,
                     nsecy = 2, nsecx=8, 
                     seed=SeedGainNonuni+self.chipID,
                     logger=self.logger)
@@ -622,9 +634,9 @@ class Chip(FocalPlane):
                 BiasCombImg.quantize()
                 BiasCombImg = galsim.ImageUS(BiasCombImg)
                 # BiasCombImg.write("%s/BiasImg_%s_%s_%s.fits" % (chip_output.subdir, BiasTag, self.chipID, i+1))
-                datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
-                date_obs = datetime_obs.strftime("%y%m%d")
-                time_obs = datetime_obs.strftime("%H%M%S")
+                # datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
+                # date_obs = datetime_obs.strftime("%y%m%d")
+                # time_obs = datetime_obs.strftime("%H%M%S")
                 timestamp_obs += 10 * 60
                 self.outputCal(
                     img=BiasCombImg,
@@ -633,10 +645,9 @@ class Chip(FocalPlane):
                     img_rot=img_rot,
                     im_type='BIAS',
                     pointing_ID=pointing_ID,
-                    date_obs=date_obs,
-                    time_obs=time_obs,
                     output_dir=chip_output.subdir,
-                    exptime=0.0)
+                    exptime=0.0,
+                    timestamp=timestamp_obs)
             del BiasCombImg
 
         # Export combined (ncombine, Vignetting + PRNU) & single vignetting flat-field file
@@ -717,7 +728,7 @@ class Chip(FocalPlane):
                     readout_noise = galsim.GaussianNoise(rng=rng_readout, sigma=self.read_noise)
                     FlatCombImg.addNoise(readout_noise)
 
-                FlatCombImg = effects.ApplyGainNonUniform16(FlatCombImg, gain=self.gain, 
+                FlatCombImg, self.gain_channel = effects.ApplyGainNonUniform16(FlatCombImg, gain=self.gain,
                     nsecy = 2, nsecx=8, 
                     seed=SeedGainNonuni+self.chipID,
                     logger=self.logger)
@@ -726,9 +737,9 @@ class Chip(FocalPlane):
                 FlatCombImg.quantize()
                 FlatCombImg = galsim.ImageUS(FlatCombImg)
                 # FlatCombImg.write("%s/FlatImg_%s_%s_%s.fits" % (chip_output.subdir, FlatTag, self.chipID, i+1))
-                datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
-                date_obs = datetime_obs.strftime("%y%m%d")
-                time_obs = datetime_obs.strftime("%H%M%S")
+                # datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
+                # date_obs = datetime_obs.strftime("%y%m%d")
+                # time_obs = datetime_obs.strftime("%H%M%S")
                 timestamp_obs += 10 * 60
                 self.outputCal(
                     img=FlatCombImg,
@@ -737,10 +748,9 @@ class Chip(FocalPlane):
                     img_rot=img_rot,
                     im_type='FLAT',
                     pointing_ID=pointing_ID,
-                    date_obs=date_obs,
-                    time_obs=time_obs,
                     output_dir=chip_output.subdir,
-                    exptime=self.flat_exptime)
+                    exptime=self.flat_exptime,
+                    timestamp=timestamp_obs)
 
             del FlatCombImg, FlatSingle, prnu_img
             # flat_img.replaceNegative(replace_value=0)
@@ -783,9 +793,9 @@ class Chip(FocalPlane):
                     cr_map[cr_map < 0] = 0
                     crmap_gsimg = galsim.Image(cr_map, dtype=np.uint16)
                     del cr_map
-                    datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
-                    date_obs = datetime_obs.strftime("%y%m%d")
-                    time_obs = datetime_obs.strftime("%H%M%S")
+                    # datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
+                    # date_obs = datetime_obs.strftime("%y%m%d")
+                    # time_obs = datetime_obs.strftime("%H%M%S")
                     self.outputCal(
                         img=crmap_gsimg,
                         ra_cen=ra_cen,
@@ -793,10 +803,9 @@ class Chip(FocalPlane):
                         img_rot=img_rot,
                         im_type='CRD',
                         pointing_ID=pointing_ID,
-                        date_obs=date_obs,
-                        time_obs=time_obs,
                         output_dir=chip_output.subdir,
-                        exptime=self.dark_exptime)
+                        exptime=self.dark_exptime,
+                        timestamp=timestamp_obs)
                     del crmap_gsimg
 
                 # Non-Linearity for Dark
@@ -839,7 +848,7 @@ class Chip(FocalPlane):
                     readout_noise = galsim.GaussianNoise(rng=rng_readout, sigma=self.read_noise)
                     DarkCombImg.addNoise(readout_noise)
 
-                DarkCombImg = effects.ApplyGainNonUniform16(
+                DarkCombImg, self.gain_channel = effects.ApplyGainNonUniform16(
                     DarkCombImg, gain=self.gain, 
                     nsecy = 2, nsecx=8, 
                     seed=SeedGainNonuni+self.chipID,
@@ -852,9 +861,9 @@ class Chip(FocalPlane):
                 DarkCombImg.quantize()
                 DarkCombImg = galsim.ImageUS(DarkCombImg)
                 # DarkCombImg.write("%s/DarkImg_%s_%s_%s.fits" % (chip_output.subdir, DarkTag, self.chipID, i+1))
-                datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
-                date_obs = datetime_obs.strftime("%y%m%d")
-                time_obs = datetime_obs.strftime("%H%M%S")
+                # datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
+                # date_obs = datetime_obs.strftime("%y%m%d")
+                # time_obs = datetime_obs.strftime("%H%M%S")
                 timestamp_obs += 10 * 60
                 self.outputCal(
                     img=DarkCombImg,
@@ -863,10 +872,9 @@ class Chip(FocalPlane):
                     img_rot=img_rot,
                     im_type='DARK',
                     pointing_ID=pointing_ID,
-                    date_obs=date_obs,
-                    time_obs=time_obs,
                     output_dir=chip_output.subdir,
-                    exptime=self.dark_exptime)
+                    exptime=self.dark_exptime,
+                    timestamp = timestamp_obs)
             del DarkCombImg
         # img = galsim.ImageUS(img)
 

ObservationSim/Instrument/Chip/Effects.py

@@ -152,6 +152,7 @@ def ApplyGainNonUniform16(GSImage, gain=1, nsecy = 2, nsecx=8, seed=202102, logg
     rg = Generator(PCG64(int(seed)))
     Random16 = (rg.random(nsecy*nsecx)-0.5)*0.04+1   # sigma~1%
     Gain16 = Random16.reshape((nsecy,nsecx))/gain
+    gain_array = np.ones(nsecy*nsecx)*gain
     if logger is not None:
         msg = str("Gain of 16 channels: " + str(Gain16))
         logger.info(msg)
@@ -163,7 +164,8 @@ def ApplyGainNonUniform16(GSImage, gain=1, nsecy = 2, nsecx=8, seed=202102, logg
     for rowi in range(nsecy):
         for coli in range(nsecx):
             GSImage.array[rowi*secsize_y:(rowi+1)*secsize_y,coli*secsize_x:(coli+1)*secsize_x] *= Gain16[rowi,coli]
-    return GSImage
+            gain_array[rowi*nsecx+coli] = 1/Gain16[rowi,coli]
+    return GSImage, gain_array
 
 
 def GainsNonUniform16(GSImage, gain=1, nsecy = 2, nsecx=8, seed=202102, logger=None):

ObservationSim/Instrument/Filter.py

@@ -2,9 +2,11 @@ import galsim
 import pylab as pl
 import os
 import numpy as np
+import gc
 
 from ObservationSim.Instrument._util import photonEnergy, calculateLimitMag
 from ObservationSim.Instrument.FilterParam import FilterParam
+from ObservationSim.Straylight import Straylight
 
 try:
     import importlib.resources as pkg_resources
@@ -38,6 +40,7 @@ class Filter(object):
         self.sky_background = filter_param.param[filter_type][5]
         self.mag_saturation = filter_param.param[filter_type][6]
         self.mag_limiting = filter_param.param[filter_type][7]
+        self.zodical_spec = None
 
     def is_too_bright(self, mag, margin=-2.5):
         return mag <= self.mag_saturation + margin
@@ -108,6 +111,24 @@ class Filter(object):
     def getSkyNoise(self, exptime, gain=1.):
         return self.sky_background * exptime / gain
 
+    def setFilterStrayLightPixel(self,jtime = 2460843., sat_pos = np.array([0,0,0]), pointing_radec = np.array([0,0]), sun_pos = np.array([0,0,0])):
+        sl = Straylight(jtime=jtime, sat_pos=sat_pos, pointing_radec=pointing_radec,sun_pos=sun_pos)
+        if self.filter_type in ["GU","GV","GI"]:
+            s_pix, spec = sl.calculateStrayLightGrating(grating = self.filter_type.upper())
+            if s_pix>0.8:
+                s_pix = 0.8
+            self.sky_background = s_pix
+            self.zodical_spec = spec
+        elif self.filter_type.lower() in ["nuv","u","g","r","i","z","y"]:
+            s_pix = sl.calculateStrayLightFilter(filter=self.filter_type.lower())
+            if s_pix>1:
+                s_pix = 1
+            self.sky_background = s_pix
+            self.zodical_spec = None
+
+        del sl
+        gc.collect()
+
     def update_limit_saturation_mags(self, exptime=150., psf_fwhm=0.1969, skyFn='sky_emiss_hubble_50_50_A.dat', chip=None):
         if self.filter_type in ["GI", "GV", "GU"]:
             return

ObservationSim/Instrument/FilterParam.py

@@ -26,7 +26,7 @@ class FilterParam(object):
         # 8) dim end magnitude
         if filter_param == None:
             filtP = {
-                    "nuv":   [2867.7,   705.4,  2470.0,   3270.0,  0.1404, 0.004, 15.7, 25.4],
+                    "NUV":   [2867.7,   705.4,  2470.0,   3270.0,  0.1404, 0.004, 15.7, 25.4],
                     "u":     [3601.1,   852.1,  3120.0,   4090.0,  0.2176, 0.021, 16.1, 25.4],
                     "g":     [4754.5,  1569.8,  3900.0,   5620.0,  0.4640, 0.164, 17.2, 26.3],
                     "r":     [6199.8,  1481.2,  5370.0,   7030.0,  0.5040, 0.207, 17.0, 26.0],

ObservationSim/Instrument/FocalPlane.py

@@ -86,15 +86,16 @@ class FocalPlane(object):
         if (xcen == None) or (ycen == None):
             xcen = self.cen_pix_x
             ycen = self.cen_pix_y
+        # dudx =  -np.cos(img_rot.rad) * pix_scale
+        # dudy =  -np.sin(img_rot.rad) * pix_scale
+        # dvdx =  -np.sin(img_rot.rad) * pix_scale
+        # dvdy =  +np.cos(img_rot.rad) * pix_scale
+
         dudx =  -np.cos(img_rot.rad) * pix_scale
-        dudy =  -np.sin(img_rot.rad) * pix_scale
+        dudy =  +np.sin(img_rot.rad) * pix_scale
         dvdx =  -np.sin(img_rot.rad) * pix_scale
-        dvdy =  +np.cos(img_rot.rad) * pix_scale
+        dvdy =  -np.cos(img_rot.rad) * pix_scale
         
-        # dudx =  +np.sin(img_rot.rad) * pix_scale
-        # dudy =  +np.cos(img_rot.rad) * pix_scale
-        # dvdx =  -np.cos(img_rot.rad) * pix_scale
-        # dvdy =  +np.sin(img_rot.rad) * pix_scale
         moscen = galsim.PositionD(x=xcen, y=ycen)
         sky_center = galsim.CelestialCoord(ra=ra*galsim.degrees, dec=dec*galsim.degrees)
         affine = galsim.AffineTransform(dudx, dudy, dvdx, dvdy, origin=moscen)

ObservationSim/Instrument/data/ccd/chip_definition.json

@@ -468,7 +468,7 @@
         "full_well": 90000
     },
     "15": {
-        "chip_name": "y-4",
+        "chip_name": "y-2",
         "pix_size": 0.01,
         "pix_scale": 0.074,
         "npix_x": 9216,
@@ -648,7 +648,7 @@
         "full_well": 90000
     },
     "25": {
-        "chip_name": "y-2",
+        "chip_name": "y-4",
         "pix_size": 0.01,
         "pix_scale": 0.074,
         "npix_x": 9216,