Skip to content
GitLab
Commit 78bf42f6 authored by Fang Yuedong's avatar Fang Yuedong
Browse files

Merge branch 'staylight_header' into 'develop' 

Straylight model

See merge request csst_sim/csst-simulation!14
parents 81589f9d 40ea0feb
Expand all Show whitespace changes
Inline Side-by-side
ObservationSim/Config/Header/ImageHeader.py
View file @ 78bf42f6
......@@ -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
......@@ -381,7 +385,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 +421,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 +437,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 +467,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 +534,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/global_header.header
View file @ 78bf42f6
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
View file @ 78bf42f6
......@@ -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
View file @ 78bf42f6
......@@ -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
......@@ -304,7 +306,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 +338,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 +468,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 +478,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:
......@@ -555,7 +566,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 +621,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 +633,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 +644,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 +727,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 +736,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 +747,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 +792,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 +802,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 +847,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 +860,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 +871,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
View file @ 78bf42f6
......@@ -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
View file @ 78bf42f6
......@@ -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
......@@ -108,6 +110,20 @@ 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())
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())
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/MockObject/SpecDisperser/SpecDisperser.py
View file @ 78bf42f6
......@@ -122,7 +122,7 @@ class SpecDisperser(object):
beam=beam)
### Account for pixel centering of the trace
yfrac_beam = ytrace_beam - floor(ytrace_beam)
yfrac_beam = ytrace_beam - floor(ytrace_beam+0.5)
ysens = lam_beam * 0
lam_index = argsort(lam_beam)
......@@ -155,7 +155,7 @@ class SpecDisperser(object):
sensitivity_beam = ysens
len_spec_x = len(dx)
len_spec_y = int(ceil(yfrac_beam[-1]) - floor(yfrac_beam[0]) + 1)
len_spec_y = int(ceil(ytrace_beam[-1]) - floor(ytrace_beam[0]) + 1)
beam_sh = (self.img_sh[0] + len_spec_y, self.img_sh[1] + len_spec_x)
modelf = zeros(product(beam_sh), dtype=float)
......@@ -165,9 +165,15 @@ class SpecDisperser(object):
dxpix = dx - dx[0] + x0[1]
dyc = cast[int](ytrace_beam)
dyc = cast[int](np.floor(ytrace_beam+0.5))
dypix = cast[int](np.floor(ytrace_beam - dyc[0] + x0[0] + 0.5))
frac_ids = yfrac_beam<0
dypix[frac_ids] = dypix[frac_ids] - 1
yfrac_beam[frac_ids] = 1+yfrac_beam[frac_ids]
dypix = dyc - dyc[0] + x0[0]
flat_index = idx[dypix, dxpix]
nonz = sensitivity_beam != 0
......
ObservationSim/MockObject/SpecDisperser/disperse_c/disperse.pyx
View file @ 78bf42f6
......@@ -86,12 +86,12 @@ def disperse_grism_object(np.ndarray[FTYPE_t, ndim=2] flam,
k1 = idxl[k]+j*shg[1]+i
if (k1 >= 0) & (k1 < nl):
flat_ids = k1*nlamb+k
full[k1] += fl_ij*yfrac[k]*flat_eff_all[flat_ids]
full[k1] += fl_ij*(1-yfrac[k])*flat_eff_all[flat_ids]
k2 = idxl[k]+(j-1)*shg[1]+i
k2 = idxl[k]+(j+1)*shg[1]+i
if (k2 >= 0) & (k2 < nl):
flat_ids = k2*nlamb+k
full[k2] += fl_ij*(1-yfrac[k])*flat_eff_all[flat_ids]
full[k2] += fl_ij*yfrac[k]*flat_eff_all[flat_ids]
else:
for i in range(0-x0[1], x0[1]):
......@@ -109,11 +109,11 @@ def disperse_grism_object(np.ndarray[FTYPE_t, ndim=2] flam,
for k in range(nk):
k1 = idxl[k]+j*shg[1]+i
if (k1 >= 0) & (k1 < nl):
full[k1] += ysens[k]*fl_ij*yfrac[k]
full[k1] += ysens[k]*fl_ij*(1-yfrac[k])
k2 = idxl[k]+(j-1)*shg[1]+i
k2 = idxl[k]+(j+1)*shg[1]+i
if (k2 >= 0) & (k2 < nl):
full[k2] += ysens[k]*fl_ij*(1-yfrac[k])
full[k2] += ysens[k]*fl_ij*yfrac[k]
return True
......
ObservationSim/MockObject/__init__.py
View file @ 78bf42f6
......@@ -4,5 +4,5 @@ from .CatalogBase import CatalogBase
from .Quasar import Quasar
from .Star import Star
from .Stamp import Stamp
from .SkybackgroundMap import *
# from .SkybackgroundMap import *
# from .CosmicRay import CosmicRay
ObservationSim/ObservationSim.py
View file @ 78bf42f6
......@@ -14,7 +14,7 @@ from ObservationSim.Config import config_dir, ChipOutput
from ObservationSim.Config.Header import generatePrimaryHeader, generateExtensionHeader
from ObservationSim.Instrument import Telescope, Filter, FilterParam, FocalPlane, Chip
from ObservationSim.Instrument.Chip import Effects
from ObservationSim.MockObject import calculateSkyMap_split_g
from ObservationSim.Straylight import calculateSkyMap_split_g
from ObservationSim.PSF import PSFGauss, FieldDistortion, PSFInterp
from ObservationSim._util import get_shear_field, makeSubDir_PointingList
from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position
......@@ -107,6 +107,12 @@ class Observation(object):
chip.img = galsim.ImageF(chip.npix_x, chip.npix_y)
chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin)
chip.img.wcs = wcs_fp
filt.setFilterStrayLightPixel(jtime = pointing.jdt, sat_pos = np.array([pointing.sat_x, pointing.sat_y, pointing.sat_z]), pointing_radec = np.array([pointing.ra,pointing.dec]), sun_pos = np.array([pointing.sun_x,pointing.sun_y,pointing.sun_z]))
print("========================sky pix========================\n")
print(filt.sky_background)
if chip.survey_type == "photometric":
sky_map = None
elif chip.survey_type == "spectroscopic":
......@@ -133,7 +139,8 @@ class Observation(object):
conf=chip.sls_conf,
pixelSize=chip.pix_scale,
isAlongY=0,
flat_cube=chip.flat_cube)
flat_cube=chip.flat_cube, zoldial_spec = filt.zodical_spec)
sky_map = sky_map+filt.sky_background
del flat_normal
if pointing.pointing_type == 'MS':
......@@ -176,7 +183,10 @@ class Observation(object):
pixel_size=chip.pix_size,
xcen=chip.x_cen,
ycen=chip.y_cen,
extName='raw')
extName='SCI',
timestamp = pointing.timestamp,
exptime = pointing.exp_time,
readoutTime = 40.)
for j in range(self.nobj):
......@@ -358,12 +368,36 @@ class Observation(object):
sat_pos=[pointing.sat_x, pointing.sat_y, pointing.sat_z],
sat_vel=[pointing.sat_vx, pointing.sat_vy, pointing.sat_vz],
chip_name=str(chip.chipID).rjust(2, '0'))
h_ext = generateExtensionHeader(
chip=chip,
xlen=chip.npix_x,
ylen=chip.npix_y,
ra=pointing.ra,
dec=pointing.dec,
pa=pointing.img_pa.deg,
gain=chip.gain,
readout=chip.read_noise,
dark=chip.dark_noise,
saturation=90000,
pixel_scale=chip.pix_scale,
pixel_size=chip.pix_size,
xcen=chip.x_cen,
ycen=chip.y_cen,
extName='SCI',
timestamp=pointing.timestamp,
exptime=pointing.exp_time,
readoutTime=40.)
chip.img = galsim.Image(chip.img.array, dtype=np.uint16)
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(chip.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(chip_output.subdir, h_prim['FILENAME'] + '.fits')
hdu1.writeto(fname, output_verify='ignore', overwrite=True)
......
ObservationSim/MockObject/SkybackgroundMap.py → ObservationSim/Straylight/SkybackgroundMap.py
View file @ 78bf42f6
......@@ -7,6 +7,7 @@ from astropy.table import Table
from scipy import interpolate
import galsim
import astropy.constants as cons
import os
......@@ -22,7 +23,7 @@ except ImportError:
###calculate sky map by sky SED
def calculateSkyMap_split_g(skyMap=None, blueLimit=4200, redLimit=6500, skyfn='sky_emiss_hubble_50_50_A.dat', conf=[''], pixelSize=0.074, isAlongY=0,
split_pos=3685, flat_cube = None):
split_pos=3685, flat_cube = None, zoldial_spec = None):
# skyMap = np.ones([yLen, xLen], dtype='float32')
#
# if isAlongY == 1:
......@@ -44,14 +45,24 @@ def calculateSkyMap_split_g(skyMap=None, blueLimit=4200, redLimit=6500, skyfn='s
fImg = galsim.Image(fimg)
try:
with pkg_resources.files('ObservationSim.MockObject.data').joinpath(skyfn) as data_path:
with pkg_resources.files('ObservationSim.Straylight.data.sky').joinpath(skyfn) as data_path:
skySpec = np.loadtxt(data_path)
except AttributeError:
with pkg_resources.path('ObservationSim.MockObject.data', skyfn) as data_path:
with pkg_resources.path('ObservationSim.Straylight.data.sky', skyfn) as data_path:
skySpec = np.loadtxt(data_path)
# skySpec = np.loadtxt(skyfn)
spec = Table(np.array([skySpec[:, 0], skySpec[:, 1]]).T, names=('WAVELENGTH', 'FLUX'))
if zoldial_spec is not None:
deltL = 0.5
lamb = np.arange(2000, 11000, deltL)
speci = interpolate.interp1d(zoldial_spec['WAVELENGTH'], zoldial_spec['FLUX'])
y = speci(lamb)
# erg/s/cm2/A --> photo/s/m2/A
s_flux = y * lamb / (cons.h.value * cons.c.value) * 1e-13
spec = Table(np.array([lamb, s_flux]).T, names=('WAVELENGTH', 'FLUX'))
if isAlongY == 0:
directParm = 0
if isAlongY ==1:
......@@ -266,10 +277,10 @@ def calculateSkyMap(xLen=9232, yLen=9126, blueLimit=4200, redLimit=6500,
fimg = np.zeros_like(skyMap)
fImg = galsim.Image(fimg)
try:
with pkg_resources.files('ObservationSim.MockObject.data').joinpath(skyfn) as data_path:
with pkg_resources.files('ObservationSim.Straylight.data.sky').joinpath(skyfn) as data_path:
skySpec = np.loadtxt(data_path)
except AttributeError:
with pkg_resources.path('ObservationSim.MockObject.data', skyfn) as data_path:
with pkg_resources.path('ObservationSim.Straylight.data.sky', skyfn) as data_path:
skySpec = np.loadtxt(data_path)
# skySpec = np.loadtxt(skyfn)
......
ObservationSim/Straylight/Straylight.py 0 → 100644
View file @ 78bf42f6
This diff is collapsed.
ObservationSim/Straylight/__init__.py 0 → 100755
View file @ 78bf42f6
from .Straylight import Straylight
from .SkybackgroundMap import *
ObservationSim/Straylight/data/Zodiacal_map1.dat 0 → 100644
View file @ 78bf42f6
0 0 5 10 15 20 25 30 45 60 75 90
0 30000 150000 7000 3140 1610 985 640 275 150 100 76
5 20000 12000 6000 2940 1540 945 625 271 150 100 76
10 16000 8000 4740 2470 1370 865 590 264 148 100 76
15 11500 6780 3440 1860 1110 755 525 251 146 100 76
20 6400 4480 2410 1410 910 635 454 237 141 99 76
25 3840 2830 1730 1100 749 545 410 223 136 97 76
30 2480 1870 1220 845 615 467 365 207 131 95 76
35 1650 1270 910 680 510 397 320 193 125 93 76
40 1180 940 700 530 416 338 282 179 120 92 76
45 910 730 555 442 356 292 250 166 116 90 76
60 505 442 352 292 243 209 183 134 104 86 76
75 338 317 269 227 196 172 151 116 93 82 76
90 259 251 225 193 166 147 132 104 86 79 76
105 212 210 197 170 150 133 119 96 82 77 76
120 188 186 177 154 138 125 113 90 77 74 76
135 179 178 166 147 134 122 110 90 77 73 76
150 179 178 165 148 137 127 116 96 79 72 76
165 196 192 179 165 151 141 131 104 82 72 76
180 230 212 195 178 163 148 134 105 83 72 76
\ No newline at end of file
ObservationSim/Straylight/data/sky/earthShine.dat 0 → 100644
View file @ 78bf42f6
#A erg/s/A/arcsec^2/cm^2
1000 2.41e-23
1100 4.38e-22
1200 4.01e-23
1300 7.41e-25
1400 4.29e-25
1500 4.16e-25
1600 2.55e-25
1700 7.89e-25
1800 9.33e-23
1900 4.39e-22
2000 1.01e-21
2100 1.6e-21
2200 7.49e-22
2300 3.32e-22
2400 2.5e-22
2500 2.39e-22
2600 5.62e-22
2700 6.77e-21
2800 2.03e-21
2900 4.32e-20
3000 9.34e-20
3100 2.07e-19
3200 3.6e-19
3300 4.27e-19
3400 6.4e-19
3500 8.2e-19
3600 1.06e-18
3700 1.22e-18
3800 1.23e-18
3900 1.52e-18
4000 2.38e-18
4250 2.38e-18
4500 2.86e-18
4750 2.79e-18
5000 2.63e-18
5250 2.67e-18
5500 2.58e-18
5750 2.54e-18
6000 2.42e-18
6250 2.26e-18
6500 2.17e-18
6750 2.07e-18
7000 1.93e-18
7250 1.85e-18
7500 1.74e-18
7750 1.63e-18
8000 1.56e-18
8250 1.48e-18
8500 1.35e-18
8750 1.31e-18
9000 1.22e-18
9250 1.15e-18
9500 1.1e-18
9750 1.04e-18
10000 1e-18
10250 9.45e-19
10500 9.04e-19
10750 8.41e-19
11000 8.03e-19
ObservationSim/Straylight/data/sky/earthshine1.dat 0 → 100644
View file @ 78bf42f6
2000 7.69E-22
2500 1.53E-21
3000 1.43E-19
3500 8.33E-19
4000 1.66E-18
4500 2.59E-18
5000 2.63E-18
5500 2.55E-18
6000 2.42E-18
7000 1.95E-18
8000 1.56E-18
9000 1.23E-18
10000 9.97E-19
11000 8.02E-19
12000 6.65E-19
13000 5.58E-19
14000 4.70E-19
15000 3.97E-19
16000 3.35E-19
17000 2.79E-19
ObservationSim/Straylight/data/sky/sky.dat 0 → 100644
View file @ 78bf42f6
2000 7.69E-22 7.94E-20 8.02E-20
2500 1.53E-21 3.83E-19 3.84E-19
3000 1.43E-19 1.63E-18 1.77E-18
3500 8.33E-19 2.72E-18 3.55E-18
4000 1.66E-18 3.12E-18 4.78E-18
4500 2.59E-18 4.97E-18 7.57E-18
5000 2.63E-18 5.07E-18 7.70E-18
5500 2.55E-18 5.17E-18 7.72E-18
6000 2.42E-18 5.14E-18 7.56E-18
7000 1.95E-18 4.48E-18 6.42E-18
8000 1.56E-18 3.82E-18 5.38E-18
9000 1.23E-18 3.18E-18 4.40E-18
10000 9.97E-19 2.70E-18 3.70E-18
11000 8.02E-19 2.26E-18 3.06E-18
12000 6.65E-19 1.94E-18 2.61E-18
13000 5.58E-19 1.68E-18 2.24E-18
14000 4.70E-19 1.46E-18 1.93E-18
15000 3.97E-19 1.26E-18 1.66E-18
16000 3.35E-19 1.09E-18 1.43E-18
17000 2.79E-19 9.27E-19 1.21E-18
\ No newline at end of file
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment