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

Merge branch 'develop' into wcs_test

parents 8649bb95 78bf42f6
Show whitespace changes
Inline Side-by-side
ObservationSim/Config/Header/ImageHeader.py
View file @ adb90dd3
......@@ -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
......@@ -382,7 +386,7 @@ def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec
h_prim['EXPTIME'] = exptime
# Define file types
file_type = {'SCI':'SCIE', 'BIAS':'BIAS', 'DARK':'DARK', 'FLAT':'FLAT', 'CRS':'cosmic_ray', 'CRD':'cosmic_ray','CALS':'CALS','CALF':'CALF'}
file_type = {'SCI':'SCIE', 'BIAS':'BIAS', 'DARK':'DARK', 'FLAT':'FLAT', 'CRS':'CRS', 'CRD':'CRD','CALS':'CALS','CALF':'CALF'}
h_prim['FILETYPE'] = file_type[im_type]
co = coord.SkyCoord(ra, dec, unit='deg')
......@@ -418,7 +422,7 @@ def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec
file_end_time = end_time_str[0:4] + end_time_str[5:7]+end_time_str[8:10] + end_time_str[11:13] + end_time_str[14:16] + end_time_str[17:19]
# h_prim['FILENAME'] = 'CSST_MSC_MS_' + im_type + '_' + file_start_time + '_' + file_end_time + '_' + OBS_id + '_' + CCDID[
# k - 1].rjust(2, '0') + '_L0_V01'
h_prim['FILENAME'] = 'CSST_MSC_MS_' + im_type + '_' + file_start_time + '_' + file_end_time + '_' + OBS_id + '_' + chip_name + '_L0_V01'
h_prim['FILENAME'] = 'CSST_MSC_MS_' + file_type[im_type] + '_' + file_start_time + '_' + file_end_time + '_' + OBS_id + '_' + chip_name + '_L0_V01'
h_prim['POSI0_X'] = sat_pos[0]
......@@ -434,12 +438,15 @@ def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec
# Get version of CSSTSim Package
from pkg_resources import get_distribution
# h_prim['SIM_VER'] = (get_distribution("CSSTSim").version, "Version of CSST MSC simulation software")
h_prim['FITSCREA'] = get_distribution("CSSTSim").version
currentDateAndTime = datetime.now()
compute_name = platform.node()
h_prim['FITSCREA'] = get_distribution("CSSTSim").version +'_' + currentDateAndTime.strftime("%Y%m%d") + '_' +compute_name
h_prim['EPOCH'] = round((Time(h_prim['EXPSTART'], format='mjd', scale='tcb')).jyear, 1)
return h_prim
def generateExtensionHeader(chip, xlen = 9216, ylen = 9232,ra = 60, dec = -40, pa = -23.433, gain = 1.0, readout = 5.0, dark = 0.02, saturation=90000, pixel_scale = 0.074, pixel_size=1e-2,
extName='SCI', row_num = None, col_num = None, xcen=None, ycen=None):
extName='SCI', row_num = None, col_num = None, xcen=None, ycen=None, timestamp = 1621915200,exptime = 150., readoutTime = 40.):
e_header_fn = os.path.split(os.path.realpath(__file__))[0] + '/extension_header.header'
f = open(os.path.split(os.path.realpath(__file__))[0] + '/filter.lst')
......@@ -461,47 +468,60 @@ def generateExtensionHeader(chip, xlen = 9216, ylen = 9232,ra = 60, dec = -40, p
# h_ext['CCDCHIP'] = CCDID[k - 1].rjust(2, '0')
# h_ext['CCDLABEL'] = filters[k-1] + '-' + filterID[k-1]
# h_ext['FILTER'] = filters[k-1]
h_ext['CCDCHIP'] = chip.chip_name
h_ext['CCDLABEL'] = chip.filter_type + '-' + str(chip.chipID).rjust(2, '0')
h_ext['CHIPID'] = str(chip.chipID).rjust(2, '0')
h_ext['CHIPLAB'] = chip.chip_name
h_ext['FILTER'] = chip.filter_type
h_ext['NAXIS1'] = xlen
h_ext['NAXIS2'] = ylen
h_ext['EXTNAME'] = extName
h_ext['GAIN1'] = gain
h_ext['GAIN2'] = gain
h_ext['GAIN3'] = gain
h_ext['GAIN4'] = gain
h_ext['GAIN5'] = gain
h_ext['GAIN6'] = gain
h_ext['GAIN7'] = gain
h_ext['GAIN8'] = gain
h_ext['GAIN9'] = gain
h_ext['GAIN10'] = gain
h_ext['GAIN11'] = gain
h_ext['GAIN12'] = gain
h_ext['GAIN13'] = gain
h_ext['GAIN14'] = gain
h_ext['GAIN15'] = gain
h_ext['GAIN16'] = gain
h_ext['RDNOIS1'] = readout
h_ext['RDNOIS2'] = readout
h_ext['RDNOIS3'] = readout
h_ext['RDNOIS4'] = readout
h_ext['RDNOIS5'] = readout
h_ext['RDNOIS6'] = readout
h_ext['RDNOIS7'] = readout
h_ext['RDNOIS8'] = readout
h_ext['RDNOIS9'] = readout
h_ext['RDNOIS10'] = readout
h_ext['RDNOIS11'] = readout
h_ext['RDNOIS12'] = readout
h_ext['RDNOIS13'] = readout
h_ext['RDNOIS14'] = readout
h_ext['RDNOIS15'] = readout
h_ext['RDNOIS16'] = readout
h_ext['GAIN01'] = chip.gain_channel[0]
h_ext['GAIN02'] = chip.gain_channel[1]
h_ext['GAIN03'] = chip.gain_channel[2]
h_ext['GAIN04'] = chip.gain_channel[3]
h_ext['GAIN05'] = chip.gain_channel[4]
h_ext['GAIN06'] = chip.gain_channel[5]
h_ext['GAIN07'] = chip.gain_channel[6]
h_ext['GAIN08'] = chip.gain_channel[7]
h_ext['GAIN09'] = chip.gain_channel[8]
h_ext['GAIN10'] = chip.gain_channel[9]
h_ext['GAIN11'] = chip.gain_channel[10]
h_ext['GAIN12'] = chip.gain_channel[11]
h_ext['GAIN13'] = chip.gain_channel[12]
h_ext['GAIN14'] = chip.gain_channel[13]
h_ext['GAIN15'] = chip.gain_channel[14]
h_ext['GAIN16'] = chip.gain_channel[15]
h_ext['RON01'] = readout
h_ext['RON02'] = readout
h_ext['RON03'] = readout
h_ext['RON04'] = readout
h_ext['RON05'] = readout
h_ext['RON06'] = readout
h_ext['RON07'] = readout
h_ext['RON08'] = readout
h_ext['RON09'] = readout
h_ext['RON10'] = readout
h_ext['RON11'] = readout
h_ext['RON12'] = readout
h_ext['RON13'] = readout
h_ext['RON14'] = readout
h_ext['RON15'] = readout
h_ext['RON16'] = readout
h_ext['PIXSCAL1'] = pixel_scale
h_ext['PIXSCAL2'] = pixel_scale
h_ext['EXPTIME'] = exptime
h_ext['DARKTIME'] = exptime + 2
datetime_obs = datetime.utcfromtimestamp(timestamp)
tstart = Time(datetime_obs)
tstart_read = Time(tstart.mjd + exptime / 86400., format="mjd")
tend_read = Time(tstart.mjd + (exptime + readoutTime) / 86400., format="mjd")
t_s1 = str(tstart_read.datetime).split()
h_ext['READT0'] = t_s1[0]+'T'+t_s1[1]
t_s2 = str(tend_read.datetime).split()
h_ext['READT1'] = t_s2[0] + 'T' + t_s2[1]
# h_ext['POS_ANG'] = pa
header_wcs = WCS_def(xlen=xlen, ylen=ylen, gapy=898.0, gapx1=534, gapx2=1309, ra_ref=ra, dec_ref=dec, pa=pa, pixel_scale=pixel_scale, pixel_size=pixel_size,
......@@ -515,8 +535,8 @@ def generateExtensionHeader(chip, xlen = 9216, ylen = 9232,ra = 60, dec = -40, p
h_ext['CD1_2'] = header_wcs['CD1_2']
h_ext['CD2_1'] = header_wcs['CD2_1']
h_ext['CD2_2'] = header_wcs['CD2_2']
h_ext['EQUINOX'] = header_wcs['EQUINOX']
h_ext['WCSDIM'] = header_wcs['WCSDIM']
# h_ext['EQUINOX'] = header_wcs['EQUINOX']
# h_ext['WCSDIM'] = header_wcs['WCSDIM']
h_ext['CTYPE1'] = header_wcs['CTYPE1']
h_ext['CTYPE2'] = header_wcs['CTYPE2']
......
ObservationSim/Config/Header/extension_header.header
View file @ adb90dd3
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 @ adb90dd3
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 @ adb90dd3
......@@ -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 @ adb90dd3
......@@ -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
......@@ -305,7 +307,10 @@ class Chip(FocalPlane):
if self.chipID == 30: confFile = ['CSST_GI8.conf', 'CSST_GI7.conf']
return confFile
def generateHeader(self, ra_cen, dec_cen, img_rot, im_type, pointing_ID, date_obs, time_obs, exptime=150.):
def generateHeader(self, ra_cen, dec_cen, img_rot, im_type, pointing_ID, exptime=150., timestamp = 1621915200):
datetime_obs = datetime.utcfromtimestamp(timestamp)
date_obs = datetime_obs.strftime("%y%m%d")
time_obs = datetime_obs.strftime("%H%M%S")
h_prim = generatePrimaryHeader(
xlen=self.npix_x,
ylen=self.npix_y,
......@@ -334,23 +339,30 @@ class Chip(FocalPlane):
pixel_size=self.pix_size,
xcen=self.x_cen,
ycen=self.y_cen,
extName='raw')
extName='SCI',
timestamp = timestamp,
exptime = exptime,
readoutTime = 40.)
return h_prim, h_ext
def outputCal(self, img, ra_cen, dec_cen, img_rot, im_type, pointing_ID, date_obs, time_obs, output_dir, exptime=150.):
def outputCal(self, img, ra_cen, dec_cen, img_rot, im_type, pointing_ID, output_dir, exptime=150., timestamp = 1621915200):
h_prim, h_ext = self.generateHeader(
ra_cen=ra_cen,
dec_cen=dec_cen,
img_rot=img_rot,
im_type=im_type,
pointing_ID=pointing_ID,
date_obs=date_obs,
time_obs=time_obs,
exptime=exptime)
exptime=exptime,
timestamp = timestamp)
hdu1 = fits.PrimaryHDU(header=h_prim)
hdu1.add_checksum()
hdu1.header.comments['CHECKSUM'] = 'HDU checksum'
hdu1.header.comments['DATASUM'] = 'data unit checksum'
hdu2 = fits.ImageHDU(img.array, header=h_ext)
hdu2.add_checksum()
hdu2.header.comments['XTENSION'] = 'extension type'
hdu2.header.comments['CHECKSUM'] = 'HDU checksum'
hdu2.header.comments['DATASUM'] = 'data unit checksum'
hdu1 = fits.HDUList([hdu1, hdu2])
fname = os.path.join(output_dir, h_prim['FILENAME']+'.fits')
hdu1.writeto(fname, output_verify='ignore', overwrite=True)
......@@ -457,9 +469,9 @@ class Chip(FocalPlane):
del cr_map
# crmap_gsimg.write("%s/CosmicRay_%s_1.fits" % (chip_output.subdir, self.chipID))
# crmap_gsimg.write("%s/CosmicRay_%s.fits" % (chip_output.subdir, self.chipID))
datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
date_obs = datetime_obs.strftime("%y%m%d")
time_obs = datetime_obs.strftime("%H%M%S")
# datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
# date_obs = datetime_obs.strftime("%y%m%d")
# time_obs = datetime_obs.strftime("%H%M%S")
self.outputCal(
img=crmap_gsimg,
ra_cen=ra_cen,
......@@ -467,10 +479,9 @@ class Chip(FocalPlane):
img_rot=img_rot,
im_type='CRS',
pointing_ID=pointing_ID,
date_obs=date_obs,
time_obs=time_obs,
output_dir=chip_output.subdir,
exptime=exptime)
exptime=exptime,
timestamp=timestamp_obs)
del crmap_gsimg
# Apply PRNU effect and output PRNU flat file:
......@@ -556,7 +567,7 @@ class Chip(FocalPlane):
else:
print(" Applying Gain (and 16 channel non-uniformity) & Quantization", flush=True)
if config["ins_effects"]["gain_16channel"] == True:
img = effects.ApplyGainNonUniform16(
img, self.gain_channel = effects.ApplyGainNonUniform16(
img, gain=self.gain,
nsecy = 2, nsecx=8,
seed=SeedGainNonuni+self.chipID,
......@@ -611,7 +622,7 @@ class Chip(FocalPlane):
if config["ins_effects"]["add_badcolumns"] == True:
BiasCombImg = effects.BadColumns(BiasCombImg-float(self.bias_level)+5, seed=SeedBadColumns, chipid=self.chipID, logger=self.logger) + float(self.bias_level)-5
BiasCombImg = effects.ApplyGainNonUniform16(BiasCombImg, gain=self.gain,
BiasCombImg, self.gain_channel = effects.ApplyGainNonUniform16(BiasCombImg, gain=self.gain,
nsecy = 2, nsecx=8,
seed=SeedGainNonuni+self.chipID,
logger=self.logger)
......@@ -623,9 +634,9 @@ class Chip(FocalPlane):
BiasCombImg.quantize()
BiasCombImg = galsim.ImageUS(BiasCombImg)
# BiasCombImg.write("%s/BiasImg_%s_%s_%s.fits" % (chip_output.subdir, BiasTag, self.chipID, i+1))
datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
date_obs = datetime_obs.strftime("%y%m%d")
time_obs = datetime_obs.strftime("%H%M%S")
# datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
# date_obs = datetime_obs.strftime("%y%m%d")
# time_obs = datetime_obs.strftime("%H%M%S")
timestamp_obs += 10 * 60
self.outputCal(
img=BiasCombImg,
......@@ -634,10 +645,9 @@ class Chip(FocalPlane):
img_rot=img_rot,
im_type='BIAS',
pointing_ID=pointing_ID,
date_obs=date_obs,
time_obs=time_obs,
output_dir=chip_output.subdir,
exptime=0.0)
exptime=0.0,
timestamp=timestamp_obs)
del BiasCombImg
# Export combined (ncombine, Vignetting + PRNU) & single vignetting flat-field file
......@@ -718,7 +728,7 @@ class Chip(FocalPlane):
readout_noise = galsim.GaussianNoise(rng=rng_readout, sigma=self.read_noise)
FlatCombImg.addNoise(readout_noise)
FlatCombImg = effects.ApplyGainNonUniform16(FlatCombImg, gain=self.gain,
FlatCombImg, self.gain_channel = effects.ApplyGainNonUniform16(FlatCombImg, gain=self.gain,
nsecy = 2, nsecx=8,
seed=SeedGainNonuni+self.chipID,
logger=self.logger)
......@@ -727,9 +737,9 @@ class Chip(FocalPlane):
FlatCombImg.quantize()
FlatCombImg = galsim.ImageUS(FlatCombImg)
# FlatCombImg.write("%s/FlatImg_%s_%s_%s.fits" % (chip_output.subdir, FlatTag, self.chipID, i+1))
datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
date_obs = datetime_obs.strftime("%y%m%d")
time_obs = datetime_obs.strftime("%H%M%S")
# datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
# date_obs = datetime_obs.strftime("%y%m%d")
# time_obs = datetime_obs.strftime("%H%M%S")
timestamp_obs += 10 * 60
self.outputCal(
img=FlatCombImg,
......@@ -738,10 +748,9 @@ class Chip(FocalPlane):
img_rot=img_rot,
im_type='FLAT',
pointing_ID=pointing_ID,
date_obs=date_obs,
time_obs=time_obs,
output_dir=chip_output.subdir,
exptime=self.flat_exptime)
exptime=self.flat_exptime,
timestamp=timestamp_obs)
del FlatCombImg, FlatSingle, prnu_img
# flat_img.replaceNegative(replace_value=0)
......@@ -784,9 +793,9 @@ class Chip(FocalPlane):
cr_map[cr_map < 0] = 0
crmap_gsimg = galsim.Image(cr_map, dtype=np.uint16)
del cr_map
datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
date_obs = datetime_obs.strftime("%y%m%d")
time_obs = datetime_obs.strftime("%H%M%S")
# datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
# date_obs = datetime_obs.strftime("%y%m%d")
# time_obs = datetime_obs.strftime("%H%M%S")
self.outputCal(
img=crmap_gsimg,
ra_cen=ra_cen,
......@@ -794,10 +803,9 @@ class Chip(FocalPlane):
img_rot=img_rot,
im_type='CRD',
pointing_ID=pointing_ID,
date_obs=date_obs,
time_obs=time_obs,
output_dir=chip_output.subdir,
exptime=self.dark_exptime)
exptime=self.dark_exptime,
timestamp=timestamp_obs)
del crmap_gsimg
# Non-Linearity for Dark
......@@ -840,7 +848,7 @@ class Chip(FocalPlane):
readout_noise = galsim.GaussianNoise(rng=rng_readout, sigma=self.read_noise)
DarkCombImg.addNoise(readout_noise)
DarkCombImg = effects.ApplyGainNonUniform16(
DarkCombImg, self.gain_channel = effects.ApplyGainNonUniform16(
DarkCombImg, gain=self.gain,
nsecy = 2, nsecx=8,
seed=SeedGainNonuni+self.chipID,
......@@ -853,9 +861,9 @@ class Chip(FocalPlane):
DarkCombImg.quantize()
DarkCombImg = galsim.ImageUS(DarkCombImg)
# DarkCombImg.write("%s/DarkImg_%s_%s_%s.fits" % (chip_output.subdir, DarkTag, self.chipID, i+1))
datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
date_obs = datetime_obs.strftime("%y%m%d")
time_obs = datetime_obs.strftime("%H%M%S")
# datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
# date_obs = datetime_obs.strftime("%y%m%d")
# time_obs = datetime_obs.strftime("%H%M%S")
timestamp_obs += 10 * 60
self.outputCal(
img=DarkCombImg,
......@@ -864,10 +872,9 @@ class Chip(FocalPlane):
img_rot=img_rot,
im_type='DARK',
pointing_ID=pointing_ID,
date_obs=date_obs,
time_obs=time_obs,
output_dir=chip_output.subdir,
exptime=self.dark_exptime)
exptime=self.dark_exptime,
timestamp = timestamp_obs)
del DarkCombImg
# img = galsim.ImageUS(img)
......
ObservationSim/Instrument/Chip/Effects.py
View file @ adb90dd3
......@@ -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 @ adb90dd3
......@@ -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 @ adb90dd3
......@@ -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 @ adb90dd3
......@@ -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 @ adb90dd3
......@@ -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 @ adb90dd3
......@@ -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.)
chip_wcs = galsim.FitsWCS(header=h_ext)
......@@ -360,12 +370,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 @ adb90dd3
......@@ -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 @ adb90dd3
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.lib').joinpath('libstraylight.so') as dllFn:
self.slcdll = ctypes.CDLL(dllFn)
except AttributeError:
with pkg_resources.path('ObservationSim.Straylight.lib', 'libstraylight.so') 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 @ adb90dd3
from .Straylight import Straylight
from .SkybackgroundMap import *
ObservationSim/Straylight/data/Zodiacal_map1.dat 0 → 100644
View file @ adb90dd3
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 @ adb90dd3
#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 @ adb90dd3
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 @ adb90dd3
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