fix discrepancy by applying saturation cut in each band, add limiting...

fix discrepancy by applying saturation cut in each band, add limiting magnitude calculator, fix randomization problem in SED assignment

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

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
@@ -15,4 +24,85 @@ def photonEnergy(lambd):
 	"""
 	nu = VC_A / lambd
 	eph = H_PLANK * nu
-	return eph
\ No newline at end of file
+	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":
@@ -196,7 +198,7 @@ class Observation(object):
 
                 # Exclude very bright/dim objects (for now)
                 # if filt.is_too_bright(mag=obj.getMagFilter(filt)):
-                if obj.getMagFilter(filt) <= 14.0:
+                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'],
     },