Merge branch 'master' of...

Merge branch 'master' of ssh://119.78.210.97:31415/data/simudata/CSSOSDataProductsSims/codes/CSST into customizable_catalog

Catalog/NGPCatalog.py

@@ -133,8 +133,6 @@ class NGPCatalog(CatalogBase):
             param['dec'] = dec_arr[igals]
             param['ra_orig'] = gals['ra_true'][igals]
             param['dec_orig'] = gals['dec_true'][igals]
-            if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200):
-                continue
             param['mag_use_normal'] = gals['mag_true_g_lsst'][igals]
             if param['mag_use_normal'] >= 26.5:
                 continue
@@ -173,6 +171,10 @@ class NGPCatalog(CatalogBase):
                 param['av'] = 0.6 * param['av'] / 3.0 # for quasar, av=[0, 0.2], 3.0=av.max-av.im
                 param['star'] = 2 # Quasar
 
+            # 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
+
             self.ids += 1
             # param['id'] = self.ids
             param['id'] = gals['galaxyID'][igals]

ObservationSim/Config/ChipOutput.py

@@ -23,12 +23,12 @@ class ChipOutput(object):
         
         # self.cat_name = 'MSC_' +  config["obs_setting"]["date_obs"] + config["obs_setting"]["time_obs"] + "_" + str(pointing_ID).rjust(7, '0') + "_" + self.chipLabel.rjust(2,'0') + ".cat"
 
-        self.cat_name = "MSC_1%s_chip_%s_filt_%s"%(str(pointing_ID).rjust(7, '0'), focal_plane.getChipLabel(chip.chipID), filt.filter_type) + ".cat"
+        self.cat_name = "MSC_1%s_chip_%s_filt_%s"%(str(pointing_ID).rjust(8, '0'), focal_plane.getChipLabel(chip.chipID), filt.filter_type) + ".cat"
 
         self.subdir = subdir
 
         # Setup logger for each chip
-        logger_filename = "MSC_1%s_chip_%s_filt_%s"%(str(pointing_ID).rjust(7, '0'), focal_plane.getChipLabel(chip.chipID), filt.filter_type) + ".log"
+        logger_filename = "MSC_1%s_chip_%s_filt_%s"%(str(pointing_ID).rjust(8, '0'), focal_plane.getChipLabel(chip.chipID), filt.filter_type) + ".log"
         self.logger = logging.getLogger()
         fh = logging.FileHandler(os.path.join(self.subdir, logger_filename), mode='w+', encoding='utf-8')
         fh.setLevel(logging.DEBUG)

ObservationSim/Instrument/Chip/Chip.py

@@ -33,6 +33,7 @@ class Chip(FocalPlane):
         self.dark_exptime = float(config["ins_effects"]['dark_exptime'])
         self.flat_exptime = float(config["ins_effects"]['flat_exptime'])
         self.readout_time = float(config["ins_effects"]['readout_time'])
+        self.full_well = int(config["ins_effects"]["full_well"])
 
         self.logger = logger
 
@@ -95,16 +96,17 @@ class Chip(FocalPlane):
         if filter_type in ['g', 'r', 'GV']: filename = 'Astro_MB.txt'
         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 == '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
+        # 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
         
         # path = os.path.join(self.ccdEffCurve_dir, filename)
         # table = Table.read(path, format='ascii')
         with pkg_resources.path('ObservationSim.Instrument.data.ccd', filename) as ccd_path:
             table = Table.read(ccd_path, format='ascii')
-        throughput = galsim.LookupTable(x=table['col1'], f=table['col2']*mirror_eff, interpolant='linear')
+        # throughput = galsim.LookupTable(x=table['col1'], f=table['col2']*mirror_eff, interpolant='linear')
+        throughput = galsim.LookupTable(x=table['col1'], f=table['col2'], interpolant='linear')
         bandpass = galsim.Bandpass(throughput, wave_type='nm')
         return bandpass
 
@@ -494,7 +496,7 @@ class Chip(FocalPlane):
         # Output images for calibration pointing
         ######################################################################################
         # Bias output
-        if config["output_setting"]["bias_output"] == True and pointing_type=='CAL':
+        if config["ins_effects"]["add_bias"] == True and config["output_setting"]["bias_output"] == True and pointing_type=='CAL':
             if self.logger is not None:
                 self.logger.info("  Output N frame Bias files")
             else:
@@ -562,7 +564,7 @@ class Chip(FocalPlane):
             del BiasCombImg
 
         # Export combined (ncombine, Vignetting + PRNU) & single vignetting flat-field file
-        if config["output_setting"]["flat_output"] == True and pointing_type=='CAL':
+        if config["ins_effects"]["flat_fielding"] == True and config["output_setting"]["flat_output"] == True and pointing_type=='CAL':
             if self.logger is not None:
                 self.logger.info("  Output N frame Flat-Field files")
             else:
@@ -671,7 +673,7 @@ class Chip(FocalPlane):
             del flat_img
 
         # Export Dark current images
-        if config["output_setting"]["dark_output"] == True and pointing_type=='CAL':
+        if config["ins_effects"]["add_dark"] == True and config["output_setting"]["dark_output"] == True and pointing_type=='CAL':
             if self.logger is not None:
                 self.logger.info("  Output N frame Dark Current files")
             else:

ObservationSim/Instrument/Filter.py

@@ -3,7 +3,7 @@ import pylab as pl
 import os
 import numpy as np
 
-from ObservationSim.Instrument._util import photonEnergy
+from ObservationSim.Instrument._util import photonEnergy, calculateLimitMag
 from ObservationSim.Instrument.FilterParam import FilterParam
 
 try:
@@ -13,7 +13,7 @@ except ImportError:
     import importlib_resources as pkg_resources
 
 class Filter(object):
-    def __init__(self, filter_id, filter_type, filter_param, ccd_bandpass):
+    def __init__(self, filter_id, filter_type, filter_param, ccd_bandpass=None):
         self.filter_id = filter_id
         self.filter_type = filter_type
         self.ccd_bandpass = ccd_bandpass
@@ -36,14 +36,15 @@ class Filter(object):
         self.efficiency = filter_param.param[filter_type][4]
         self.sky_background = filter_param.param[filter_type][5]
         self.mag_saturation = filter_param.param[filter_type][6]
-        self.mag_dim = filter_param.param[filter_type][7]
+        self.mag_limiting = filter_param.param[filter_type][7]
         # self.filter_dir = filter_param.filter_dir
 
     def is_too_bright(self, mag):
-        return mag <= self.mag_saturation - 1.0
+        # return mag <= self.mag_saturation - 1.0
+        return mag <= 14.0
 
     def is_too_dim(self, mag):
-        return mag >= self.mag_dim + 1.0
+        return mag >= self.mag_limiting + 1.0
 
     def _get_bandpasses(self, filter_dir=None, unit='A'):
         if self.filter_id < 7:  # Photometric
@@ -51,8 +52,10 @@ class Filter(object):
             # filter_file = os.path.join(filter_dir, self.filter_type+".dat")
             # bandpass_full = galsim.Bandpass(filter_file, wave_type=unit)
             with pkg_resources.path('ObservationSim.Instrument.data.filters', self.filter_type.lower() + '.txt') as filter_file:
+                self.filter_bandpass = galsim.Bandpass(str(filter_file), wave_type=unit)
+            with pkg_resources.path('ObservationSim.Instrument.data.throughputs', self.filter_type.lower() + '_throughput.txt') as filter_file:
                 bandpass_full = galsim.Bandpass(str(filter_file), wave_type=unit)
-            bandpass_full = bandpass_full * self.ccd_bandpass
+            # bandpass_full = bandpass_full * self.ccd_bandpass
 
             # Get sub-bandpasses
             bandpass_sub_list = []
@@ -95,3 +98,18 @@ class Filter(object):
 
     def getSkyNoise(self, exptime, gain=1.):
         return self.sky_background * exptime / gain
+
+    def update_limit_saturation_mags(self, exptime=150., psf_fwhm=0.1969, skyFn='sky_emiss_hubble_50_50_A.dat', chip=None):
+        if chip is not None:
+            pix_scale = chip.pix_scale
+            read_noise = chip.read_noise
+            dark_noise = chip.dark_noise
+            full_well = chip.full_well
+        else:
+            pix_scale = 0.074
+            read_noise = 5.0
+            dark_noise = 0.02
+            full_well = 90000
+        throughput_file = self.filter_type.lower() + '_throughput.txt'
+
+        self.mag_limiting, self.mag_saturation = calculateLimitMag(psf_fwhm=psf_fwhm, pixelSize=pix_scale, throughputFn=throughput_file, readout=5.0, skyFn=skyFn, darknoise=dark_noise, exTime=exptime, fw=full_well)

ObservationSim/Instrument/_util.py

 import numpy as np
 import os
+import math
+from pylab import *
+from scipy import interpolate
+
+try:
+    import importlib.resources as pkg_resources
+except ImportError:
+    # Try backported to PY<37 'importlib_resources'
+    import importlib_resources as pkg_resources
 
 VC_A = 2.99792458e+18  # speed of light: A/s
 VC_M = 2.99792458e+8   # speed of light: m/s
@@ -16,3 +25,84 @@ def photonEnergy(lambd):
 	nu = VC_A / lambd
 	eph = H_PLANK * nu
 	return eph
+
+'''
+description: 
+param {*} aperture: unit m, default 2 m
+param {*} psf_fwhm: psf fwhm, default 0.1969"
+param {*} pixelSize: pixel size, default 0.074"
+param {*} pmRation: the ratio of souce flux in the limit mag calculation
+param {*} throughputFn: throuput file name
+param {*} readout: unit, e-/pixel
+param {*} skyFn: sky sed file name, average of hst, 'sky_emiss_hubble_50_50_A.dat'
+param {*} darknoise: unit, e-/pixel/s
+param {*} exTime: exposure time one time, default 150s
+param {*} exNum: exposure number, defautl 1
+param {*} fw, full well value( or saturation value),default 90000e-/pixel
+return {*} limit mag and saturation mag
+'''
+
+def calculateLimitMag(aperture = 2.0, psf_fwhm = 0.1969,pixelSize = 0.074, pmRation = 0.8, throughputFn = 'i_throughput.txt', readout = 5.0, skyFn= 'sky_emiss_hubble_50_50_A.dat', darknoise = 0.02,exTime = 150, exNum = 1, fw = 90000):
+	with pkg_resources.path('ObservationSim.Instrument.data.throughputs', throughputFn) as data_file:
+		throughput_f = np.loadtxt(data_file)
+	thr_i = interpolate.interp1d(throughput_f[:,0]/10, throughput_f[:,1]); # wavelength in anstrom
+	f_s = 200
+	f_e = 1100
+	delt_f = 0.5
+
+	data_num = int((f_e-f_s)/delt_f+1)
+
+	eff = np.zeros([data_num,2])
+	eff[:,0] = np.arange(f_s,f_e+delt_f,delt_f)
+	eff[:,1] = thr_i(eff[:,0])
+
+	wave = np.arange(f_s,f_e+delt_f,delt_f)
+	wavey = np.ones(wave.shape[0])
+
+	with pkg_resources.path('ObservationSim.Instrument.data.throughputs', skyFn) as data_file:
+		skydata = np.loadtxt(data_file)
+	skydatai = interpolate.interp1d(skydata[:,0]/10, skydata[:,1]*10)
+
+	sky_data = np.zeros([data_num,2])
+	sky_data[:,0] = np.arange(f_s,f_e+delt_f,delt_f)
+	sky_data[:,1] = skydatai(sky_data[:,0])
+
+	flux_sky = trapz((sky_data[:,1])*eff[:,1],sky_data[:,0])
+	skyPix = flux_sky*pixelSize*pixelSize*pi*(aperture*aperture/4)
+
+	
+	###limit mag
+
+	r_pix = psf_fwhm*0.7618080243778568/pixelSize # radius RE80, pixel
+	cnum = math.pi * r_pix * r_pix
+	sn = 5
+
+	d = skyPix*exTime*exNum*cnum + darknoise*exTime*exNum*cnum+readout*readout*cnum*exNum
+	a=1
+	b=-sn*sn
+	c=-sn*sn*d
+
+	flux = (-b+sqrt(b*b-4*a*c))/(2*a)/pmRation
+	limitMag = -2.5*log10(flux/(54799275581.04437 * trapz(wavey*eff[:,1]/wave,wave, 0.1)*exTime*exNum*pi*(aperture/2)*(aperture/2)))
+
+
+	### saturation mag
+
+	from astropy.modeling.models import Gaussian2D
+	m_size = int(20 * psf_fwhm/pixelSize)
+	if m_size%2 == 0:
+		m_size + 1
+	m_cen = m_size//2
+	psf_sigma = psf_fwhm/2.355/pixelSize
+
+	gaussShape = Gaussian2D(1, m_cen, m_cen, psf_sigma, psf_sigma)
+	yp, xp = np.mgrid[0:m_size, 0:m_size]
+	psfMap = gaussShape(xp, yp)
+	maxRatio = np.amax(psfMap)/np.sum(psfMap)
+	print(maxRatio)
+
+	flux_sat = fw/maxRatio*exNum
+	satMag = -2.5*log10(flux_sat/(54799275581.04437 * trapz(wavey*eff[:,1]/wave,wave, 0.1)*exTime*exNum*pi*(aperture/2)*(aperture/2)));
+
+
+	return limitMag , satMag
\ No newline at end of file

ObservationSim/ObservationSim.py

@@ -47,8 +47,7 @@ class Observation(object):
             filter_id, filter_type = chip.getChipFilter()
             filt = Filter(filter_id=filter_id, 
                 filter_type=filter_type, 
-                filter_param=self.filter_param, 
-                ccd_bandpass=chip.effCurve)
+                filter_param=self.filter_param)
             self.chip_list.append(chip)
             self.filter_list.append(filt)
 
@@ -76,6 +75,9 @@ class Observation(object):
         chip_output.logger.info('Chip : %d' % chip.chipID)
         chip_output.logger.info(':::::::::::::::::::::::::::END:::::::::::::::::::::::::::::::::::')
 
+        # Update the limiting magnitude using exposure time in pointing
+        filt.update_limit_saturation_mags(exptime=pointing.exp_time, chip=chip)
+
         if self.config["psf_setting"]["psf_model"] == "Gauss":
             psf_model = PSFGauss(chip=chip)
         elif self.config["psf_setting"]["psf_model"] == "Interp":
@@ -195,7 +197,8 @@ class Observation(object):
                     continue
 
                 # Exclude very bright/dim objects (for now)
-                if filt.is_too_bright(mag=obj.getMagFilter(filt)):
+                # if filt.is_too_bright(mag=obj.getMagFilter(filt)):
+                if filt.is_too_bright(mag=obj.mag_use_normal):
                     # print("obj too birght!!", flush=True)
                     if obj.type != 'galaxy':
                         bright_obj += 1

setup.py

@@ -47,6 +47,7 @@ setup(name='CSSTSim',
         'ObservationSim.Instrument.data': ['*.txt', '*.dat'], 
         'ObservationSim.Instrument.data.ccd': ['*.txt'], 
         'ObservationSim.Instrument.data.filters': ['*.txt', '*.list', '*.dat'], 
+        'ObservationSim.Instrument.data.throughputs': ['*.txt', '*.dat'], 
         'ObservationSim.Instrument.data.sls_conf': ['*.conf', '*.fits'],
         'Catalog.data': ['*.fits'],
     },