Skip to content
GitLab
Commit 00a96bcb authored by Zhang Xin's avatar Zhang Xin
Browse files

add model: straylight; fix header file consistent with doc:KSC-00-JK-0001-02.01.pdf

parent 81589f9d
Show whitespace changes
Inline Side-by-side
ObservationSim/Config/Header/ImageHeader.py
View file @ 00a96bcb
...@@ -18,6 +18,10 @@ from astropy.wcs.utils import fit_wcs_from_points ...@@ -18,6 +18,10 @@ from astropy.wcs.utils import fit_wcs_from_points
from astropy.time import Time from astropy.time import Time
from astropy import wcs from astropy import wcs
from datetime import datetime
# import socket
import platform
def chara2digit(char): def chara2digit(char):
""" Function to judge and convert characters to digitals """ Function to judge and convert characters to digitals
...@@ -381,7 +385,7 @@ def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec ...@@ -381,7 +385,7 @@ def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec
h_prim['EXPTIME'] = exptime h_prim['EXPTIME'] = exptime
# Define file types # 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':'cosmic_ray', 'CRD':'CRD','CALS':'CALS','CALF':'CALF'}
h_prim['FILETYPE'] = file_type[im_type] h_prim['FILETYPE'] = file_type[im_type]
co = coord.SkyCoord(ra, dec, unit='deg') co = coord.SkyCoord(ra, dec, unit='deg')
...@@ -417,7 +421,7 @@ def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec ...@@ -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] 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[ # h_prim['FILENAME'] = 'CSST_MSC_MS_' + im_type + '_' + file_start_time + '_' + file_end_time + '_' + OBS_id + '_' + CCDID[
# k - 1].rjust(2, '0') + '_L0_V01' # 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] h_prim['POSI0_X'] = sat_pos[0]
...@@ -433,12 +437,15 @@ def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec ...@@ -433,12 +437,15 @@ def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec
# Get version of CSSTSim Package # Get version of CSSTSim Package
from pkg_resources import get_distribution from pkg_resources import get_distribution
# h_prim['SIM_VER'] = (get_distribution("CSSTSim").version, "Version of CSST MSC simulation software") # 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 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, 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' 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') 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 ...@@ -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['CCDCHIP'] = CCDID[k - 1].rjust(2, '0')
# h_ext['CCDLABEL'] = filters[k-1] + '-' + filterID[k-1] # h_ext['CCDLABEL'] = filters[k-1] + '-' + filterID[k-1]
# h_ext['FILTER'] = filters[k-1] # h_ext['FILTER'] = filters[k-1]
h_ext['CCDCHIP'] = chip.chip_name h_ext['CHIPID'] = str(chip.chipID).rjust(2, '0')
h_ext['CCDLABEL'] = chip.filter_type + '-' + str(chip.chipID).rjust(2, '0') h_ext['CHIPLAB'] = chip.chip_name
h_ext['FILTER'] = chip.filter_type h_ext['FILTER'] = chip.filter_type
h_ext['NAXIS1'] = xlen h_ext['NAXIS1'] = xlen
h_ext['NAXIS2'] = ylen h_ext['NAXIS2'] = ylen
h_ext['EXTNAME'] = extName h_ext['EXTNAME'] = extName
h_ext['GAIN1'] = gain h_ext['GAIN01'] = chip.gain_channel[0]
h_ext['GAIN2'] = gain h_ext['GAIN02'] = chip.gain_channel[1]
h_ext['GAIN3'] = gain h_ext['GAIN03'] = chip.gain_channel[2]
h_ext['GAIN4'] = gain h_ext['GAIN04'] = chip.gain_channel[3]
h_ext['GAIN5'] = gain h_ext['GAIN05'] = chip.gain_channel[4]
h_ext['GAIN6'] = gain h_ext['GAIN06'] = chip.gain_channel[5]
h_ext['GAIN7'] = gain h_ext['GAIN07'] = chip.gain_channel[6]
h_ext['GAIN8'] = gain h_ext['GAIN08'] = chip.gain_channel[7]
h_ext['GAIN9'] = gain h_ext['GAIN09'] = chip.gain_channel[8]
h_ext['GAIN10'] = gain h_ext['GAIN10'] = chip.gain_channel[9]
h_ext['GAIN11'] = gain h_ext['GAIN11'] = chip.gain_channel[10]
h_ext['GAIN12'] = gain h_ext['GAIN12'] = chip.gain_channel[11]
h_ext['GAIN13'] = gain h_ext['GAIN13'] = chip.gain_channel[12]
h_ext['GAIN14'] = gain h_ext['GAIN14'] = chip.gain_channel[13]
h_ext['GAIN15'] = gain h_ext['GAIN15'] = chip.gain_channel[14]
h_ext['GAIN16'] = gain h_ext['GAIN16'] = chip.gain_channel[15]
h_ext['RDNOIS1'] = readout h_ext['RON01'] = readout
h_ext['RDNOIS2'] = readout h_ext['RON02'] = readout
h_ext['RDNOIS3'] = readout h_ext['RON03'] = readout
h_ext['RDNOIS4'] = readout h_ext['RON04'] = readout
h_ext['RDNOIS5'] = readout h_ext['RON05'] = readout
h_ext['RDNOIS6'] = readout h_ext['RON06'] = readout
h_ext['RDNOIS7'] = readout h_ext['RON07'] = readout
h_ext['RDNOIS8'] = readout h_ext['RON08'] = readout
h_ext['RDNOIS9'] = readout h_ext['RON09'] = readout
h_ext['RDNOIS10'] = readout h_ext['RON10'] = readout
h_ext['RDNOIS11'] = readout h_ext['RON11'] = readout
h_ext['RDNOIS12'] = readout h_ext['RON12'] = readout
h_ext['RDNOIS13'] = readout h_ext['RON13'] = readout
h_ext['RDNOIS14'] = readout h_ext['RON14'] = readout
h_ext['RDNOIS15'] = readout h_ext['RON15'] = readout
h_ext['RDNOIS16'] = readout h_ext['RON16'] = readout
h_ext['PIXSCAL1'] = pixel_scale h_ext['PIXSCAL1'] = pixel_scale
h_ext['PIXSCAL2'] = 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 # 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, 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 ...@@ -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['CD1_2'] = header_wcs['CD1_2']
h_ext['CD2_1'] = header_wcs['CD2_1'] h_ext['CD2_1'] = header_wcs['CD2_1']
h_ext['CD2_2'] = header_wcs['CD2_2'] h_ext['CD2_2'] = header_wcs['CD2_2']
h_ext['EQUINOX'] = header_wcs['EQUINOX'] # h_ext['EQUINOX'] = header_wcs['EQUINOX']
h_ext['WCSDIM'] = header_wcs['WCSDIM'] # h_ext['WCSDIM'] = header_wcs['WCSDIM']
h_ext['CTYPE1'] = header_wcs['CTYPE1'] h_ext['CTYPE1'] = header_wcs['CTYPE1']
h_ext['CTYPE2'] = header_wcs['CTYPE2'] h_ext['CTYPE2'] = header_wcs['CTYPE2']
......
ObservationSim/Config/Header/extension_header.header
View file @ 00a96bcb
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
View file @ 00a96bcb
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 @ 00a96bcb
...@@ -3,16 +3,18 @@ import galsim ...@@ -3,16 +3,18 @@ import galsim
from astropy.time import Time from astropy.time import Time
class Pointing(object): 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.id = id
self.ra = ra self.ra = ra
self.dec = dec self.dec = dec
self.img_pa = img_pa * galsim.degrees self.img_pa = img_pa * galsim.degrees
self.timestamp = timestamp self.timestamp = timestamp
self.sat_x, self.sat_y, self.sat_z = sat_x, sat_y, sat_z 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.sat_vx, self.sat_vy, self.sat_vz = sat_vx, sat_vy, sat_vz
self.exp_time = exp_time self.exp_time = exp_time
self.pointing_type = pointing_type self.pointing_type = pointing_type
self.jdt = 0.
def read_pointing_columns(self, columns, id=0, t=1621915200, pointing_type='MS'): def read_pointing_columns(self, columns, id=0, t=1621915200, pointing_type='MS'):
self.id = id self.id = id
...@@ -24,10 +26,14 @@ class Pointing(object): ...@@ -24,10 +26,14 @@ class Pointing(object):
if col_len > 5: if col_len > 5:
jdt = np.double(columns[5]) jdt = np.double(columns[5])
t_temp = Time(jdt, format='jd') t_temp = Time(jdt, format='jd')
self.jdt = jdt
self.timestamp = t_temp.unix self.timestamp = t_temp.unix
self.sat_x = float(columns[6]) self.sat_x = float(columns[6])
self.sat_y = float(columns[7]) self.sat_y = float(columns[7])
self.sat_z = float(columns[8]) 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_vx = float(columns[15])
self.sat_vy = float(columns[16]) self.sat_vy = float(columns[16])
self.sat_vz = float(columns[17]) self.sat_vz = float(columns[17])
......
ObservationSim/Instrument/Chip/Chip.py
View file @ 00a96bcb
...@@ -27,7 +27,9 @@ class Chip(FocalPlane): ...@@ -27,7 +27,9 @@ class Chip(FocalPlane):
# self.npix_x = 9216 # self.npix_x = 9216
# self.npix_y = 9232 # self.npix_y = 9232
# self.pix_scale = 0.074 # pixel scale # 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._set_attributes_from_config(config)
self.logger = logger self.logger = logger
...@@ -304,7 +306,10 @@ class Chip(FocalPlane): ...@@ -304,7 +306,10 @@ class Chip(FocalPlane):
if self.chipID == 30: confFile = ['CSST_GI8.conf', 'CSST_GI7.conf'] if self.chipID == 30: confFile = ['CSST_GI8.conf', 'CSST_GI7.conf']
return confFile 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( h_prim = generatePrimaryHeader(
xlen=self.npix_x, xlen=self.npix_x,
ylen=self.npix_y, ylen=self.npix_y,
...@@ -333,23 +338,30 @@ class Chip(FocalPlane): ...@@ -333,23 +338,30 @@ class Chip(FocalPlane):
pixel_size=self.pix_size, pixel_size=self.pix_size,
xcen=self.x_cen, xcen=self.x_cen,
ycen=self.y_cen, ycen=self.y_cen,
extName='raw') extName='SCI',
timestamp = timestamp,
exptime = exptime,
readoutTime = 40.)
return h_prim, h_ext 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( h_prim, h_ext = self.generateHeader(
ra_cen=ra_cen, ra_cen=ra_cen,
dec_cen=dec_cen, dec_cen=dec_cen,
img_rot=img_rot, img_rot=img_rot,
im_type=im_type, im_type=im_type,
pointing_ID=pointing_ID, pointing_ID=pointing_ID,
date_obs=date_obs, exptime=exptime,
time_obs=time_obs, timestamp = timestamp)
exptime=exptime)
hdu1 = fits.PrimaryHDU(header=h_prim) hdu1 = fits.PrimaryHDU(header=h_prim)
hdu1.add_checksum() 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 = fits.ImageHDU(img.array, header=h_ext)
hdu2.add_checksum() 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]) hdu1 = fits.HDUList([hdu1, hdu2])
fname = os.path.join(output_dir, h_prim['FILENAME']+'.fits') fname = os.path.join(output_dir, h_prim['FILENAME']+'.fits')
hdu1.writeto(fname, output_verify='ignore', overwrite=True) hdu1.writeto(fname, output_verify='ignore', overwrite=True)
...@@ -456,9 +468,9 @@ class Chip(FocalPlane): ...@@ -456,9 +468,9 @@ class Chip(FocalPlane):
del cr_map del cr_map
# crmap_gsimg.write("%s/CosmicRay_%s_1.fits" % (chip_output.subdir, self.chipID)) # 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)) # crmap_gsimg.write("%s/CosmicRay_%s.fits" % (chip_output.subdir, self.chipID))
datetime_obs = datetime.utcfromtimestamp(timestamp_obs) # datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
date_obs = datetime_obs.strftime("%y%m%d") # date_obs = datetime_obs.strftime("%y%m%d")
time_obs = datetime_obs.strftime("%H%M%S") # time_obs = datetime_obs.strftime("%H%M%S")
self.outputCal( self.outputCal(
img=crmap_gsimg, img=crmap_gsimg,
ra_cen=ra_cen, ra_cen=ra_cen,
...@@ -466,10 +478,9 @@ class Chip(FocalPlane): ...@@ -466,10 +478,9 @@ class Chip(FocalPlane):
img_rot=img_rot, img_rot=img_rot,
im_type='CRS', im_type='CRS',
pointing_ID=pointing_ID, pointing_ID=pointing_ID,
date_obs=date_obs,
time_obs=time_obs,
output_dir=chip_output.subdir, output_dir=chip_output.subdir,
exptime=exptime) exptime=exptime,
timestamp=timestamp_obs)
del crmap_gsimg del crmap_gsimg
# Apply PRNU effect and output PRNU flat file: # Apply PRNU effect and output PRNU flat file:
...@@ -555,7 +566,7 @@ class Chip(FocalPlane): ...@@ -555,7 +566,7 @@ class Chip(FocalPlane):
else: else:
print(" Applying Gain (and 16 channel non-uniformity) & Quantization", flush=True) print(" Applying Gain (and 16 channel non-uniformity) & Quantization", flush=True)
if config["ins_effects"]["gain_16channel"] == True: if config["ins_effects"]["gain_16channel"] == True:
img = effects.ApplyGainNonUniform16( img, self.gain_channel = effects.ApplyGainNonUniform16(
img, gain=self.gain, img, gain=self.gain,
nsecy = 2, nsecx=8, nsecy = 2, nsecx=8,
seed=SeedGainNonuni+self.chipID, seed=SeedGainNonuni+self.chipID,
...@@ -610,7 +621,7 @@ class Chip(FocalPlane): ...@@ -610,7 +621,7 @@ class Chip(FocalPlane):
if config["ins_effects"]["add_badcolumns"] == True: 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.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, nsecy = 2, nsecx=8,
seed=SeedGainNonuni+self.chipID, seed=SeedGainNonuni+self.chipID,
logger=self.logger) logger=self.logger)
...@@ -622,9 +633,9 @@ class Chip(FocalPlane): ...@@ -622,9 +633,9 @@ class Chip(FocalPlane):
BiasCombImg.quantize() BiasCombImg.quantize()
BiasCombImg = galsim.ImageUS(BiasCombImg) BiasCombImg = galsim.ImageUS(BiasCombImg)
# BiasCombImg.write("%s/BiasImg_%s_%s_%s.fits" % (chip_output.subdir, BiasTag, self.chipID, i+1)) # BiasCombImg.write("%s/BiasImg_%s_%s_%s.fits" % (chip_output.subdir, BiasTag, self.chipID, i+1))
datetime_obs = datetime.utcfromtimestamp(timestamp_obs) # datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
date_obs = datetime_obs.strftime("%y%m%d") # date_obs = datetime_obs.strftime("%y%m%d")
time_obs = datetime_obs.strftime("%H%M%S") # time_obs = datetime_obs.strftime("%H%M%S")
timestamp_obs += 10 * 60 timestamp_obs += 10 * 60
self.outputCal( self.outputCal(
img=BiasCombImg, img=BiasCombImg,
...@@ -633,10 +644,9 @@ class Chip(FocalPlane): ...@@ -633,10 +644,9 @@ class Chip(FocalPlane):
img_rot=img_rot, img_rot=img_rot,
im_type='BIAS', im_type='BIAS',
pointing_ID=pointing_ID, pointing_ID=pointing_ID,
date_obs=date_obs,
time_obs=time_obs,
output_dir=chip_output.subdir, output_dir=chip_output.subdir,
exptime=0.0) exptime=0.0,
timestamp=timestamp_obs)
del BiasCombImg del BiasCombImg
# Export combined (ncombine, Vignetting + PRNU) & single vignetting flat-field file # Export combined (ncombine, Vignetting + PRNU) & single vignetting flat-field file
...@@ -717,7 +727,7 @@ class Chip(FocalPlane): ...@@ -717,7 +727,7 @@ class Chip(FocalPlane):
readout_noise = galsim.GaussianNoise(rng=rng_readout, sigma=self.read_noise) readout_noise = galsim.GaussianNoise(rng=rng_readout, sigma=self.read_noise)
FlatCombImg.addNoise(readout_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, nsecy = 2, nsecx=8,
seed=SeedGainNonuni+self.chipID, seed=SeedGainNonuni+self.chipID,
logger=self.logger) logger=self.logger)
...@@ -726,9 +736,9 @@ class Chip(FocalPlane): ...@@ -726,9 +736,9 @@ class Chip(FocalPlane):
FlatCombImg.quantize() FlatCombImg.quantize()
FlatCombImg = galsim.ImageUS(FlatCombImg) FlatCombImg = galsim.ImageUS(FlatCombImg)
# FlatCombImg.write("%s/FlatImg_%s_%s_%s.fits" % (chip_output.subdir, FlatTag, self.chipID, i+1)) # FlatCombImg.write("%s/FlatImg_%s_%s_%s.fits" % (chip_output.subdir, FlatTag, self.chipID, i+1))
datetime_obs = datetime.utcfromtimestamp(timestamp_obs) # datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
date_obs = datetime_obs.strftime("%y%m%d") # date_obs = datetime_obs.strftime("%y%m%d")
time_obs = datetime_obs.strftime("%H%M%S") # time_obs = datetime_obs.strftime("%H%M%S")
timestamp_obs += 10 * 60 timestamp_obs += 10 * 60
self.outputCal( self.outputCal(
img=FlatCombImg, img=FlatCombImg,
...@@ -737,10 +747,9 @@ class Chip(FocalPlane): ...@@ -737,10 +747,9 @@ class Chip(FocalPlane):
img_rot=img_rot, img_rot=img_rot,
im_type='FLAT', im_type='FLAT',
pointing_ID=pointing_ID, pointing_ID=pointing_ID,
date_obs=date_obs,
time_obs=time_obs,
output_dir=chip_output.subdir, output_dir=chip_output.subdir,
exptime=self.flat_exptime) exptime=self.flat_exptime,
timestamp=timestamp_obs)
del FlatCombImg, FlatSingle, prnu_img del FlatCombImg, FlatSingle, prnu_img
# flat_img.replaceNegative(replace_value=0) # flat_img.replaceNegative(replace_value=0)
...@@ -783,9 +792,9 @@ class Chip(FocalPlane): ...@@ -783,9 +792,9 @@ class Chip(FocalPlane):
cr_map[cr_map < 0] = 0 cr_map[cr_map < 0] = 0
crmap_gsimg = galsim.Image(cr_map, dtype=np.uint16) crmap_gsimg = galsim.Image(cr_map, dtype=np.uint16)
del cr_map del cr_map
datetime_obs = datetime.utcfromtimestamp(timestamp_obs) # datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
date_obs = datetime_obs.strftime("%y%m%d") # date_obs = datetime_obs.strftime("%y%m%d")
time_obs = datetime_obs.strftime("%H%M%S") # time_obs = datetime_obs.strftime("%H%M%S")
self.outputCal( self.outputCal(
img=crmap_gsimg, img=crmap_gsimg,
ra_cen=ra_cen, ra_cen=ra_cen,
...@@ -793,10 +802,9 @@ class Chip(FocalPlane): ...@@ -793,10 +802,9 @@ class Chip(FocalPlane):
img_rot=img_rot, img_rot=img_rot,
im_type='CRD', im_type='CRD',
pointing_ID=pointing_ID, pointing_ID=pointing_ID,
date_obs=date_obs,
time_obs=time_obs,
output_dir=chip_output.subdir, output_dir=chip_output.subdir,
exptime=self.dark_exptime) exptime=self.dark_exptime,
timestamp=timestamp_obs)
del crmap_gsimg del crmap_gsimg
# Non-Linearity for Dark # Non-Linearity for Dark
...@@ -839,7 +847,7 @@ class Chip(FocalPlane): ...@@ -839,7 +847,7 @@ class Chip(FocalPlane):
readout_noise = galsim.GaussianNoise(rng=rng_readout, sigma=self.read_noise) readout_noise = galsim.GaussianNoise(rng=rng_readout, sigma=self.read_noise)
DarkCombImg.addNoise(readout_noise) DarkCombImg.addNoise(readout_noise)
DarkCombImg = effects.ApplyGainNonUniform16( DarkCombImg, self.gain_channel = effects.ApplyGainNonUniform16(
DarkCombImg, gain=self.gain, DarkCombImg, gain=self.gain,
nsecy = 2, nsecx=8, nsecy = 2, nsecx=8,
seed=SeedGainNonuni+self.chipID, seed=SeedGainNonuni+self.chipID,
...@@ -852,9 +860,9 @@ class Chip(FocalPlane): ...@@ -852,9 +860,9 @@ class Chip(FocalPlane):
DarkCombImg.quantize() DarkCombImg.quantize()
DarkCombImg = galsim.ImageUS(DarkCombImg) DarkCombImg = galsim.ImageUS(DarkCombImg)
# DarkCombImg.write("%s/DarkImg_%s_%s_%s.fits" % (chip_output.subdir, DarkTag, self.chipID, i+1)) # DarkCombImg.write("%s/DarkImg_%s_%s_%s.fits" % (chip_output.subdir, DarkTag, self.chipID, i+1))
datetime_obs = datetime.utcfromtimestamp(timestamp_obs) # datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
date_obs = datetime_obs.strftime("%y%m%d") # date_obs = datetime_obs.strftime("%y%m%d")
time_obs = datetime_obs.strftime("%H%M%S") # time_obs = datetime_obs.strftime("%H%M%S")
timestamp_obs += 10 * 60 timestamp_obs += 10 * 60
self.outputCal( self.outputCal(
img=DarkCombImg, img=DarkCombImg,
...@@ -863,10 +871,9 @@ class Chip(FocalPlane): ...@@ -863,10 +871,9 @@ class Chip(FocalPlane):
img_rot=img_rot, img_rot=img_rot,
im_type='DARK', im_type='DARK',
pointing_ID=pointing_ID, pointing_ID=pointing_ID,
date_obs=date_obs,
time_obs=time_obs,
output_dir=chip_output.subdir, output_dir=chip_output.subdir,
exptime=self.dark_exptime) exptime=self.dark_exptime,
timestamp = timestamp_obs)
del DarkCombImg del DarkCombImg
# img = galsim.ImageUS(img) # img = galsim.ImageUS(img)
......
ObservationSim/Instrument/Chip/Effects.py
View file @ 00a96bcb
...@@ -152,6 +152,7 @@ def ApplyGainNonUniform16(GSImage, gain=1, nsecy = 2, nsecx=8, seed=202102, logg ...@@ -152,6 +152,7 @@ def ApplyGainNonUniform16(GSImage, gain=1, nsecy = 2, nsecx=8, seed=202102, logg
rg = Generator(PCG64(int(seed))) rg = Generator(PCG64(int(seed)))
Random16 = (rg.random(nsecy*nsecx)-0.5)*0.04+1 # sigma~1% Random16 = (rg.random(nsecy*nsecx)-0.5)*0.04+1 # sigma~1%
Gain16 = Random16.reshape((nsecy,nsecx))/gain Gain16 = Random16.reshape((nsecy,nsecx))/gain
gain_array = np.ones(nsecy*nsecx)*gain
if logger is not None: if logger is not None:
msg = str("Gain of 16 channels: " + str(Gain16)) msg = str("Gain of 16 channels: " + str(Gain16))
logger.info(msg) logger.info(msg)
...@@ -163,7 +164,8 @@ def ApplyGainNonUniform16(GSImage, gain=1, nsecy = 2, nsecx=8, seed=202102, logg ...@@ -163,7 +164,8 @@ def ApplyGainNonUniform16(GSImage, gain=1, nsecy = 2, nsecx=8, seed=202102, logg
for rowi in range(nsecy): for rowi in range(nsecy):
for coli in range(nsecx): for coli in range(nsecx):
GSImage.array[rowi*secsize_y:(rowi+1)*secsize_y,coli*secsize_x:(coli+1)*secsize_x] *= Gain16[rowi,coli] 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): def GainsNonUniform16(GSImage, gain=1, nsecy = 2, nsecx=8, seed=202102, logger=None):
......
ObservationSim/Instrument/Filter.py
View file @ 00a96bcb
...@@ -5,6 +5,7 @@ import numpy as np ...@@ -5,6 +5,7 @@ import numpy as np
from ObservationSim.Instrument._util import photonEnergy, calculateLimitMag from ObservationSim.Instrument._util import photonEnergy, calculateLimitMag
from ObservationSim.Instrument.FilterParam import FilterParam from ObservationSim.Instrument.FilterParam import FilterParam
from ObservationSim.Straylight import Straylight
try: try:
import importlib.resources as pkg_resources import importlib.resources as pkg_resources
...@@ -108,6 +109,17 @@ class Filter(object): ...@@ -108,6 +109,17 @@ class Filter(object):
def getSkyNoise(self, exptime, gain=1.): def getSkyNoise(self, exptime, gain=1.):
return self.sky_background * exptime / gain 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
else:
s_pix = sl.calculateStrayLightFilter(filter=self.filter_type.lower())
self.sky_background = s_pix
self.zodical_spec = None
def update_limit_saturation_mags(self, exptime=150., psf_fwhm=0.1969, skyFn='sky_emiss_hubble_50_50_A.dat', chip=None): 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"]: if self.filter_type in ["GI", "GV", "GU"]:
return return
......
ObservationSim/MockObject/SkybackgroundMap.py
View file @ 00a96bcb
...@@ -7,6 +7,7 @@ from astropy.table import Table ...@@ -7,6 +7,7 @@ from astropy.table import Table
from scipy import interpolate from scipy import interpolate
import galsim import galsim
import astropy.constants as cons
import os import os
...@@ -22,7 +23,7 @@ except ImportError: ...@@ -22,7 +23,7 @@ except ImportError:
###calculate sky map by sky SED ###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, 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') # skyMap = np.ones([yLen, xLen], dtype='float32')
# #
# if isAlongY == 1: # if isAlongY == 1:
...@@ -52,6 +53,16 @@ def calculateSkyMap_split_g(skyMap=None, blueLimit=4200, redLimit=6500, skyfn='s ...@@ -52,6 +53,16 @@ def calculateSkyMap_split_g(skyMap=None, blueLimit=4200, redLimit=6500, skyfn='s
# skySpec = np.loadtxt(skyfn) # skySpec = np.loadtxt(skyfn)
spec = Table(np.array([skySpec[:, 0], skySpec[:, 1]]).T, names=('WAVELENGTH', 'FLUX')) 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: if isAlongY == 0:
directParm = 0 directParm = 0
if isAlongY ==1: if isAlongY ==1:
......
ObservationSim/ObservationSim.py
View file @ 00a96bcb
...@@ -107,6 +107,12 @@ class Observation(object): ...@@ -107,6 +107,12 @@ class Observation(object):
chip.img = galsim.ImageF(chip.npix_x, chip.npix_y) chip.img = galsim.ImageF(chip.npix_x, chip.npix_y)
chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin) chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin)
chip.img.wcs = wcs_fp 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": if chip.survey_type == "photometric":
sky_map = None sky_map = None
elif chip.survey_type == "spectroscopic": elif chip.survey_type == "spectroscopic":
...@@ -133,7 +139,8 @@ class Observation(object): ...@@ -133,7 +139,8 @@ class Observation(object):
conf=chip.sls_conf, conf=chip.sls_conf,
pixelSize=chip.pix_scale, pixelSize=chip.pix_scale,
isAlongY=0, 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 del flat_normal
if pointing.pointing_type == 'MS': if pointing.pointing_type == 'MS':
...@@ -176,7 +183,10 @@ class Observation(object): ...@@ -176,7 +183,10 @@ class Observation(object):
pixel_size=chip.pix_size, pixel_size=chip.pix_size,
xcen=chip.x_cen, xcen=chip.x_cen,
ycen=chip.y_cen, ycen=chip.y_cen,
extName='raw') extName='SCI',
timestamp = pointing.timestamp,
exptime = pointing.exp_time,
readoutTime = 40.)
for j in range(self.nobj): for j in range(self.nobj):
...@@ -358,12 +368,36 @@ class Observation(object): ...@@ -358,12 +368,36 @@ class Observation(object):
sat_pos=[pointing.sat_x, pointing.sat_y, pointing.sat_z], sat_pos=[pointing.sat_x, pointing.sat_y, pointing.sat_z],
sat_vel=[pointing.sat_vx, pointing.sat_vy, pointing.sat_vz], sat_vel=[pointing.sat_vx, pointing.sat_vy, pointing.sat_vz],
chip_name=str(chip.chipID).rjust(2, '0')) 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) chip.img = galsim.Image(chip.img.array, dtype=np.uint16)
hdu1 = fits.PrimaryHDU(header=h_prim) hdu1 = fits.PrimaryHDU(header=h_prim)
hdu1.add_checksum() 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 = fits.ImageHDU(chip.img.array, header=h_ext)
hdu2.add_checksum() 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]) hdu1 = fits.HDUList([hdu1, hdu2])
fname = os.path.join(chip_output.subdir, h_prim['FILENAME'] + '.fits') fname = os.path.join(chip_output.subdir, h_prim['FILENAME'] + '.fits')
hdu1.writeto(fname, output_verify='ignore', overwrite=True) hdu1.writeto(fname, output_verify='ignore', overwrite=True)
......
ObservationSim/Straylight/Straylight.py 0 → 100644
View file @ 00a96bcb
import ctypes
import numpy as np
import astropy.constants as cons
from scipy import interpolate
import math
from astropy.table import Table
import astropy.coordinates as coord
from astropy import units as u
import sys
try:
import importlib.resources as pkg_resources
except ImportError:
# Try backported to PY<37 'importlib_resources'
import importlib_resources as pkg_resources
filterPivotWave = {'nuv':2875.5,'u':3629.6,'g':4808.4,'r':6178.2, 'i':7609.0, 'z':9012.9,'y':9627.9}
filterIndex = {'nuv':0,'u':1,'g':2,'r':3, 'i':4, 'z':5,'y':6}
filterCCD = {'nuv':'UV0','u':'UV0','g':'Astro_MB','r':'Astro_MB', 'i':'Basic_NIR', 'z':'Basic_NIR','y':'Basic_NIR'}
bandRange = {'nuv':[2504.0,3230.0],'u':[3190.0,4039.0],'g':[3989.0,5498.0],'r':[5438.0,6956.0], 'i':[6886.0,8469.0], 'z':[8379.0,10855.0],'y':[9217.0, 10900.0], 'GU':[2550, 4000],'GV':[4000, 6200],'GI':[6200,10000]}
# Instrument_dir = '/Users/zhangxin/Work/SlitlessSim/CSST_SIM/CSST_C6/straylight/straylight/Instrument/'
SpecOrder = ['-2','-1','0','1','2']
filterMirrorEff = {'nuv':0.54,'u':0.68,'g':0.8,'r':0.8, 'i':0.8, 'z':0.8,'y':0.8}
def transRaDec2D(ra, dec):
x1 = np.cos(dec / 57.2957795) * np.cos(ra / 57.2957795);
y1 = np.cos(dec / 57.2957795) * np.sin(ra / 57.2957795);
z1 = np.sin(dec / 57.2957795);
return np.array([x1, y1, z1])
def getAngle132(x1 = 0, y1 = 0, z1 = 0, x2 = 0, y2 = 0, z2 = 0, x3 = 0, y3 = 0, z3 = 0):
cosValue = 0;
angle = 0;
x11 = x1-x3;
y11 = y1-y3;
z11 = z1-z3;
x22 = x2-x3;
y22 = y2-y3;
z22 = z2-z3;
tt = np.sqrt((x11*x11 + y11*y11 + z11* z11) * (x22*x22 + y22*y22 + z22*z22));
if(tt==0):
return 0;
cosValue = (x11*x22+y11*y22+z11*z22)/tt;
if (cosValue > 1):
cosValue = 1;
if (cosValue < -1):
cosValue = -1;
angle = math.acos(cosValue);
return angle * 360 / (2 * math.pi);
def calculateAnglePwithEarth(sat = np.array([0,0,0]), pointing = np.array([0,0,0]), sun = np.array([0,0,0])):
modSat = np.sqrt(sat[0]*sat[0] + sat[1]*sat[1]+sat[2]*sat[2])
modPoint = np.sqrt(pointing[0]*pointing[0] + pointing[1]*pointing[1] + pointing[2]*pointing[2])
withLocalZenithAngle = (pointing[0] * sat[0] + pointing[1] * sat[1] + pointing[2] * sat[2]) / (modPoint*modSat)
innerM_sat_sun = sat[0] * sun[0] + sat[1] * sun[1] + sat[2] * sun[2]
cosAngle = innerM_sat_sun / (modSat * cons.au.value/1000)
isInSunSide = 1
if (cosAngle < -0.3385737): #cos109.79
isInSunSide = -1;
elif cosAngle >= -0.3385737 and cosAngle <= 0.3385737:
isInSunSide = 0;
return math.acos(withLocalZenithAngle)*180/math.pi,isInSunSide
# /**
# * *eCoor = ra, *eCoor+1 = dec
# */
def Cartesian2Equatorial(carCoor = np.array([0,0,0])):
eCoor = np.zeros(2)
if (carCoor[0] > 0 and carCoor[1] >= 0):
eCoor[0] = math.atan(carCoor[1] / carCoor[0]) * 360 / (2 * math.pi)
elif (carCoor[0] < 0):
eCoor[0] = (math.atan(carCoor[1] / carCoor[0]) + math.pi) * 360 / (2 * math.pi)
elif (carCoor[0] > 0 and carCoor[1] < 0):
eCoor[0] = (math.atan(carCoor[1] / carCoor[0]) + 2 * math.pi) * 360 / (2 * math.pi)
elif (carCoor[0] == 0 and carCoor[1] < 0):
eCoor[0] = 270
elif (carCoor[0] == 0 and carCoor[1] > 0):
eCoor[0] = 90
eCoor[1] = math.atan(carCoor[2] / np.sqrt(carCoor[0] * carCoor[0] + carCoor[1] * carCoor[1])) * 360 / (2 * math.pi)
return eCoor
class Straylight(object):
def __init__(self, jtime = 2460843., sat_pos = np.array([0,0,0]), pointing_radec = np.array([0,0]), sun_pos = np.array([0,0,0])):
self.jtime = jtime
self.sat = sat_pos
self.sun_pos = sun_pos
self.equator = coord.SkyCoord(pointing_radec[0]*u.degree, pointing_radec[1]*u.degree,frame='icrs')
self.ecliptic = self.equator.transform_to('barycentrictrueecliptic')
self.pointing = transRaDec2D(pointing_radec[0], pointing_radec[1])
platForm = sys.platform
if platForm == 'darwin':
try:
with pkg_resources.files('ObservationSim.Straylight.lib').joinpath('libstraylight.dylib') as dllFn:
self.slcdll = ctypes.CDLL(dllFn)
except AttributeError:
with pkg_resources.path('ObservationSim.Straylight.lib', 'libstraylight.dylib') as dllFn:
self.slcdll = ctypes.CDLL(dllFn)
elif platForm == 'linux':
try:
with pkg_resources.files('ObservationSim.Straylight.so').joinpath('libstraylight.dylib') as dllFn:
self.slcdll = ctypes.CDLL(dllFn)
except AttributeError:
with pkg_resources.path('ObservationSim.Straylight.so', 'libstraylight.dylib') as dllFn:
self.slcdll = ctypes.CDLL(dllFn)
# self.slcdll=ctypes.CDLL('./libstraylight.dylib')
self.slcdll.Calculate.argtypes = [ctypes.c_double, ctypes.POINTER(ctypes.c_double),
ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_double),
ctypes.POINTER(ctypes.c_double), ctypes.c_char_p]
self.slcdll.PointSource.argtypes = [ctypes.c_double, ctypes.POINTER(ctypes.c_double),
ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_double),
ctypes.POINTER(ctypes.c_double), ctypes.c_char_p]
self.slcdll.EarthShine.argtypes = [ctypes.c_double, ctypes.POINTER(ctypes.c_double),
ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_double),
ctypes.POINTER(ctypes.c_double)]
self.slcdll.Zodiacal.argtypes = [ctypes.c_double, ctypes.POINTER(ctypes.c_double),
ctypes.POINTER(ctypes.c_double)]
self.slcdll.ComposeY.argtypes = [ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_double),
ctypes.POINTER(ctypes.c_double)]
self.slcdll.Init.argtypes = [ctypes.c_char_p, ctypes.c_char_p, ctypes.c_char_p, ctypes.c_char_p]
try:
with pkg_resources.files('ObservationSim.Straylight.lib').joinpath('DE405') as tFn:
self.deFn = tFn.as_uri()[7:]
except AttributeError:
with pkg_resources.path('ObservationSim.Straylight.lib', 'DE405') as tFn:
self.deFn = tFn.as_uri()[7:]
try:
with pkg_resources.files('ObservationSim.Straylight.lib').joinpath('PST') as tFn:
self.PSTFn = tFn.as_uri()[7:]
except AttributeError:
with pkg_resources.path('ObservationSim.Straylight.lib', 'PST') as tFn:
self.PSTFn = tFn.as_uri()[7:]
try:
with pkg_resources.files('ObservationSim.Straylight.lib').joinpath('R') as tFn:
self.RFn = tFn.as_uri()[7:]
except AttributeError:
with pkg_resources.path('ObservationSim.Straylight.lib', 'R') as tFn:
self.RFn = tFn.as_uri()[7:]
try:
with pkg_resources.files('ObservationSim.Straylight.lib').joinpath('Zodiacal') as tFn:
self.ZolFn = tFn.as_uri()[7:]
except AttributeError:
with pkg_resources.path('ObservationSim.Straylight.lib', 'Zodiacal') as tFn:
self.ZolFn = tFn.as_uri()[7:]
try:
with pkg_resources.files('ObservationSim.Straylight.lib').joinpath('BrightGaia_with_csst_mag') as tFn:
self.brightStarTabFn = tFn.as_uri()[7:]
except AttributeError:
with pkg_resources.path('ObservationSim.Straylight.lib', 'BrightGaia_with_csst_mag') as tFn:
self.brightStarTabFn = tFn.as_uri()[7:]
print(self.deFn)
self.slcdll.Init(str.encode(self.deFn), str.encode(self.PSTFn), str.encode(self.RFn), str.encode(self.ZolFn))
def getFilterAndCCD_Q(self, filter = 'i'):
try:
with pkg_resources.files('ObservationSim.Instrument.data.ccd').joinpath(filterCCD[filter] + '.txt') as ccd_fn:
q_ccd_f = np.loadtxt(ccd_fn)
except AttributeError:
with pkg_resources.path('ObservationSim.Instrument.data.ccd', filterCCD[filter] + '.txt') as ccd_fn:
q_ccd_f = np.loadtxt(ccd_fn)
try:
with pkg_resources.files('ObservationSim.Instrument.data.filters').joinpath(filter + '.txt') as filter_fn:
q_fil_f = np.loadtxt(filter_fn)
except AttributeError:
with pkg_resources.path('ObservationSim.Instrument.data.filters', filter + '.txt') as filter_fn:
q_fil_f = np.loadtxt(filter_fn)
band_s = 2000
band_e = 11000
q_ccd = np.zeros([q_ccd_f.shape[0]+2,q_ccd_f.shape[1]])
q_ccd[1:-1,:] = q_ccd_f
q_ccd[0] = [band_s,0]
q_ccd[-1] = [band_e,0]
q_fil = np.zeros([q_fil_f.shape[0]+2,q_fil_f.shape[1]])
q_fil[1:-1,:] = q_fil_f
q_fil[0] = [band_s,0]
q_fil[-1] = [band_e,0]
q_fil_i = interpolate.interp1d(q_fil[:,0], q_fil[:,1])
q_ccd_i = interpolate.interp1d(q_ccd[:,0], q_ccd[:,1])
bands = np.arange(bandRange[filter][0], bandRange[filter][1],0.5)
q_ccd_fil = q_fil_i(bands)*q_ccd_i(bands)
return np.trapz(q_ccd_fil, bands)/(bandRange[filter][1]-bandRange[filter][0])
def calculateEarthShineFilter(self, filter = 'i', pixel_size_phy = 10 ):
sat = (ctypes.c_double*3)()
sat[:] = self.sat
ob = (ctypes.c_double*3)()
ob[:]=self.pointing
py1 = (ctypes.c_double*3)()
py2 = (ctypes.c_double*3)()
self.slcdll.ComposeY(ob,py1,py2)
earth_e1 = (ctypes.c_double*7)()
self.slcdll.EarthShine(self.jtime,sat,ob,py1,earth_e1)
earth_e2 = (ctypes.c_double*7)()
self.slcdll.EarthShine(self.jtime,sat,ob,py2,earth_e2)
band_earth_e1 = earth_e1[:][filterIndex[filter]]
band_earth_e2 = earth_e2[:][filterIndex[filter]]
q=self.getFilterAndCCD_Q(filter=filter)
p_lambda = filterPivotWave[filter]
c = cons.c.value
h = cons.h.value
pix_earth_e1 = band_earth_e1/(h*c/(p_lambda*1e-10))*pixel_size_phy*1e-6*pixel_size_phy*1e-6*q
pix_earth_e2 = band_earth_e2/(h*c/(p_lambda*1e-10))*pixel_size_phy*1e-6*pixel_size_phy*1e-6*q
if pix_earth_e1< pix_earth_e2:
return pix_earth_e1, py1[:]
else:
return pix_earth_e2, py2[:]
"""
calculate zodiacal call c++ program, seems to have some problem
"""
def calculateZodiacalFilter1(self, filter = 'i', pixel_size_phy = 10 ):
sat = (ctypes.c_double*3)()
sat[:] = self.sat
ob = (ctypes.c_double*3)()
ob[:]=self.pointing
zodical_e = (ctypes.c_double*7)()
self.slcdll.Zodiacal(self.jtime,ob,zodical_e)
ob1 = (ctypes.c_double*2)()
ob1[:] = np.array([self.ecliptic.lon.value, self.ecliptic.lat.value])
zodical_e1 = (ctypes.c_double*7)()
self.slcdll.Zodiacal1(ob1,zodical_e1)
band_zodical_e = zodical_e[:][filterIndex[filter]]
q=self.getFilterAndCCD_Q(filter=filter)
p_lambda = filterPivotWave[filter]
c = cons.c.value
h = cons.h.value
pix_zodical_e = band_zodical_e/(h*c/(p_lambda*1e-10))*pixel_size_phy*1e-6*pixel_size_phy*1e-6*q
return pix_zodical_e, band_zodical_e
"""
calculate zodiacal use python
"""
def calculateZodiacalFilter2(self,filter = 'i', aper = 2, pixelsize = 0.074, sun_pos = np.array([0,0,0])):
spec, v_mag = self.calculateZodicalSpec(longitude = self.ecliptic.lon.value, latitude = self.ecliptic.lat.value, sun_pos = sun_pos)
# spec = self.calculateZodicalSpec(longitude = lon, latitude = lat)
try:
with pkg_resources.files('ObservationSim.Instrument.data.throughputs').joinpath(filter + '_throughput.txt') as throughputFn:
throughput = np.loadtxt(throughputFn)
except AttributeError:
with pkg_resources.path('ObservationSim.Instrument.data.throughputs', filter + '_throughput.txt') as throughputFn:
throughput = np.loadtxt(throughputFn)
deltL = 0.5
lamb = np.arange(bandRange[filter][0], bandRange[filter][1], deltL)
speci = interpolate.interp1d(spec['WAVELENGTH'], spec['FLUX'])
y = speci(lamb)
# erg/s/cm2/A --> photo/s/m2/A
flux = y * lamb / (cons.h.value * cons.c.value) * 1e-13
throughput_i = interpolate.interp1d(throughput[:, 0], throughput[:, 1])
throughput_ = throughput_i(lamb)
sky_pix = np.trapz(flux*throughput_, lamb) * math.pi * aper*aper/4 * pixelsize * pixelsize
# sky_pix_e = np.trapz(y, lamb) * math.pi * aper*aper/4 * pixelsize * pixelsize/(10*10*1e-6*1e-6)*1e-7*1e4
return sky_pix, v_mag#, sky_pix_e
def calculateStarLightFilter(self, filter = 'i', pointYaxis = np.array([1,1,1]), pixel_size_phy = 10 ):
sat = (ctypes.c_double*3)()
sat[:] = self.sat
ob = (ctypes.c_double*3)()
ob[:]=self.pointing
py = (ctypes.c_double*3)()
py[:] = pointYaxis
q=self.getFilterAndCCD_Q(filter=filter)
p_lambda = filterPivotWave[filter]
c = cons.c.value
h = cons.h.value
star_e1 = (ctypes.c_double*7)()
self.slcdll.PointSource(self.jtime,sat,ob,py,star_e1, str.encode(self.brightStarTabFn))
band_star_e1 = star_e1[:][filterIndex[filter]]
pix_star_e1 = band_star_e1/(h*c/(p_lambda*1e-10))*pixel_size_phy*1e-6*pixel_size_phy*1e-6*q
return pix_star_e1
def calculateEarthshineGrating(self, grating = 'GU', pixel_size_phy = 10, normFilter = 'g', aper = 2, pixelsize = 0.074):
sat = (ctypes.c_double*3)()
sat[:] = self.sat
ob = (ctypes.c_double*3)()
ob[:]=self.pointing
py1 = (ctypes.c_double*3)()
py2 = (ctypes.c_double*3)()
self.slcdll.ComposeY(ob,py1,py2)
earth_e1 = (ctypes.c_double*7)()
self.slcdll.EarthShine(self.jtime,sat,ob,py1,earth_e1)
earth_e2 = (ctypes.c_double*7)()
self.slcdll.EarthShine(self.jtime,sat,ob,py2,earth_e2)
# zodical_e = (ctypes.c_double*7)()
# self.slcdll.Zodiacal(self.jtime,ob,zodical_e)
band_earth_e1 = earth_e1[:][filterIndex[normFilter]]
band_earth_e2 = earth_e2[:][filterIndex[normFilter]]
band_earth_e = band_earth_e2
py = py2[:]
if band_earth_e1<band_earth_e2:
band_earth_e = band_earth_e1
py = py1[:]
# band_earth_e = np.min([band_earth_e1, band_earth_e2])
# band_earth_e1 = 0
# band_earth_e2 = 0
# band_zodical_e = zodical_e[:][filterIndex[normFilter]]
q=self.getFilterAndCCD_Q(filter=normFilter)
p_lambda = filterPivotWave[normFilter]
c = cons.c.value
h = cons.h.value
pix_earth_e = band_earth_e/(h*c/(p_lambda*1e-10))*pixel_size_phy*1e-6*pixel_size_phy*1e-6*q
# pix_earth_e2 = band_earth_e2/(h*c/(p_lambda*1e-10))*pixel_size_phy*1e-6*pixel_size_phy*1e-6*q
# pix_zodical_e = band_zodical_e/(h*c/(p_lambda*1e-10))*pixel_size_phy*1e-6*pixel_size_phy*1e-6*q
# pix_earth_e = np.min([pix_earth_e1, pix_earth_e2])
earthshine_v, earthshine_spec = self.calculatEarthshineByHSTSpec(filter = normFilter, aper = aper, pixelsize = pixelsize)
lamb_earth = earthshine_spec['WAVELENGTH']
flux_earth = earthshine_spec['FLUX']*pix_earth_e/earthshine_v
# print(pix_earth_e,earthshine_v)
earth_v_grating = 0
for s_order in SpecOrder:
try:
with pkg_resources.files('ObservationSim.Instrument.data.sls_conf').joinpath(
grating + '.Throughput.' + s_order + 'st.fits') as thpFn:
thp_ = Table.read(thpFn)
except AttributeError:
with pkg_resources.path('ObservationSim.Instrument.data.sls_conf',
grating + '.Throughput.' + s_order + 'st.fits') as thpFn:
thp_ = Table.read(thpFn)
thpFn_i = interpolate.interp1d(thp_['WAVELENGTH'], thp_['SENSITIVITY'])
thp = thpFn_i(lamb_earth)
beamsEarth = np.trapz(flux_earth*thp,lamb_earth)* math.pi*aper*aper/4 * pixelsize * pixelsize
earth_v_grating = earth_v_grating + beamsEarth
# print(beamsEarth)
# print(earthshine_v, pix_earth_e, earth_v_grating)
return earth_v_grating, py
def calculateStarLightGrating(self, grating = 'GU', pointYaxis = np.array([1,1,1]), pixel_size_phy = 10 ):
sat = (ctypes.c_double*3)()
sat[:] = self.sat
ob = (ctypes.c_double*3)()
ob[:]=self.pointing
py = (ctypes.c_double*3)()
py[:] = pointYaxis
# q=self.getFilterAndCCD_Q(filter=filter)
# p_lambda = filterPivotWave[filter]
c = cons.c.value
h = cons.h.value
star_e1 = (ctypes.c_double*7)()
self.slcdll.PointSource(self.jtime,sat,ob,py,star_e1, str.encode(self.brightStarTabFn))
filterPivotWaveList = np.zeros(7)
bandRangeList = np.zeros(7)
filterMirrorEffList = np.zeros(7)
filterNameList = list(filterPivotWave.keys())
for i in np.arange(7):
filterPivotWaveList[i] = filterPivotWave[filterNameList[i]]
filterMirrorEffList[i] = filterMirrorEff[filterNameList[i]]
brange = bandRange[filterNameList[i]]
bandRangeList[i] = brange[1] - brange[0]
filterFlux_lamb = star_e1[:]/bandRangeList/filterMirrorEffList/(h*c/(filterPivotWaveList*1e-10))
filterFlux_lambi = interpolate.interp1d(filterPivotWaveList,filterFlux_lamb,fill_value="extrapolate")
lamb_g = np.arange(bandRange[grating][0], bandRange[grating][1],1)
flux_g = filterFlux_lambi(lamb_g)
# flux_total_g = np.trapz(flux_g,lamb_g)
starLight_grating = 0
for s_order in SpecOrder:
try:
with pkg_resources.files('ObservationSim.Instrument.data.sls_conf').joinpath(
grating + '.Throughput.' + s_order + 'st.fits') as thpFn:
thp_ = Table.read(thpFn)
except AttributeError:
with pkg_resources.path('ObservationSim.Instrument.data.sls_conf',
grating + '.Throughput.' + s_order + 'st.fits') as thpFn:
thp_ = Table.read(thpFn)
thpFn_i = interpolate.interp1d(thp_['WAVELENGTH'], thp_['SENSITIVITY'])
thp = thpFn_i(lamb_g)
beamsZol = np.trapz(flux_g*thp,lamb_g)*pixel_size_phy*1e-6*pixel_size_phy*1e-6
starLight_grating = starLight_grating + beamsZol
# print(beamsZol)
# band_star_e1 = star_e1[:][filterIndex[filter]]
# pix_star_e1 = band_star_e1/(h*c/(p_lambda*1e-10))*pixel_size_phy*1e-6*pixel_size_phy*1e-6*q
return starLight_grating
def calculatEarthshineByHSTSpec(self, filter = 'g', aper = 2, pixelsize = 0.074, s = 2000, e = 11000):
try:
with pkg_resources.files('ObservationSim.Straylight.data.sky').joinpath('earthShine.dat') as specFn:
spec = np.loadtxt(specFn)
except AttributeError:
with pkg_resources.path('ObservationSim.Straylight.data.sky',
'earthShine.dat') as specFn:
spec = np.loadtxt(specFn)
try:
with pkg_resources.files('ObservationSim.Instrument.data.throughputs').joinpath(filter + '_throughput.txt') as throughputFn:
throughput = np.loadtxt(throughputFn)
except AttributeError:
with pkg_resources.path('ObservationSim.Instrument.data.throughputs', filter + '_throughput.txt') as throughputFn:
throughput = np.loadtxt(throughputFn)
deltL = 0.5
lamb = np.arange(bandRange[filter][0], bandRange[filter][1], deltL)
speci = interpolate.interp1d(spec[:, 0], spec[:, 1])
y = speci(lamb)
# erg/s/cm2/A --> photo/s/m2/A
flux = y * lamb / (cons.h.value * cons.c.value) * 1e-13
throughput_i = interpolate.interp1d(throughput[:, 0], throughput[:, 1])
throughput_ = throughput_i(lamb)
sky_pix = np.trapz(flux*throughput_, lamb) * math.pi * aper*aper/4 * pixelsize * pixelsize
lamb = np.arange(s, e, deltL)
speci = interpolate.interp1d(spec[:, 0], spec[:, 1])
y = speci(lamb)
# erg/s/cm2/A --> photo/s/m2/A
flux = y * lamb / (cons.h.value * cons.c.value) * 1e-13
return sky_pix, Table(np.array([lamb, flux]).T,names=('WAVELENGTH', 'FLUX'))
def calculateZodicalSpec(self,longitude = 50, latitude = 60, sun_pos = np.array([0,0,0])):
from scipy.interpolate import interp2d
from scipy.interpolate import griddata
try:
with pkg_resources.files('ObservationSim.Straylight.data').joinpath('Zodiacal_map1.dat') as z_map_fn:
ZL = np.loadtxt(z_map_fn)
except AttributeError:
with pkg_resources.path('ObservationSim.Straylight.data',
'Zodiacal_map1.dat') as z_map_fn:
ZL = np.loadtxt(z_map_fn)
# zl_sh = ZL.shape
# x = np.arange(0,zl_sh[1],1)
# y = np.arange(0,zl_sh[0],1)
x = ZL[0,1:]
y = ZL[1:,0]
X,Y = np.meshgrid(x,y)
# f_sur = interp2d(X,Y,ZL,kind='linear')
sun_radec = Cartesian2Equatorial(sun_pos)
sun_eclip = coord.SkyCoord(sun_radec[0]*u.degree, sun_radec[1]*u.degree,frame='icrs')
sun_equtor = sun_eclip.transform_to('barycentrictrueecliptic')
longitude = longitude - (sun_equtor.lon*u.degree).value
longitude = np.abs(longitude)
# print((sun_equtor.lon*u.degree).value)
if (longitude > 180):
longitude = 360 - longitude
latitude = np.abs(latitude)
lo = longitude
la = latitude
zl = griddata((X.flatten(),Y.flatten()),ZL[1:,1:].flatten(),(la,lo), method='cubic').min()
zl = zl*(math.pi*math.pi)/(180*180)/(3600*3600)*1e-4*1e7*1e-8*1e-4
# print(zl , '\n')
try:
with pkg_resources.files('ObservationSim.Straylight.data.sky').joinpath('zodiacal.dat') as zodical_fn:
spec = np.loadtxt(zodical_fn)
except AttributeError:
with pkg_resources.path('ObservationSim.Straylight.data.sky',
'zodiacal.dat') as zodical_fn:
spec = np.loadtxt(zodical_fn)
speci = interpolate.interp1d(spec[:, 0], spec[:, 1])
flux5000 = speci(5000)
f_ration = zl/flux5000
v_mag = np.log10(f_ration)*(-2.5)+22.1
# print("factor:", v_mag, lo, la)
return Table(np.array([spec[:,0], spec[:,1]*f_ration]).T,names=('WAVELENGTH', 'FLUX')), v_mag
def calculateStrayLightFilter(self, filter = 'i', pixel_size_phy = 10, pixel_scale = 0.074):
e1,py = self.calculateEarthShineFilter(filter = filter, pixel_size_phy = pixel_size_phy)
e2, _ = self.calculateZodiacalFilter2(filter = 'i', sun_pos=self.sun_pos, pixelsize = pixel_scale)
e3 = self.calculateStarLightFilter(filter = 'i',pointYaxis = py, pixel_size_phy = pixel_size_phy)
return e1+e2+e3
def calculateStrayLightGrating(self, grating = 'GI', pixel_size_phy = 10, normFilter_es = 'g'):
e1,py = self.calculateEarthshineGrating(grating = grating, pixel_size_phy = pixel_size_phy, normFilter = normFilter_es)
e2 = self.calculateStarLightGrating(grating = grating, pointYaxis = py)
spec, _ = self.calculateZodicalSpec(longitude = self.ecliptic.lon.value, latitude = self.ecliptic.lat.value, sun_pos = self.sun_pos)
return e1+e2, spec
def testZodiacal(lon = 285.04312526255366, lat = 30.):
c_eclip = coord.SkyCoord(lon*u.degree, lat*u.degree,frame='barycentrictrueecliptic')
c_equtor = c_eclip.transform_to('icrs')
sl = Straylight(jtime = 2459767.00354975, sat = np.array([]), radec = np.array([(c_equtor.ra*u.degree).value, (c_equtor.dec*u.degree).value]))
e_zol, v_mag = sl.calculateZodiacalFilter2(filter = 'i', sun_pos=np.array([-3.70939436e+07, 1.35334903e+08, 5.86673104e+07]))
print(e_zol)
# ju=2.4608437604166665e+06
# sat = (ctypes.c_double*3)()
# sat[:] = np.array([5873.752, -1642.066, 2896.744])
# ob = (ctypes.c_double*3)()
# ob[:]=np.array([0.445256,0.76061,-0.47246])
# sl = StrayLight(jtime = ju, sat = np.array([5873.752, -1642.066, 2896.744]), pointing = np.array([-0.445256,-0.76061,0.47246]))
# fn = '/Users/zhangxin/Work/SurveyPlan/point/csst_survey_sim_20211028/E17.5_b17.5_beta_11.6_opt_transtime_1_CMG_1_dp_2_0.25_da_10_Texp_1.5_DEC60_500_0.1_800_1000_+5deg.dat'
#
# surveylist = np.loadtxt(fn)
# sky_pix = np.zeros([surveylist.shape[0],7])
#
#
# i = 693438
# c_eclip = coord.SkyCoord(surveylist[:,2]*u.degree, surveylist[:,1]*u.degree,frame='barycentrictrueecliptic')
# c_equtor = c_eclip.transform_to('icrs')
#
#
# # pointing = transRaDec2D((c_equtor[i].ra*u.degree).value, (c_equtor[i].dec*u.degree).value)
# # # print(ju, pointing, surveylist[i,3:9])
# # ju = surveylist[i,0]
# # sl = StrayLight(jtime = ju, sat = surveylist[i,3:6], pointing = pointing)
# # sl.caculateStrayLightGrating(grating = 'GI', pixel_size_phy = 10, normFilter = 'g')
#
# for i in np.arange(surveylist.shape[0]):
# print(i)
# if i > 300:
# break
# # if i != 300:
# # continue
# # if i != 693438:
# # continue
# ju = surveylist[i,0]
# pointing = transRaDec2D((c_equtor[i].ra*u.degree).value, (c_equtor[i].dec*u.degree).value)
# # print(ju, pointing, surveylist[i,3:9])
# sl = StrayLight(jtime = ju, sat = surveylist[i,3:6], radec = np.array([(c_equtor[i].ra*u.degree).value, (c_equtor[i].dec*u.degree).value]))
# # strayl_i,s_zoldical ,s_earth, s_earth1 = sl.caculateStrayLightFilter(filter = 'i')
# # print(i,strayl_i,s_zoldical,s_earth, s_earth1)
# p_cart= transRaDec2D((c_equtor[i].ra*u.degree).value, (c_equtor[i].dec*u.degree).value)
# sky_pix[i,6] = getAngle132(x1 = surveylist[i,6], y1 = surveylist[i,7], z1 = surveylist[i,8], x2 = p_cart[0], y2 = p_cart[1], z2 = p_cart[2], x3 = 0, y3 = 0, z3 = 0)
#
# earthZenithAngle,isInSunSide = calculateAnglePwithEarth(sat = surveylist[i,3:6], pointing = pointing, sun = surveylist[i,6:9])
# sky_pix[i,4] = earthZenithAngle
# sky_pix[i,5] = isInSunSide
#
# e1_,py = sl.caculateEarthShineFilter(filter = 'i')
# # e2, e2_ = sl.calculateZodiacalFilter1(filter = 'i')
# e3, v_mag = sl.calculateZodiacalFilter2(filter = 'i', sun_pos=surveylist[i,6:9])
# # e4 = sl.caculateStarLightFilter(filter = 'i',pointYaxis = py)
# # e4 = 0
#
# e1,py = sl.caculateEarthshineGrating(grating = 'GI', pixel_size_phy = 10, normFilter = 'g')
#
# # e2 = sl.caculateStarLightGrating(grating = 'GV', pointYaxis = py)
# e2 = sl.caculateStarLightGrating(grating = 'GI', pointYaxis = py)
# e4 = sl.caculateStarLightFilter(filter = 'i',pointYaxis = py)
#
# e5=sl.caculateStrayLightFilter(filter = 'i', pixel_size_phy = 10, pixel_scale = 0.074, sun_pos = surveylist[i,6:9])
# e6,_=sl.caculateStrayLightGrating(grating = 'GI', normFilter_es = 'g', sun_pos = surveylist[i,6:9])
#
# sky_pix[i,0] = e1
# sky_pix[i,1] = e2
# sky_pix[i,2] = e3
# sky_pix[i,3] = e4
# print(e1+e2,e1_+e3+e4,e5,e6)
#
# # print(e1,e2,e3,e4)
ObservationSim/Straylight/__init__.py 0 → 100755
View file @ 00a96bcb
from .Straylight import Straylight
ObservationSim/Straylight/data/Zodiacal_map1.dat 0 → 100644
View file @ 00a96bcb
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
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