diff --git a/csst_mci_sim/csst_mci_sim.py b/csst_mci_sim/csst_mci_sim.py new file mode 100644 index 0000000000000000000000000000000000000000..acc18e6b5fef28fe624ada6ca37cda0a9ab4ea80 --- /dev/null +++ b/csst_mci_sim/csst_mci_sim.py @@ -0,0 +1,5813 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Mon Oct 18 09:38:35 2021 + + +@author: yan, zhaojun""" + + +""" +The CSST MCI Image Simulator +============================================= + +This file contains an image simulator for the CSST MCI. + +The approximate sequence of events in the simulator is as follows: + + #. Read in a configuration file, which defines for example, + detector characteristics (bias, dark and readout noise, gain, + plate scale and pixel scale, exposure time etc.). + #. Read in another file containing charge trap definitions (for CTI modelling). + #. Read in a file defining the cosmic rays (trail lengths and cumulative distributions). + + + #. Loop over the number of exposures to co-add and for each object in the object catalog: + + + + #. Apply a multiplicative flat-field map to emulate pixel-to-pixel non-uniformity [optional]. + + #. Add cosmic ray tracks onto the CCD with random positions but known distribution [optional]. + #. Apply detector charge bleeding in column direction [optional]. + #. Add constant dark current [optional]. + + #. Add photon (Poisson) noise [optional] + #. Add cosmetic defects from an input file [optional]. + #. Add overscan regions in the serial direction [optional]. + #. Apply the CDM03 radiation damage model [optional]. + #. Apply non-linearity model to the pixel data [optional]. + #. Add readout noise selected [optional]. + #. Convert from electrons to ADUs using a given gain factor. + #. Add a given bias level and discretise the counts (the output is going to be in 16bit unsigned integers). + #. Finally the simulated image is converted to a FITS file, a WCS is assigned + and the output is saved to the current working directory. + + # add galaxy lensing effect to self.information['lensing'], set this parameter in .config file; update@24.10.16 + +.. Warning:: The code is still work in progress and new features are being added. + The code has been tested, but nevertheless bugs may be lurking in corners, so + please report any weird or inconsistent simulations to the author. + + +Contact Information: zhaojunyan@shao.ac.cn + +------- + +""" +#from scipy.integrate import simps +import math +import sys +import os +from astropy import units as u +from astropy.io import fits +from astropy.table import Table +from tqdm import tqdm +import numpy as np +from astropy.coordinates import get_sun +from astropy.time import Time +from datetime import datetime, timedelta +import pandas as pd +import galsim +import configparser as ConfigParser +from scipy import ndimage +from astropy.coordinates import SkyCoord +from scipy import interpolate +from scipy.signal import fftconvolve +import pandas +import ctypes +import astropy.coordinates as coord +from scipy.interpolate import interp1d +import ipy.spatial as spatial +import sc +import +#from matplotlib import pyplot as plt + +###################### +try: + from support import logger as lg + from support import cosmicrays + from support import shao + from support import sed + from support import MCIinstrumentModel + from mci_so import cdm03bidir + +except: + from .support import logger as lg + from .support import cosmicrays + from .support import shao + from .support import sed + from .support import MCIinstrumentModel + from .mci_so import cdm03bidir + +#################################################### +#from joblib import Parallel, delayed + +########################### functions ######################### + + +def get_file_extension(filename): + # 使用os.path.splitext()函数分割文件名和扩展名 + file_name, extension = os.path.splitext(filename) + return extension + + +########################################################## + +""" +Charge Transfer Inefficiency +============================ + +This file contains a simple class to run a CDM03 CTI model developed by Alex Short (ESA). + +This now contains both the official CDM03 and a new version that allows different trap +parameters in parallel and serial direction. + +:requires: NumPy +:requires: CDM03 (FORTRAN code, f2py -c -m cdm03bidir cdm03bidir.f90) + +:version: 0.35 +""" + + +# CDM03bidir +class CDM03bidir(): + """ + Class to run CDM03 CTI model, class Fortran routine to perform the actual CDM03 calculations. + + :param settings: input parameters + :type settings: dict + :param data: input data to be radiated + :type data: ndarray + :param log: instance to Python logging + :type log: logging instance + """ + + def __init__(self, settings, data, log=None): + """ + Class constructor. + + :param settings: input parameters + :type settings: dict + :param data: input data to be radiated + :type data: ndarray + :param log: instance to Python logging + :type log: logging instance + """ + self.data = data + self.values = dict(quads=(0, 1, 2, 3), xsize=2048, + ysize=2066, dob=0.0, rdose=8.0e9) + self.values.update(settings) + self.log = log + self._setupLogger() + + # default CDM03 settings + self.params = dict(beta_p=0.6, beta_s=0.6, fwc=200000., vth=1.168e7, vg=6.e-11, t=20.48e-3, + sfwc=730000., svg=1.0e-10, st=5.0e-6, parallel=1., serial=0.0) + # update with inputs + self.params.update(self.values) + + # read in trap information + trapdata = np.loadtxt( + self.values['dir_path']+self.values['paralleltrapfile']) + if trapdata.ndim > 1: + self.nt_p = trapdata[:, 0] + self.sigma_p = trapdata[:, 1] + self.taur_p = trapdata[:, 2] + else: + # only one trap species + self.nt_p = [trapdata[0],] + self.sigma_p = [trapdata[1],] + self.taur_p = [trapdata[2],] + + trapdata = np.loadtxt( + self.values['dir_path']+self.values['serialtrapfile']) + if trapdata.ndim > 1: + self.nt_s = trapdata[:, 0] + self.sigma_s = trapdata[:, 1] + self.taur_s = trapdata[:, 2] + else: + # only one trap species + self.nt_s = [trapdata[0],] + self.sigma_s = [trapdata[1],] + self.taur_s = [trapdata[2],] + + # scale thibaut's values + if 'thibaut' in self.values['parallelTrapfile']: + self.nt_p /= 0.576 # thibaut's values traps / pixel + self.sigma_p *= 1.e4 # thibaut's values in m**2 + if 'thibaut' in self.values['serialTrapfile']: + self.nt_s *= 0.576 # thibaut's values traps / pixel #should be division? + self.sigma_s *= 1.e4 # thibaut's values in m**2 + + def _setupLogger(self): + """ + Set up the logger. + """ + self.logger = True + # if self.log is None: + # self.logger = False + + def applyRadiationDamage(self, data, iquadrant=0): + """ + Apply radian damage based on FORTRAN CDM03 model. The method assumes that + input data covers only a single quadrant defined by the iquadrant integer. + + :param data: imaging data to which the CDM03 model will be applied to. + :type data: ndarray + + :param iquandrant: number of the quadrant to process + :type iquandrant: int + + cdm03 - Function signature:: + + sout = cdm03(sinp,iflip,jflip,dob,rdose,in_nt,in_sigma,in_tr,[xdim,ydim,zdim]) + Required arguments: + sinp : input rank-2 array('d') with bounds (xdim,ydim) + iflip : input int + jflip : input int + dob : input float + rdose : input float + in_nt : input rank-1 array('d') with bounds (zdim) + in_sigma : input rank-1 array('d') with bounds (zdim) + in_tr : input rank-1 array('d') with bounds (zdim) + Optional arguments: + xdim := shape(sinp,0) input int + ydim := shape(sinp,1) input int + zdim := len(in_nt) input int + Return objects: + sout : rank-2 array('d') with bounds (xdim,ydim) + + .. Note:: Because Python/NumPy arrays are different row/column based, one needs + to be extra careful here. NumPy.asfortranarray will be called to get + an array laid out in Fortran order in memory. Before returning the + array will be laid out in memory in C-style (row-major order). + + :return: image that has been run through the CDM03 model + :rtype: ndarray + """"" + # return data + + iflip = iquadrant / 2 + jflip = iquadrant % 2 + + params = [self.params['beta_p'], self.params['beta_s'], self.params['fwc'], self.params['vth'], + self.params['vg'], self.params['t'], self.params['sfwc'], self.params['svg'], self.params['st'], + self.params['parallel'], self.params['serial']] + + if self.logger: + self.log.info('nt_p=' + str(self.nt_p)) + self.log.info('nt_s=' + str(self.nt_s)) + self.log.info('sigma_p= ' + str(self.sigma_p)) + self.log.info('sigma_s= ' + str(self.sigma_s)) + self.log.info('taur_p= ' + str(self.taur_p)) + self.log.info('taur_s= ' + str(self.taur_s)) + self.log.info('dob=%f' % self.values['dob']) + self.log.info('rdose=%e' % self.values['rdose']) + self.log.info('xsize=%i' % data.shape[1]) + self.log.info('ysize=%i' % data.shape[0]) + self.log.info('quadrant=%i' % iquadrant) + self.log.info('iflip=%i' % iflip) + self.log.info('jflip=%i' % jflip) + +################################################################################# + + CTIed = cdm03bidir.cdm03(np.asfortranarray(data), + jflip, iflip, + self.values['dob'], self.values['rdose'], + self.nt_p, self.sigma_p, self.taur_p, + self.nt_s, self.sigma_s, self.taur_s, + params, + [data.shape[0], data.shape[1], len(self.nt_p), len(self.nt_s), len(self.params)]) + + return np.asanyarray(CTIed) + +################################################################################################################# +################################################################################################################# + + +def transRaDec2D(ra, dec): + # radec转为竞天程序里的ob, 赤道坐标系下的笛卡尔三维坐标xyz. + 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 ill2flux(E, path): + + # use template from sky_bkg (background_spec_hst.dat) + filename = path+'MCI_inputData/refs/background_spec_hst.dat' + cat_spec = pd.read_csv(filename, sep='\s+', header=None, comment='#') + wave0 = cat_spec[0].values # A + spec0 = cat_spec[2].values # erg/s/cm^2/A/arcsec^2 + + # convert erg/s/cm^2/A/arcsec^2 to erg/s/cm^2/A/sr + flux1 = spec0 / (1/4.25452e10) + # convert erg/s/cm^2/A/sr to W/m^2/sr/um + flux2 = flux1 * 10 + + # 对接收面积积分，输出单位 W/m^2/nm + D = 2 # meter + f = 28 # meter, 焦距，转换关系来源于王维notes. + flux3 = flux2 * np.pi * D**2 / 4 / f**2 / 10**3 + + f = interp1d(wave0, flux3) + wave_range = np.arange(3800, 7800) + flux3_mean = f(wave_range) + delta_lamba = 0.1 # nm + E0 = np.sum(flux3_mean * delta_lamba) + + factor = E / E0 + spec_scaled = factor * spec0 + + return wave0, spec_scaled + +############################################################## + + +def earth_angle(time_jd, x_sat, y_sat, z_sat, ra_obj, dec_obj): + + ra_sat = np.arctan2(y_sat, x_sat) / np.pi * 180 + dec_sat = np.arctan2(z_sat, np.sqrt(x_sat**2+y_sat**2)) / np.pi * 180 + radec_sat = SkyCoord(ra=ra_sat*u.degree, + dec=dec_sat*u.degree, frame='gcrs') + lb_sat = radec_sat.transform_to('geocentrictrueecliptic') + + # get the obj location + radec_obj = SkyCoord(ra=ra_obj*u.degree, + dec=dec_obj*u.degree, frame='gcrs') + lb_obj = radec_obj.transform_to('geocentrictrueecliptic') + + # calculate the angle between sub-satellite point and the earth side + earth_radius = 6371 # km + sat_height = np.sqrt(x_sat**2 + y_sat**2 + z_sat**2) + angle_a = np.arcsin(earth_radius/sat_height) / np.pi * 180 + + # calculate the angle between satellite position and the target position + angle_b = lb_sat.separation(lb_obj) + + # calculat the earth angle + angle = 180 - angle_a - angle_b.degree + + return angle +####################################### + + +def CCDnonLinearityModel(data, beta=6e-7): + """ + + The non-linearity is modelled based on the results presented. + :param data: data to which the non-linearity model is being applied to + :type data: ndarray + + :return: input data after conversion with the non-linearity model + :rtype: float or ndarray + """ + out = data-beta*data**2 + + return out +# 33333 + + +class StrayLight(object): + + def __init__(self, path, jtime=2460843., sat=np.array([0, 0, 0]), radec=np.array([0, 0])): + self.path = path + self.jtime = jtime + self.sat = sat + self.equator = coord.SkyCoord( + radec[0]*u.degree, radec[1]*u.degree, frame='icrs') + self.ecliptic = self.equator.transform_to('barycentrictrueecliptic') + self.pointing = transRaDec2D(radec[0], radec[1]) + self.slcdll = ctypes.CDLL( + self.path+'MCI_inputData/refs/libstraylight.so') # 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] + self.deFn = self.path+"MCI_inputData/refs/DE405" + self.PSTFn = self.path+"MCI_inputData/refs/PST" + self.RFn = self.path+"MCI_inputData/refs/R" + self.ZolFn = self.path+"MCI_inputData/refs/Zodiacal" + self.brightStarTabFn = self.path+"MCI_inputData/refs/BrightGaia_with_csst_mag" + + self.slcdll.Init(str.encode(self.deFn), str.encode( + self.PSTFn), str.encode(self.RFn), str.encode(self.ZolFn)) + + def caculateStarLightFilter(self, filter='i'): + + filterIndex = {'nuv': 0, 'u': 1, 'g': 2, + 'r': 3, 'i': 4, 'z': 5, 'y': 6} + + 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) + + star_e1 = (ctypes.c_double*7)() + self.slcdll.PointSource(self.jtime, sat, ob, py1, + star_e1, str.encode(self.brightStarTabFn)) + star_e2 = (ctypes.c_double*7)() + self.slcdll.PointSource(self.jtime, sat, ob, py2, + star_e2, str.encode(self.brightStarTabFn)) + + band_star_e1 = star_e1[:][filterIndex[filter]] + band_star_e2 = star_e2[:][filterIndex[filter]] + + return max(band_star_e1, band_star_e2) + + def caculateEarthShineFilter(self, filter='i'): + + filterIndex = {'nuv': 0, 'u': 1, 'g': 2, + 'r': 3, 'i': 4, 'z': 5, 'y': 6} + 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) # e[7]代表7个波段的照度 + + 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]] + + return max(band_earth_e1, band_earth_e2) + +############################################################################### +############################################################################### +############################################################################### + + +def make_c_coor(fov, step): + """ + Draw the mesh grids for a fov * fov box with given step . + """ + + nc = int(fov/step) + + bs = fov + ds = bs / nc + xx01 = np.linspace(-bs / 2.0, bs / 2.0 - ds, nc) + 0.5 * ds + xx02 = np.linspace(-bs / 2.0, bs / 2.0 - ds, nc) + 0.5 * ds + xg2, xg1 = np.meshgrid(xx01, xx02) + return xg1, xg2 + +############################################################################## + +############################################################################## + +# datetime to mjd + + +def time2mjd(dateT): + t0 = datetime(1858, 11, 17, 0, 0, 0, 0) # + mjd = (dateT-t0).days + mjd_s = dateT.hour*3600.0+dateT.minute*60.0 + \ + dateT.second+dateT.microsecond/1000000.0 + return mjd+mjd_s/86400.0 + + +def time2jd(dateT): + t0 = datetime(1858, 11, 17, 0, 0, 0, 0) # + mjd = (dateT-t0).days + mjd_s = dateT.hour*3600.0+dateT.minute*60.0 + \ + dateT.second+dateT.microsecond/1000000.0 + return mjd+mjd_s/86400.0++2400000.5 + +# mjd to datetime + + +def mjd2time(mjd): + t0 = datetime(1858, 11, 17, 0, 0, 0, 0) + return t0+timedelta(days=mjd) + + +def jd2time(jd): + mjd = jd2mjd(jd) + return mjd2time(mjd) + +# mjd to jd +# def mjd2jd(mjd): +# return mjd+2400000.5 + + +def jd2mjd(jd): + return jd-2400000.5 + + +def dt2hmd(dt): + ## dt is datetime + hour = dt.hour + minute = dt.minute + second = dt.second + if hour < 10: + str_h = '0'+str(hour) + else: + str_h = str(hour) + + if minute < 10: + str_m = '0'+str(minute) + else: + str_m = str(minute) + + if second < 10: + str_d = '0'+str(second+dt.microsecond*1e-6) + else: + str_d = str(second+dt.microsecond*1e-6) + + return str_h+':'+str_m+':'+str_d + +############################################################################### + + +def deg2HMS(ra0, dec0): + '''convert deg to ra's HMS and dec's DMS''' + c = SkyCoord(ra=ra0*u.degree, dec=dec0*u.degree) + ss = c.to_string('hmsdms') + + # print(ss) + + for k in range(5, len(ss)-1): + + if ss[k] == 's': + s_idx = k + + if ss[k] == 'd': + d_idx = k + + if ss[k] == 'm': + m_idx = k + + if c.ra.hms.s < 10 and c.ra.hms.s > 0.1: + temp = str(c.ra.hms.s) + ra_ss = '0'+temp[:3] + if c.ra.hms.s > 10: + temp = str(c.ra.hms.s) + ra_ss = temp[:4] + + if c.ra.hms.s < 0.1: + temp = str(c.ra.hms.s) + ra_ss = '00.0' + + dms_d = ss[s_idx+2:d_idx] + dms_m = ss[d_idx+1:m_idx] + dms_s = ss[m_idx+1:m_idx+3] + + hhmmss = ss[0:2]+ss[3:5]+ra_ss + return hhmmss+dms_d+dms_m+dms_s + +############################################################################### + +############################################################################### +############################################################################# +# distortion function + + +def distortField(ra, dec, ch): + # % ra ,dec are the idea position in arcsec , and the center position is ra=0, dec=0 + '''MCI geometry distortion effect in channel G R and I''' + distortA = dict() + + distortB = dict() + + distortA['g'] = np.array([0.000586224851462036, 0.999778507355332, -0.000377391397200834, -4.40403573019393e-06, 0.000147444528630966, -1.93281825050268e-06, 5.73520238714447e-07, -1.05194725549731e-08, 7.55124671251098e-07, -3.85623216175251e-09, - + 1.36254798567168e-11, 1.36983654024573e-09, -4.53104347730230e-11, 1.50641840169747e-09, -1.05226890727366e-11, 1.06471556810228e-12, -1.42299485385165e-16, 5.90008722878028e-12, 2.36749977009538e-17, 4.11061448730916e-12, 3.39287277469805e-17]) + distortB['g'] = np.array([1.00113878750680, -0.000495107770623867, 0.999162436209357, 9.67138672907941e-05, -3.88808001621361e-06, 0.000267940195644359, -3.03504629416955e-09, 7.59508802931395e-07, -1.13952684052500e-08, 9.58672893217047e-07, + 2.96397490595616e-10, -3.01506961685930e-11, 2.22570315950441e-09, -4.26077028191289e-11, 2.07336026666512e-09, 1.71450079597168e-16, 2.94455108664049e-12, 8.18246797239659e-17, 8.32235684990184e-12, 1.05203957047210e-17, 5.58541337682320e-12]) + + distortA['r'] = np.array([0.000530712822898294, 0.999814397089242, -0.000366783811357609, -1.08650906916023e-05, 0.000146500108063480, -4.48994291741769e-06, 5.66992328511032e-07, -2.10826363791224e-08, 7.51050898843171e-07, -7.74459106063614e-09, - + 5.78712436084704e-11, 1.36389366137393e-09, -9.50057131576629e-11, 1.49666815592076e-09, -2.16074939825761e-11, 2.07124539578673e-12, -1.32787329448528e-13, 5.78542850850562e-12, -5.82328830750271e-14, 4.04881815088026e-12, 2.46095156355282e-15]) + distortB['r'] = np.array([1.00110972320988, -0.000483253233767592, 0.999181377618578, 9.60316928622784e-05, -9.06574382424422e-06, 0.000266147076486786, -6.64513480754565e-09, 7.53794491639207e-07, -2.25340150955077e-08, 9.53032549996219e-07, + 2.93844965469408e-10, -6.80521155195455e-11, 2.21272371816576e-09, -9.13097728847483e-11, 2.06225316601308e-09, 9.41096324034730e-14, 3.00253083795091e-12, -1.40935245750903e-13, 8.27122551593984e-12, -2.41500808188601e-13, 5.52074803508046e-12]) + + distortA['i'] = np.array([0.000532645905256854, 0.999813271238129, -0.000366606839338804, -1.06782246147986e-05, 0.000146529811926289, -4.54488847969004e-06, 5.68028930846942e-07, -2.11055358580630e-08, 7.51187628288423e-07, -7.93885390157909e-09, - + 5.63104333653064e-11, 1.36559362569725e-09, -9.64353839395137e-11, 1.50231201562648e-09, -2.04159956020294e-11, 1.91408535007684e-12, -7.46369323634455e-14, 5.86071470639700e-12, -1.65328163262914e-13, 4.07648238374705e-12, 3.40880674871652e-14]) + distortB['i'] = np.array([1.00111276400037, -0.000484310769918440, 0.999180286636823, 9.61240720951134e-05, -8.76526511019577e-06, 0.000266192433565878, -6.59462433108999e-09, 7.54391611102465e-07, -2.30957980169170e-08, 9.53612393886641e-07, + 2.95558631849358e-10, -7.08307092409029e-11, 2.21952307845185e-09, -9.03816603147003e-11, 2.06450141393973e-09, -3.77836686311826e-14, 3.13358046393416e-12, -5.20417845240954e-14, 8.24487152802918e-12, -1.17846469609206e-13, 5.35830020655191e-12]) + + x = ra/0.05*10/1000 # convert ra in arcsec to mm + + y = dec/0.05*10/1000 # convert ra in arcsec to mm + + dd = np.array([1, x, y, x*x, x*y, y*y, x**3, x**2*y, x*y**2, y**3, x**4, x**3*y, + x**2*y**2, x*y**3, y**4, x**5, x**4*y, x**3*y**2, x**2*y**3, x*y**4, y**5]) + + xc = np.dot(distortA[ch], dd) + yc = np.dot(distortB[ch], dd) + + distortRa = xc*1000/10*0.05-0.00265 + + distortDec = yc*1000/10*0.05-5.0055 + + return distortRa, distortDec + +##################################################################################### +# def world_to_pixel(sra, +# sdec, +# rotsky, +# tra, +# tdec, +# x_center_pixel, +# y_center_pixel, +# pixelsize=0.05): +# """_summary_ + +# Parameters +# ---------- +# sra : np.array +# star ra such as np.array([22,23,24]),unit:deg +# sdec : np.array +# stat dec such as np.array([10,20,30]),unit:deg +# rotsky : float +# rotation angel of the telescope,unit:deg +# tra : float +# telescope pointing ra,such as 222.2 deg +# rdec : float +# telescope pointing dec,such as 33.3 deg +# x_center_pixel : float +# x refer point of MCI,usually the center of the CCD,which start from (0,0) +# y_center_pixel : float +# y refer point of MCI,usually the center of the CCD,which start from(0,0) +# pixelsize : float +# pixelsize for MCI ccd, default :0.05 arcsec / pixel + +# Returns +# ------- +# pixel_position:list with two np.array +# such as [array([124.16937605, 99. , 99. ]), +# array([ 97.52378661, 99. , 100.50014483])] +# """ +# theta_r = rotsky / 180 * np.pi +# w = WCS(naxis=2) +# w.wcs.crpix = [x_center_pixel + 1, y_center_pixel + 1] # pixel from (1, 1) +# w.wcs.cd = np.array([[ +# -np.cos(-theta_r) * pixelsize / 3600, +# np.sin(-theta_r) * pixelsize / 3600 +# ], +# [ +# np.sin(-theta_r) * pixelsize / 3600, +# np.cos(-theta_r) * pixelsize / 3600 +# ]]) +# w.wcs.crval = [tra, tdec] +# w.wcs.ctype = ["RA---TAN", "DEC--TAN"] +# pixel_position = w.all_world2pix(sra, sdec, 0) +# return pixel_position + +############################################################################### + + +def cal_pos(center_ra, center_dec, rotTel, rotSky, sRa, sDec): + # ''' center_ra: telescope pointing in ra direction, unit: degree + # center_dec: telescope pointing in dec direction, unit: degree + # rotTel: telescope totation in degree + # rotSky: telescope rotation in degree + # sRa: source ra in arcsec + # sDec: source dec in arcsec + + boresight = galsim.CelestialCoord( + ra=-center_ra*galsim.degrees, dec=center_dec*galsim.degrees) + + q = rotTel - rotSky + cq, sq = np.cos(q), np.sin(q) + jac = [cq, -sq, sq, cq] + affine = galsim.AffineTransform(*jac) + wcs = galsim.TanWCS(affine, boresight, units=galsim.arcsec) + objCoord = galsim.CelestialCoord( + ra=-sRa*galsim.arcsec, dec=sDec*galsim.arcsec) + field_pos = wcs.toImage(objCoord) + return -field_pos.x, field_pos.y # return the field position + +##################################################################################### + +# def krebin(a, sample): +# """ Fast Rebinning with flux conservation +# New shape must be an integer divisor of the current shape. +# This algorithm is much faster than rebin_array +# Parameters +# ---------- +# a : array_like +# input array +# sample : rebinned sample 2 or 3 or 4 or other int value +# """ +# # Klaus P's fastrebin from web +# size=a.shape +# shape=np.zeros(2,dtype=int) +# shape[0]=int(size[0]/sample) +# shape[1]=int(size[1]/sample) + +# sh = shape[0], a.shape[0] // shape[0], shape[1], a.shape[1] // shape[1] + +# return a.reshape(sh).sum(-1).sum(1) + +##################################################################### +######################################################### + + +def centroid(data): + ''' + calculate the centroid of the input two-dimentional image data + + Parameters + ---------- + data : input image. + + Returns + ------- + cx: the centroid column number, in horizontal direction definet in python image show + cy: the centroid row number , in vertical direction + + ''' + ### + from scipy import ndimage + cy, cx = ndimage.center_of_mass(data) + return cx, cy + +############################################################################### + +############################################################################### + + +''' +PSF interpolation for MCI_sim +''' + + +###find neighbors-KDtree ### + +def findNeighbors(tx, ty, px, py, dn=5): + """ + find nearest neighbors by 2D-KDTree + + Parameters: + tx, ty (float, float): a given position + px, py (numpy.array, numpy.array): position data for tree + dr (float-optional): distance in arcsec + dn (int-optional): nearest-N + OnlyDistance (bool-optional): only use distance to find neighbors. Default: True + + Returns: + indexq (numpy.array): index + """ + datax = px + datay = py + tree = spatial.KDTree(list(zip(datax.ravel(), datay.ravel()))) + indexq = [] + dist, indexq = tree.query([tx, ty], dn) + return indexq + +############################################################################### +###PSF-IDW### + + +def psfMaker_IDW(px, py, PSFMat, cen_col, cen_row, dn=5, IDWindex=3, OnlyNeighbors=True): + """ + psf interpolation by IDW + + Parameters: + px, py (float, float): position of the target + PSFMat (numpy.array): PSF matrix array at given position and wavelength + cen_col, cen_row (numpy.array, numpy.array): potions of the psf centers in arcsec + IDWindex (int-optional): the power index of IDW + OnlyNeighbors (bool-optional): only neighbors are used for psf interpolation + + Returns: + psfMaker (numpy.array) + """ + + minimum_psf_weight = 1e-8 + ref_col = px # target position in x + ref_row = py # target position in y + + ngy, ngx = PSFMat[:, :, 0, 0].shape # PSF data size + + #################################### + ####### my code ###### + psfWeight = np.zeros([dn]) + + if OnlyNeighbors == True: + + neigh = findNeighbors(px, py, cen_col, cen_row, dn) + + ####################################################################### + for ipsf in range(len(neigh)): + + idx = neigh[ipsf] + dist = np.sqrt((ref_col - cen_col[idx]) + ** 2 + (ref_row - cen_row[idx])**2) + + if IDWindex == 1: + psfWeight[ipsf] = dist + if IDWindex == 2: + psfWeight[ipsf] = dist**2 + if IDWindex == 3: + psfWeight[ipsf] = dist**3 + if IDWindex == 4: + psfWeight[ipsf] = dist**4 + if IDWindex == 5: + psfWeight[ipsf] = dist**5 + + psfWeight[ipsf] = max(psfWeight[ipsf], minimum_psf_weight) + psfWeight[ipsf] = 1./psfWeight[ipsf] + + # 3 + + psfWeight /= np.sum(psfWeight) + + psfMaker = np.zeros([ngy, ngx, 7], dtype=np.float32) + + ######################################################### + for waven in range(7): + + for ipsf in range(len(neigh)): + + idx = neigh[ipsf] + + iPSFMat = PSFMat[:, :, waven, idx].copy() + + ipsfWeight = psfWeight[ipsf] + + psfMaker[:, :, waven] += iPSFMat * ipsfWeight + + psfMaker[:, :, waven] /= np.nansum(psfMaker[:, :, waven]) + + return psfMaker + +############################################################################### +############################################################################### +############################################################################### + + +class MCIsimulator(): + """ + CSST MCI Simulator + + The image that is being build is in:: + + self.image + self.image_g g channel + self.image_r r channel + self.image_i i channel + + + :param opts: OptionParser instance + :type opts: OptionParser instance + """ + + def __init__(self, configfile): + """ + Class Constructor. + + :param opts: OptionParser instance + :type opts: OptionParser instance + """ + self.configfile = configfile + self.section = 'TEST' # simulation section; + self.information = MCIinstrumentModel.MCIinformation() + + # update settings with defaults + self.information.update(dict(quadrant=int(0), + ccdx=int(0), + ccdy=int(0), + ccdxgap=1.643, + ccdygap=8.116, + xsize=2000, + ysize=2000, + fullwellcapacity=90000, + pixel_size=0.05, + dark=0.001, + exptime=300.0, + readouttime=4., + rdose=8.0e9, + ghostCutoff=22.0, + ghostratio=5.e-5, + coveringfraction=0.05, # CR: + # cosmicraylengths =self.information['dir_path']+'MCI_inputData/data/cdf_cr_length.dat', + # cosmicraydistance=self.information['dir_path']+'MCI_inputData/data/cdf_cr_total.dat', + # parallelTrapfile =self.information['dir_path']+'MCI_inputData/data/cdm_euclid_parallel.dat', + # serialTrapfile =self.information['dir_path']+'MCI_inputData/data/cdm_euclid_serial.dat', + mode='same')) + +#################################################### + +############################################################################### + def readConfigs(self, simnumber, source): + """ + Reads the config file information using configParser and sets up a logger. + """ + self.config = ConfigParser.RawConfigParser() + + if simnumber == 1: + print('beging config file name : ') + print(self.configfile) + + # self.config.readfp(open(self.configfile)) + self.config.read_file(open(self.configfile)) + + +########################################################################### + + + def processConfigs(self): + """ + Processes configuration information and save the information to a dictionary self.information. + For explanation of each field, see /data/test.config. Note that if an input field does not exist, + then the values are taken from the default instrument model as described in + support.MCIinstrumentModel.VISinformation(). Any of the defaults can be overwritten by providing + a config file with a correct field name. + """ + # parse options and update the information dictionary + options = self.config.options(self.section) + settings = {} + for option in options: + try: + settings[option] = self.config.getint(self.section, option) + except ValueError: + try: + settings[option] = self.config.getfloat( + self.section, option) + except ValueError: + settings[option] = self.config.get(self.section, option) + + self.information.update(settings) + + # ghost ratio can be in engineering format, so getfloat does not capture it... + self.information['ghostratio'] = float( + self.config.get(self.section, 'ghostRatio')) + + ################################################################################################### + ################################################################################################## + + # booleans to control the flow + + self.cosmicRays = self.config.getboolean(self.section, 'cosmicRays') + self.darknoise = self.config.getboolean(self.section, 'darknoise') + self.cosmetics = self.config.getboolean(self.section, 'cosmetics') + self.radiationDamage = self.config.getboolean( + self.section, 'radiationDamage') + self.bleeding = self.config.getboolean(self.section, 'bleeding') + self.overscans = self.config.getboolean(self.section, 'overscans') + self.nonlinearity = self.config.getboolean( + self.section, 'nonlinearity') + self.readoutNoise = self.config.getboolean( + self.section, 'readoutNoise') + self.skyback = self.config.getboolean(self.section, 'skyback') + self.TianceEffect = self.config.getboolean( + self.section, 'TianceEffect') + self.intscale = self.config.getboolean(self.section, 'intscale') + self.ghosts = self.config.getboolean(self.section, 'ghosts') + self.shutterEffect = self.config.getboolean( + self.section, 'shutterEffect') + self.flatfieldM = self.config.getboolean(self.section, 'flatfieldm') + self.PRNUeffect = self.config.getboolean(self.section, 'PRNUeffect') + self.appFatt = self.config.getboolean(self.section, 'appFatt') + self.sky_shift_rot = self.config.getboolean( + self.section, 'sky_shift_rot') + self.distortion = self.config.getboolean(self.section, 'distortion') + self.sim_star = self.config.getboolean(self.section, 'sim_star') + self.sim_galaxy = self.config.getboolean(self.section, 'sim_galaxy') + + self.save_starpsf = self.config.getboolean( + self.section, 'save_starpsf') + self.save_cosmicrays = self.config.getboolean( + self.section, 'save_cosmicrays') + self.lensing = self.config.getboolean(self.section, 'lensing') + + ###############################################3################################# + self.booleans = dict(cosmicRays=self.cosmicRays, + darknoise=self.darknoise, + cosmetics=self.cosmetics, + radiationDamage=self.radiationDamage, + bleeding=self.bleeding, + overscans=self.overscans, + nonlinearity=self.nonlinearity, + readoutNoise=self.readoutNoise, + skyback=self.skyback, + TianceEffect=self.TianceEffect, + intscale=self.intscale, + ghosts=self.ghosts, + shutterEffect=self.shutterEffect, + flatfieldM=self.flatfieldM, + PRNUeffect=self.PRNUeffect, + appFatt=self.appFatt, + sky_shift_rot=self.sky_shift_rot, + distortion=self.distortion, + sim_star=self.sim_star, + sim_galaxy=self.sim_galaxy, + save_starpsf=self.save_starpsf, + save_cosmicrays=self.save_cosmicrays, + lensing=self.lensing) + ##################################################################### + + return + + ############################################################################## +############################################################################### + + def _createEmpty(self): + """ + Creates and empty array of a given x and y size full of zeros. + add g r i channel images; + + Creates lensing parameters; + + """ + self.image_g = np.zeros( + (self.information['ysize'], self.information['xsize']), dtype=float) + self.image_r = np.zeros( + (self.information['ysize'], self.information['xsize']), dtype=float) + self.image_i = np.zeros( + (self.information['ysize'], self.information['xsize']), dtype=float) + return + +############################################################################### +############################################################################## + def _loadGhostModel(self): + """ + Reads in a ghost model from a FITS file and stores the data to self.ghostModel. + + Currently assumes that the ghost model has already been properly scaled and that the pixel + scale of the input data corresponds to the nominal VIS pixel scale. Futhermore, assumes that the + distance to the ghost from y=0 is appropriate (given current knowledge, about 750 VIS pixels). + """ + + self.ghostOffsetX = self.information['ghostoffsetx'] + self.ghostOffsetY = self.information['ghostoffsety'] + + self.ghostMax = self.information['ghostratio'] + self.log.info('Maximum in the ghost model %e' % self.ghostMax) + return + ############################################### + + def readCosmicRayInformation(self): + """ + Reads in the cosmic ray track information from two input files. + + Stores the information to a dictionary called cr. + """ + self.log.info('Reading in cosmic ray information from %s and %s' % (self.information['cosmicraylengths'], + self.information['cosmicraydistance'])) + + crLengths = np.loadtxt( + self.information['dir_path']+self.information['cosmicraylengths']) + crDists = np.loadtxt( + self.information['dir_path']+self.information['cosmicraydistance']) + + self.cr = dict(cr_u=crLengths[:, 0], cr_cdf=crLengths[:, 1], cr_cdfn=np.shape(crLengths)[0], + cr_v=crDists[:, 0], cr_cde=crDists[:, 1], cr_cden=np.shape(crDists)[0]) + + return + +############################################################################### + +############################################################################### + def configure(self, simnumber, sourcein, dir_path, result_path): + """ + Configures the simulator with input information and creates and empty array to which the final image will + be build on. + """ + self.readConfigs(simnumber, sourcein) + + self.processConfigs() + + self._createEmpty() + + self.information['dir_path'] = dir_path + + self.information['result_path'] = result_path + + self.source = sourcein + + ##print('print information:', self.information) + + ########################################################################### + + now = datetime.utcnow() + # data_time=now.strftime("%Y-%m-%d-%H-%M-%S") + result_day = now.strftime("%Y-%m-%d") + + # if self.information['dir_path'] =='/nfsdata/share/simulation-unittest/mci_sim/': + # self.result_path = self.information['dir_path'] +'mci_sim_result/'+result_day + # else: + + # home_path = os.environ['HOME'] + + # if home_path == '/home/yan': + + # self.result_path = '../MCI_simResult/'+self.source+"_"+result_day + # else: + # self.result_path = '/data/mcisimdata/'+result_day + + print(self.information['result_path']) + + self.result_path = self.information['result_path'] + \ + self.source+"_"+result_day + + if os.path.isdir(self.result_path) == False: + os.mkdir(self.result_path) + os.mkdir(self.result_path+'/cali_Data') + os.mkdir(self.result_path+'/log_Data') + os.mkdir(self.result_path+'/ori_Cali') + os.mkdir(self.result_path+'/ori_Sky') + os.mkdir(self.result_path+'/sky_Data') + os.mkdir(self.result_path+'/PSF_Data') + + now = datetime.utcnow() + + data_time = now.strftime("%Y-%m-%d-%H-%M-%S") + + print('simnumber', simnumber) + + self.log = lg.setUpLogger(self.result_path+'/log_Data/MCIsim_' + + self.source+'_'+data_time+'_Num_'+str(simnumber)+'.log') + + self.log.info('-------STARTING A NEW SIMULATION------------') + + for key, value in self.information.items(): + self.log.info('%s = %s' % (key, value)) + + self.log.info('The exposure order is %i ' % simnumber) + + # load instrument model, these values are also stored in the FITS header + + self.information['G_filters'] = ["F275W", "F280N", "NUV", + "WU", "CBU", "F343N", "u", "F373N", "F395N", "F336W"] + self.information['R_filters'] = ["F487N", "F502N", "CBV", "r", + "F656N", "F658N", "F467M", "F555W", "F606W", "F673N"] + self.information['I_filters'] = ["z", "y", "F815N", + "CBI", "F925N", "F960M", "F968N", "F845M", "F850LP", "F814W"] + + # load telescope efficiency data + + self.tel_eff = np.load( + self.information['dir_path']+'MCI_inputData/tel_eff/tel_eff.npy', allow_pickle=True).item() + + # load MCI filter data + self.filterP = np.load( + self.information['dir_path']+'MCI_inputData/MCI_filters/mci_filterPWTC.npy', allow_pickle=True).item() + + return + +################################################################################ + def load_filter_PSF(self): + # load filter PSF in three channels ,G , R and I + ### + # load filter name in three channels from information + self.filter_g = self.information['filter_g'] + self.filter_r = self.information['filter_r'] + self.filter_i = self.information['filter_i'] + + self.filter_psf = dict() + + filtername = self.filter_g + self.filter_psf['g'] = np.load( + self.information['dir_path']+'MCI_inputData/PSF/'+filtername+'_PSF.npy', allow_pickle=True).item() + + filtername = self.filter_r + self.filter_psf['r'] = np.load( + self.information['dir_path']+'MCI_inputData/PSF/'+filtername+'_PSF.npy', allow_pickle=True).item() + + filtername = self.filter_i + self.filter_psf['i'] = np.load( + self.information['dir_path']+'MCI_inputData/PSF/'+filtername+'_PSF.npy', allow_pickle=True).item() + + return + +############################################################################################### + + def get_PSF(self, px, py, ch): + """ + Get the PSF at a given image position + + Parameters: + px,py : target position + filternmae : the selected filter name + ch : the MCI channle , g or r or i + Returns: + PSF: PSF array at given 7 wavelength and at target position + """ + PSFMat = self.filter_psf[ch]['psf_mat'] + cen_col = self.filter_psf[ch]['psf_field_X'] + cen_row = self.filter_psf[ch]['psf_field_Y'] + imPSF = psfMaker_IDW(px, py, PSFMat, cen_col, cen_row, IDWindex=2) + + return imPSF + +############################################################################### + +############################################################################### + + def cal_SED_photons(self, flux): + # flux_arr, input flux array with unit of 1e-17 erg/s/A/cm^2 + + # set the wavelenth for MCI from 250nm to 1100nm + wave = np.linspace(2500, 10000, 8501) + wave = wave/10 + + # calcutle photons from original spec cube data, + # data: input data, two-dimentional, unit : 1e-17 erg/s/A/cm^2 + + planckh = 6.62620*10**-27 # % erg s; + cc = 2.99792458*10**17 # in nm/s + telarea = 3.1415926*100*100 # in cm^2，望远镜聚光面积�? + # convert original unit to unit of erg/s/A/cm^2 + fluxlam = 1e-17*(flux) + # wave in nm ;; + # single photon energy in erg/photon; cc与lambda单位需要一致； + ephoton = planckh*cc/wave + + Nphoton = fluxlam/ephoton*telarea * \ + self.information['exptime'] # in unit of photons + + # load telescope efficiency data + teleff = dict() + teleff['g'] = self.tel_eff['G'] # fisrt colunmm is wavelength in nm + teleff['r'] = self.tel_eff['R'] # second colunmm is efficiency + teleff['i'] = self.tel_eff['I'] + ##### + + ft = dict() + ft['g'] = self.filter_g + ft['r'] = self.filter_r + ft['i'] = self.filter_i + + # load filter name in three channels + ft_wave = dict() + ft_wave['g'] = self.filterP[ft['g']]['wave_nm'] + ft_wave['r'] = self.filterP[ft['r']]['wave_nm'] + ft_wave['i'] = self.filterP[ft['i']]['wave_nm'] + + ft_eff = dict() + ft_eff['g'] = self.filterP[ft['g']]['throughout'] + ft_eff['r'] = self.filterP[ft['r']]['throughout'] + ft_eff['i'] = self.filterP[ft['i']]['throughout'] + + chlist = ['g', 'r', 'i'] + + Sumphoton = dict() + PSF_eff_weight = dict() + + for k in range(3): + ch = chlist[k] + tel_wavearr = teleff[ch][:, 0] + tel_effarr = teleff[ch][:, 1] + tel_effnew = np.interp(wave, tel_wavearr, tel_effarr) + + ft_wavearr = ft_wave[ch] + ft_effarr = ft_eff[ch] + ft_effnew = np.interp(wave, ft_wavearr, ft_effarr) + + # SED flux*filter_flux*tel_eff + eff_arr = Nphoton*ft_effnew*tel_effnew + + Sumphoton[ch] = np.sum(eff_arr) + + # seven wavelength for the seleced filter + iwave = self.filter_psf[ch]['psf_iwave'][:] + + psf_coeff = np.interp(iwave, wave, eff_arr) + + if psf_coeff.sum() == 0: + psf_coeff = 1.0/7*np.ones(7) + else: + psf_coeff = psf_coeff/psf_coeff.sum() + + PSF_eff_weight[ch] = psf_coeff + + return Sumphoton, PSF_eff_weight + + ################################################## +############################################################################### + + def cal_sky_noise(self): + ####### add earthshine ####### + ####### add earthshine ####### + # self.earthshine_wave # A + # self.earthshine_flux # erg/s/cm^2/A/arcsec^2 + + # set the wavelenth for MCI from 250nm to 1100nm + wave = np.linspace(2500, 11000, 8501) + + skynoise_flux = self.earthshine_flux+self.zodiacal_flux # erg/s/cm^2/A/arcsec^2 + + wave = wave/10 + + # data: input data, two-dimentional, unit : 1e-17 erg/s/A/cm^2 + planckh = 6.62620*10**-27 # % erg s; + cc = 2.99792458*10**17 # in nm/s + telarea = 3.1415926*100*100 # in cm^2， + ephoton = planckh*cc/wave # in erg/photon, + # each pixel will have the photons in unit of photons/s/A + Nphoton_flux = skynoise_flux/ephoton*telarea * \ + 0.05*0.05*self.information['exptime'] + + # load telescope efficiency data + teleff = dict() + teleff['g'] = self.tel_eff['G'] # fisrt colunmm is wavelength in nm + teleff['r'] = self.tel_eff['R'] # second colunmm is efficiency + teleff['i'] = self.tel_eff['I'] + ##### + + ft = dict() + ft['g'] = self.filter_g + ft['r'] = self.filter_r + ft['i'] = self.filter_i + + # load filter name in three channels + ft_wave = dict() + ft_wave['g'] = self.filterP[ft['g']]['wave_nm'] + ft_wave['r'] = self.filterP[ft['r']]['wave_nm'] + ft_wave['i'] = self.filterP[ft['i']]['wave_nm'] + + ft_eff = dict() + ft_eff['g'] = self.filterP[ft['g']]['throughout'] + ft_eff['r'] = self.filterP[ft['r']]['throughout'] + ft_eff['i'] = self.filterP[ft['i']]['throughout'] + + chlist = ['g', 'r', 'i'] + + self.Sky_Noise = dict() + + for k in range(3): + ch = chlist[k] + tel_wavearr = teleff[ch][:, 0] + tel_effarr = teleff[ch][:, 1] + tel_effnew = np.interp(wave, tel_wavearr, tel_effarr) + + ft_wavearr = ft_wave[ch] + ft_effarr = ft_eff[ch] + ft_effnew = np.interp(wave, ft_wavearr, ft_effarr) + + self.Sky_Noise[ch] = np.sum(Nphoton_flux*ft_effnew*tel_effnew) + + return self.Sky_Noise + + ################################################## + ############################################################################## + + def cal_StarSED_img(self): + ##################### + + pixelscale = self.information['pixel_size'] # arcsec, pixel scale size + ghostOffsetX = self.information['ghostoffsetx'] + ghostOffsetY = self.information['ghostoffsety'] + ghostMx = self.information['ghostratio'] + + rotTelPos = self.information['rotTelPos'] + rotSkyPos = self.information['rotSkyPos'] + + theta = rotTelPos - rotSkyPos + + ############ load star data catlog ##################### + # 'selection_MACSJ0744.9 3927_20230517_concat.csv' + starcat = self.information['starcat'] + + if get_file_extension(starcat) == '.csv': + + # starcat='selection_20230517_concat.fits' + ################################################### + + self.log.info('Stat catlog file name is %s' % (starcat)) + + ########################################## + df2 = pandas.read_csv( + self.information['dir_path']+'MCI_inputData/star_input/'+starcat) + + df3 = df2[(abs(df2['ra_gaia']-df2['ra_gaia'].mean()) < 400/3600.0) + & (abs(df2['dec_gaia']-df2['dec_gaia'].mean()) < 400/3600.0)] + + df3.index = range(len(df3)) + + ################################################### + + # np.save('umag.npy', df3['umag'] ) + self.star = df3 + + del df2, df3 + + self.information['ra_obj'] = self.star['ra_gaia'].mean() + self.information['dec_obj'] = self.star['dec_gaia'].mean() + nsrcs = len(self.star['ra_gaia']) + + if get_file_extension(starcat) == '.fits': + + self.log.info('Stat catlog file name is %s' % (starcat)) + + df0 = fits.open( + self.information['dir_path']+'/MCI_inputData/star_input/'+starcat) + + self.star = df0[1].data + + self.information['ra_obj'] = self.star['ra'].mean() + self.information['dec_obj'] = self.star['dec'].mean() + nsrcs = len(self.star['ra']) + + ################################## + + self.information['star_ra'] = self.information['ra_obj'] + self.information['star_dec'] = self.information['dec_obj'] + + center_ra = self.information['T_disRa'] + self.information['ra_obj'] + center_dec = self.information['T_disDec'] + self.information['dec_obj'] + + self.information['ra_pnt0'] = center_ra + self.information['dec_pnt0'] = center_dec + + channel = ['g', 'r', 'i'] + ####################################################################### + #################### cal_PSF_array ################################# + + if self.save_starpsf: + + self.log.info('calculate and save star PSF data......') + + primary_g = fits.PrimaryHDU() + PSF_g = fits.HDUList([primary_g]) + + primary_r = fits.PrimaryHDU() + PSF_r = fits.HDUList([primary_r]) + + primary_i = fits.PrimaryHDU() + PSF_i = fits.HDUList([primary_i]) + + fov = 0.05*min(self.information['xsize'], + self.information['ysize']) + + dec_arr, ra_arr = make_c_coor(fov, 25) + + k = 0 + + ####################################################################### + for ii in range(len(ra_arr)): + for jj in range(len(ra_arr)): + + ################################################################## + + galRa = ra_arr[ii, jj] + center_ra * \ + 3600 # ra of PFS, arcsecond + galDec = dec_arr[ii, jj]+center_dec * \ + 3600 # dec of PSF, arcsecond + + fsx, fsy = cal_pos(center_ra, center_dec, + rotTelPos, rotSkyPos, galRa, galDec) + + ############# do field distottion ########## + + if self.distortion: + + for i in range(3): + ch = channel[i] + fpx, fpy = distortField(fsx, fsy, ch) + + else: + + fpx = fsx + fpy = fsy + + ################################################################### + + psf = dict() + + for i in range(3): + + ch = channel[i] + + psfmat = self.get_PSF(fpx, fpy, ch) + + temp = 0 + + for iwave in range(7): + temp = temp+1/7.0*psfmat[:, :, iwave] + + temp = temp/temp.sum() + + # rotate the PSF data + if abs(theta.deg) > 0: + # here we choose reshape=False, the rotated image will + psf[ch] = ndimage.rotate( + temp, theta.deg, order=1, reshape=False) + else: + psf[ch] = temp + ################################### + + hdu_g = fits.ImageHDU(psf['g']) + hdu_g.header['ra'] = ra_arr[ii, jj] + hdu_g.header['dec'] = dec_arr[ii, jj] + hdu_g.header['rot_deg'] = -theta.deg + PSF_g.append(hdu_g) + ################################### + hdu_r = fits.ImageHDU(psf['r']) + hdu_r.header['ra'] = ra_arr[ii, jj] + hdu_r.header['dec'] = dec_arr[ii, jj] + hdu_r.header['rot_deg'] = -theta.deg + PSF_r.append(hdu_r) + ################################### + hdu_i = fits.ImageHDU(psf['i']) + hdu_i.header['ra'] = ra_arr[ii, jj] + hdu_i.header['dec'] = dec_arr[ii, jj] + hdu_i.header['rot_deg'] = -theta.deg + PSF_i.append(hdu_i) + ################################### + del hdu_g + del hdu_r + del hdu_i + k = k+1 + + ############################################ + file_g = self.result_path+'/PSF_Data/'+'PSF_sim_No.' + \ + str(self.information['simnumber'])+'_C1.fits' + PSF_g.writeto(file_g, overwrite=True) + + file_r = self.result_path+'/PSF_Data/'+'PSF_sim_No.' + \ + str(self.information['simnumber'])+'_C2.fits' + PSF_r.writeto(file_r, overwrite=True) + + file_i = self.result_path+'/PSF_Data/'+'PSF_sim_No.' + \ + str(self.information['simnumber'])+'_C3.fits' + PSF_i.writeto(file_i, overwrite=True) + + del PSF_g + del PSF_r + del PSF_i + self.log.info('finish save star PSF data......') + + ########## finish save PSF fits ########################################## + ############################################################################ + + self.log.info('load star catlog successfully') + + ################################################################## + obj = dict() + obj['g'] = np.zeros(3) + obj['r'] = np.zeros(3) + obj['i'] = np.zeros(3) + + fullimg = dict() + star_input = dict() + star_input['g'] = np.zeros((nsrcs, 3)) + star_input['r'] = np.zeros((nsrcs, 3)) + star_input['i'] = np.zeros((nsrcs, 3)) + + star_output = dict() + star_output['g'] = np.zeros((nsrcs, 3)) + star_output['r'] = np.zeros((nsrcs, 3)) + star_output['i'] = np.zeros((nsrcs, 3)) + + final_image = dict() + + final_image['g'] = galsim.ImageF( + int(self.information['xsize']), int(self.information['ysize'])) + final_image['r'] = galsim.ImageF( + int(self.information['xsize']), int(self.information['ysize'])) + final_image['i'] = galsim.ImageF( + int(self.information['xsize']), int(self.information['ysize'])) + + final_image['g'].setOrigin(0, 0) + final_image['r'].setOrigin(0, 0) + final_image['i'].setOrigin(0, 0) + ####################################################################### + nlayccd = 0 + + ############################################ + + ra = self.information['ra_pnt0'] + dec = self.information['dec_pnt0'] + + time_jd = time2jd(self.dt) + + x_sat = float(self.orbit_pars[self.orbit_exp_num, 1]) + y_sat = float(self.orbit_pars[self.orbit_exp_num, 2]) + z_sat = float(self.orbit_pars[self.orbit_exp_num, 3]) + + wave0, zodi0 = self.zodiacal( + ra, dec, self.TianCe_day) # erg/s/cm^2/A/arcsec^2 + + self.log.info('dir_path:') + self.log.info(self.information['dir_path']) + + # EarthShine from straylight + sl = StrayLight(self.information['dir_path'], jtime=time_jd, sat=np.array( + [x_sat, y_sat, z_sat]), radec=np.array([(ra*u.degree).value, (dec*u.degree).value])) + earth_e = sl.caculateEarthShineFilter(filter='r') + star_e = sl.caculateStarLightFilter(filter='r') + angle_earth = earth_angle(time_jd, x_sat, y_sat, z_sat, ra, dec) + + if angle_earth < 0: + earth_e = 0 + + earthshine_wave0, earthshine_flux0 = ill2flux( + earth_e+star_e, self.information['dir_path']) + + # sample as mci wavelength + wave_mci = np.linspace(2500, 11000, 8501) # np.arange(2500, 11000, 1) + f1 = interp1d(wave0, zodi0) + zodi_mci = f1(wave_mci) + + f2 = interp1d(earthshine_wave0, earthshine_flux0) + earthshine_mci = f2(wave_mci) + + self.zodiacal_wave = wave_mci # in A + self.zodiacal_flux = zodi_mci + + self.earthshine_wave = wave_mci # A + self.earthshine_flux = earthshine_mci + + ####################################################################################### + ######################################################################################### + + self.cal_sky_noise() + + ################################# + + self.information['target'] = deg2HMS( + self.information['ra_obj'], self.information['dec_obj']) + + self.information['CRVAL1'] = self.information['ra_pnt0'] + self.information['CRVAL2'] = self.information['dec_pnt0'] + + self.information['CRPIX1'] = self.information['xsize']/2-0.5 + self.information['CRPIX2'] = self.information['ysize']/2-0.5 + + self.information['CD1_1'] = - \ + np.cos(-theta)*self.information['pixel_size']/3600.0 + self.information['CD1_2'] = np.sin(-theta) * \ + self.information['pixel_size']/3600.0 + + self.information['CD2_1'] = np.sin(-theta) * \ + self.information['pixel_size']/3600.0 + self.information['CD2_2'] = np.cos(-theta) * \ + self.information['pixel_size']/3600.0 + + ####################################################################### + if self.TianceEffect: + if get_file_extension(starcat) == '.fits': + + ra_list = self.star['ra'].tolist() + dec_list = self.star['dec'].tolist() + + pmra_list = [0.0 for i in range(len(ra_list))] + pmdec_list = [0.0 for i in range(len(ra_list))] + parallax_list = [0.0 for i in range(len(ra_list))] + rv_list = [0.0 for i in range(len(ra_list))] + else: + + ra_list = self.star['ra_gaia'].tolist() + dec_list = self.star['dec_gaia'].tolist() + pmra_list = self.star['pmra_gaia'].tolist() + pmdec_list = self.star['pmdec_gaia'].tolist() + parallax_list = self.star['parallax_gaia'].tolist() + rv_list = [0.0 for i in range(len(ra_list))] + + ################################################ + newRa, newDec = shao.onOrbitObsPosition(self.information['dir_path'], ra_list, dec_list, pmra_list, + pmdec_list, rv_list, parallax_list, len( + ra_list), + self.information['pos_x'], self.information['pos_y'], self.information['pos_z'], self.information['velocity_x'], self.information['velocity_y'], self.information['velocity_z'], "J2000", self.TianCe_day, self.TianCe_exp_start) + else: + + if get_file_extension(starcat) == '.fits': + newRa = self.star['ra'] + newDec = self.star['dec'] + + else: + + newRa = self.star['ra_gaia'] + newDec = self.star['dec_gaia'] + ###################################################### + + st_ra = [] + st_dec = [] + st_phot_C1 = [] + st_phot_C2 = [] + st_phot_C3 = [] + + st_posX_C1 = [] + st_posX_C2 = [] + st_posX_C3 = [] + + st_posY_C1 = [] + st_posY_C2 = [] + st_posY_C3 = [] + + st_magu = [] + st_magg = [] + st_magr = [] + st_magi = [] + st_magz = [] + + #################### generate star image ########## + + for j in tqdm(range(nsrcs)): # nsrcs length + + if j == 0: + self.log.info('begin iteration........') + # 33 + # if self.debug: + #self.log.info('j = %i' % (j)) + + starRa = 3600.0*(newRa[j]) # ra of star, arcsecond + starDec = 3600.0*(newDec[j]) # dec of star, arcsecond + + ################################################################### + ################################################################### + fsx, fsy = cal_pos(center_ra, center_dec, + rotTelPos, rotSkyPos, starRa, starDec) + # row number on CCD image + row = fsy/pixelscale+self.information['ysize']/2-0.5 + # col number on CCD image + col = fsx/pixelscale+self.information['xsize']/2-0.5 + + # print(col,row) + + if row >= self.information['ysize']+100 or row <= -100 or col >= self.information['xsize']+100 or col <= -100: + continue + + nlayccd = nlayccd+1 + + if self.debug: + if nlayccd > 1: + break + + ################################################################ + for i in range(3): + ### + ch = channel[i] + star_input[ch][int(nlayccd)-1, 0] = fsx # ra + star_input[ch][int(nlayccd)-1, 1] = fsy # dec + + ############# do field distottion ########## + + if self.distortion: + + for i in range(3): + + ch = channel[i] + fpx, fpy = distortField(fsx, fsy, ch) + obj[ch][0] = fpx + obj[ch][1] = fpy + + star_output[ch][int(nlayccd)-1, 0] = fpx # ra + star_output[ch][int(nlayccd)-1, 1] = fpy # dec + + else: + + fpx = fsx + fpy = fsy + + for i in range(3): + + ch = channel[i] + + obj[ch][0] = fpx + obj[ch][1] = fpy + + star_output[ch][int(nlayccd)-1, 0] = fpx # ra + star_output[ch][int(nlayccd)-1, 1] = fpy # dec + + ####################################################################### + + ####### use SED_code to generate star SED ####### + if get_file_extension(starcat) == '.fits': + umag = self.star['gmag'][j] + gmag = self.star['gmag'][j] + rmag = self.star['rmag'][j] + imag = self.star['imag'][j] + zmag = self.star['zmag'][j] + + # SED of j-th star + else: + umag = self.star['umag'][j] + gmag = self.star['gmag'][j] + rmag = self.star['rmag'][j] + imag = self.star['imag'][j] + zmag = self.star['zmag'][j] + + flag = (umag > 0 and gmag > 0 and rmag > + 0 and imag > 0 and zmag > 0) + + if flag == False: + continue + + wave = np.linspace(2500, 11000, 8501) + # Input redshift + redshift = 0 + + # Loading galaxy SED template + t = sed.Gal_Temp(self.information['dir_path']) + # Calculating the magnitude (u, g, r, i, z) of each template + t.toMag(redshift=redshift) + # Calculating flux + ugriz = np.array([umag, gmag, rmag, imag, zmag]) + star_flux = sed.Model_Galaxy_SED( + wave, ugriz, redshift, t, self.information['dir_path']) + + # cal_SED_photons(self, flux_arr): + intscales, PSF_eff_weight = self.cal_SED_photons(star_flux) + + gx = dict() + gy = dict() + + #self.log.info('magu =%f, magg=%f, magr=%f, magi=%f,magz=%f' % (umag,gmag,rmag,imag,zmag)) + for i in range(3): + ch = channel[i] + gy[ch] = obj[ch][1]/pixelscale+self.information['ysize'] / \ + 2-0.5 # row number on CCD image + gx[ch] = obj[ch][0]/pixelscale+self.information['xsize'] / \ + 2-0.5 # col number on CCD image + + #self.log.info('Channel in =%s, PosX(in pixel)=%f, PosY(in pixel)=%f' % (channel[i], gx[ch], gy[ch])) + + if self.debug: + + self.log.info('Star number on CCD is = %i' % (nlayccd)) + + if nlayccd % 10 == 0: + self.log.info('Star number on CCD is = %i' % (nlayccd)) + + ###### + + st_ra.append(np.float64(starRa/3600.0)) + st_dec.append(np.float64(starDec/3600.0)) + + st_phot_C1.append(intscales['g']) + st_phot_C2.append(intscales['r']) + st_phot_C3.append(intscales['i']) + + st_posX_C1.append(gx['g']) + st_posX_C2.append(gx['r']) + st_posX_C3.append(gx['i']) + + st_posY_C1.append(gy['g']) + st_posY_C2.append(gy['r']) + st_posY_C3.append(gy['i']) + + st_magu.append(umag) + st_magg.append(gmag) + st_magr.append(rmag) + st_magi.append(imag) + st_magz.append(zmag) + + ## self.log.info('begin cal PSF') + ###################################################################### + psf = dict() + + for i in range(3): + + ch = channel[i] + ####### + psfmat = self.get_PSF(fsx, fsy, ch) + + for iwave in range(7): + if iwave == 0: + temp = PSF_eff_weight[ch][iwave]*psfmat[:, :, iwave] + else: + temp = temp + \ + PSF_eff_weight[ch][iwave]*psfmat[:, :, iwave] + + temp = temp/temp.sum() + + # rotate the PSF data + if abs(theta.deg > 0): + # here we choose reshape=False, the rotated image will + psf[ch] = ndimage.rotate( + temp, theta.deg, order=1, reshape=True) + else: + psf[ch] = temp + + conv = psf[ch] + conv = conv/conv.sum() + + conv = conv*intscales[ch] + ################################################# + + stamp_img = galsim.Image(conv, copy=True) + stamp_img.setOrigin(0, 0) + stamp_img.scale = 0.025 + + # self.log.info('begining galsim.PhotonArray.makeFromImage') + + # self.log.info('stamp_img.array.max()= %i' % (stamp_img.array.max())) + + ###print('stamp_img.array.max(): ', stamp_img.array.max() ) + + ################################################# + + photons = galsim.PhotonArray.makeFromImage( + stamp_img, max_flux=max(1.0, stamp_img.array.max()/100000.0)) + + cx0, cy0 = centroid(stamp_img.array) + + ##self.log.info('appfat= %s'%(self.appFatt)) + + if self.appFatt: + ###self.log.info('begin appfat ........') + + # apply treering and bright fatter and diffusion; + SimpleTreeRing = galsim.SiliconSensor().simple_treerings( + amplitude=self.information['treering']) + + cx, cy = centroid(conv) + + # set subimge center as the pos + pos = galsim.PositionD(x=int(cx), y=int(cy)) + + siliconsensor = galsim.SiliconSensor( + strength=self.information['fatter'], diffusion_factor=self.information['difusion'], treering_func=SimpleTreeRing, treering_center=pos) + + image = galsim.Image(conv.shape[0], conv.shape[1]) + image.scale = 0.025 + image.setOrigin(0, 0) + + photons.x = photons.x/image.scale + photons.y = photons.y/image.scale + + siliconsensor.accumulate(photons, image) + + photons = galsim.PhotonArray.makeFromImage( + image, max_flux=max(1.0, stamp_img.array.max()/100000.0)) + + cx0, cy0 = centroid(image.array) + + ############################################ + ############################################ + photons.x = photons.x-cx0*0.025 + photons.y = photons.y-cy0*0.025 + + ########################################################## + ######################################################## + photons.x = photons.x/0.05+gx[ch] + photons.y = photons.y/0.05+gy[ch] + + photons.addTo(final_image[ch]) + + ############################### + # debug code + # fi= galsim.ImageF(int(self.information['xsize']), int(self.information['ysize'])) + # fi.setOrigin(0,0) + # photons.addTo(fi) + # icx,i cy=centroid(fi.array) + # print(icx,icy) + ############################################################### + + #### add ghost image #### + if self.ghosts: + ghostphotons = galsim.PhotonArray( + len(photons.x), x=photons.x, y=photons.y, flux=photons.flux) + ghostphotons.flux = ghostMx*ghostphotons.flux + ghostphotons.x = photons.x+ghostOffsetY + ghostphotons.y = photons.y+ghostOffsetX + + if ghostphotons.flux.max() > 0.01: + ghostphotons.addTo(final_image[ch]) + + ############################################################### + print('total star nummber on CCD is: ', nlayccd) + + self.log.info('total star nummber on CCD = %i' % (nlayccd)) + + tab = Table() + tab['st_ra'] = np.asarray(st_ra) + tab['st_dec'] = np.asarray(st_dec) + + tab['phot_C1'] = np.asarray(st_phot_C1) + tab['phot_C2'] = np.asarray(st_phot_C2) + tab['phot_C3'] = np.asarray(st_phot_C3) + + tab['posX_C1'] = np.asarray(st_posX_C1) + tab['posX_C2'] = np.asarray(st_posX_C2) + tab['posX_C3'] = np.asarray(st_posX_C3) + + tab['posY_C1'] = np.asarray(st_posY_C1) + tab['posY_C2'] = np.asarray(st_posY_C2) + tab['posY_C3'] = np.asarray(st_posY_C3) + + tab['st_magu'] = np.asarray(st_magu) + tab['st_magg'] = np.asarray(st_magg) + tab['st_magr'] = np.asarray(st_magr) + tab['st_magi'] = np.asarray(st_magi) + tab['st_magz'] = np.asarray(st_magz) + + # 33 + + file = self.result_path+'/log_Data/star_info' + \ + '_SN_'+str(self.information['simnumber'])+'.fits' + tab.write(file, overwrite=True) + + ################################################################# + ################################################################# + + fullimg['g'] = final_image['g'].array + fullimg['r'] = final_image['r'].array + fullimg['i'] = final_image['i'].array + + return fullimg + +################################################################################ +################################################################################ + + def zodiacal(self, ra, dec, time): + """ + For given RA, DEC and TIME, return the interpolated zodical spectrum in Leinert-1998. + + :param ra: RA in unit of degree, ICRS frame + :param dec: DEC in unit of degree, ICRS frame + :param time: the specified string that in ISO format i.e., yyyy-mm-dd. + :return: + wave_A: wavelength of the zodical spectrum + spec_mjy: flux of the zodical spectrum, in unit of MJy/sr + spec_erg: flux of the zodical spectrum, in unit of erg/s/cm^2/A/sr + + """ + + # get solar position + dt = datetime.fromisoformat(time) + ###jd = julian.to_jd(dt, fmt='jd') + + jd = time2jd(dt) + t = Time(jd, format='jd', scale='utc') + + astro_sun = get_sun(t) + ra_sun, dec_sun = astro_sun.gcrs.ra.deg, astro_sun.gcrs.dec.deg + + radec_sun = SkyCoord(ra=ra_sun*u.degree, + dec=dec_sun*u.degree, frame='gcrs') + lb_sun = radec_sun.transform_to('geocentrictrueecliptic') + + # get offsets between the target and sun. + radec_obj = SkyCoord(ra=ra*u.degree, dec=dec*u.degree, frame='icrs') + lb_obj = radec_obj.transform_to('geocentrictrueecliptic') + + beta = abs(lb_obj.lat.degree) + lamda = abs(lb_obj.lon.degree - lb_sun.lon.degree) + + # interpolated zodical surface brightness at 0.5 um + zodi = pd.read_csv( + self.information['dir_path']+'MCI_inputData/refs/zodi_map.dat', sep='\s+', header=None, comment='#') + beta_angle = np.array([0, 5, 10, 15, 20, 25, 30, 45, 60, 75]) + lamda_angle = np.array([0, 5, 10, 15, 20, 25, 30, 35, 40, 45, + 60, 75, 90, 105, 120, 135, 150, 165, 180]) + xx, yy = np.meshgrid(beta_angle, lamda_angle) + f = interpolate.interp2d(xx, yy, zodi, kind='linear') + zodi_obj = f(beta, lamda) # 10^�? W m�? sr�? um�? + + # read the zodical spectrum in the ecliptic + cat_spec = pd.read_csv( + self.information['dir_path']+'MCI_inputData/refs/solar_spec.dat', sep='\s+', header=None, comment='#') + wave = cat_spec[0].values # A + spec0 = cat_spec[1].values # 10^-8 W m^�? sr^�? μm^�? + zodi_norm = 252 # 10^-8 W m^�? sr^�? μm^�? + + spec = spec0 * (zodi_obj / zodi_norm) * 1e-8 # W m^�? sr^�? μm^�? + + # convert to the commonly used unit of MJy/sr, erg/s/cm^2/A/sr + wave_A = wave # A + # spec_mjy = spec * 0.1 * wave_A**2 / 3e18 * 1e23 * 1e-6 # MJy/sr + spec_erg = spec * 0.1 # erg/s/cm^2/A/sr + # erg/s/cm^2/A/arcsec^2 + spec_erg2 = spec_erg / 4.25452e10 + + self.zodiacal_wave = wave_A # in A + + self.zodiacal_flux = spec_erg2 + + return wave_A, spec_erg2 + ################################################################################## + + ########################################################################## + def cal_Lensing_galaxy_img(self): + ##################### + + # arcsec, pixel scale size; + self.dsx_arc = self.information['pixel_size'] + # arcsec, field size; + self.bsz_arc = max( + self.information['ysize'], self.information['xsize'])*self.dsx_arc + self.bsz_deg = self.bsz_arc/3600. + + ############################################################# + # galaxy catlog information + ra = 74.0563449685417 + dec = -34.155241212932374 + + # halo_id = 229600100382 + if self.sim_star == False: + self.information['star_ra'] = ra + self.information['star_dec'] = dec + + self.information['gal_ra'] = ra + self.information['gal_dec'] = dec + + self.information['ra_obj'] = self.information['star_ra'] + self.information['dec_obj'] = self.information['star_dec'] + + losscale = self.information['pixel_size'] + rotTelPos = self.information['rotTelPos'] + rotSkyPos = self.information['rotSkyPos'] + + theta = rotTelPos - rotSkyPos + + center_ra = self.information['T_disRa'] + self.information['star_ra'] + center_dec = self.information['T_disDec'] + \ + self.information['star_dec'] + + self.information['ra_pnt0'] = center_ra + self.information['dec_pnt0'] = center_dec + #################################################################### + ################################################################################## + + time_jd = time2jd(self.dt) + + x_sat = float(self.orbit_pars[self.orbit_exp_num, 1]) + y_sat = float(self.orbit_pars[self.orbit_exp_num, 2]) + z_sat = float(self.orbit_pars[self.orbit_exp_num, 3]) + + wave0, zodi0 = self.zodiacal( + self.information['ra_pnt0'], self.information['dec_pnt0'], self.TianCe_day) # erg/s/cm^2/A/arcsec^2 + + # EarthShine from straylight + sl = StrayLight(self.information['dir_path'], jtime=time_jd, sat=np.array([x_sat, y_sat, z_sat]), + radec=np.array([(ra*u.degree).value, (dec*u.degree).value])) + + earth_e = sl.caculateEarthShineFilter(filter='r') + star_e = sl.caculateStarLightFilter(filter='r') + angle_earth = earth_angle(time_jd, x_sat, y_sat, z_sat, ra, dec) + + if angle_earth < 0: + earth_e = 0 + + earthshine_wave0, earthshine_flux0 = ill2flux( + earth_e+star_e, self.information['dir_path']) + + # sample as ifs wavelength + wave_mci = np.linspace(2500, 11000, 8501) # np.arange(2500, 11000, 1) + f1 = interp1d(wave0, zodi0) + zodi_mci = f1(wave_mci) + + f2 = interp1d(earthshine_wave0, earthshine_flux0) + earthshine_mci = f2(wave_mci) + + # df = pd.DataFrame({'wave': wave_ifs, + # 'zodi': zodi_ifs, + # 'earthshine': earthshine_ifs}) + + self.zodiacal_wave = wave_mci # in A + self.zodiacal_flux = zodi_mci + + self.earthshine_wave = wave_mci # A + self.earthshine_flux = earthshine_mci + + ######################################################################################## + + ################################################################################## + self.cal_sky_noise() + + ############ load galaxy data with SED ############################ + + self.information['target'] = deg2HMS( + self.information['star_ra'], self.information['star_dec']) + + self.information['CRVAL1'] = self.information['ra_pnt0'] + self.information['CRVAL2'] = self.information['dec_pnt0'] + + self.information['CRPIX1'] = self.information['xsize']/2-0.5 + self.information['CRPIX2'] = self.information['ysize']/2-0.5 + + self.information['CD1_1'] = - \ + np.cos(-theta)*self.information['pixel_size']/3600.0 + self.information['CD1_2'] = np.sin(-theta) * \ + self.information['pixel_size']/3600.0 + + self.information['CD2_1'] = np.sin(-theta) * \ + self.information['pixel_size']/3600.0 + self.information['CD2_2'] = np.cos(-theta) * \ + self.information['pixel_size']/3600.0 + + obj = dict() + obj['g'] = np.zeros(3) + obj['r'] = np.zeros(3) + obj['i'] = np.zeros(3) + + fullimg = dict() + galaxy_input = dict() + + nsrcs = 72000 + + galaxy_input['g'] = np.zeros((nsrcs, 3)) + galaxy_input['r'] = np.zeros((nsrcs, 3)) + galaxy_input['i'] = np.zeros((nsrcs, 3)) + + galaxy_output = dict() + galaxy_output['g'] = np.zeros((nsrcs, 3)) + galaxy_output['r'] = np.zeros((nsrcs, 3)) + galaxy_output['i'] = np.zeros((nsrcs, 3)) + + channel = ['g', 'r', 'i'] + final_image = dict() + + final_image['g'] = galsim.ImageF( + int(self.information['xsize']), int(self.information['ysize'])) + final_image['r'] = galsim.ImageF( + int(self.information['xsize']), int(self.information['ysize'])) + final_image['i'] = galsim.ImageF( + int(self.information['xsize']), int(self.information['ysize'])) + + final_image['g'].setOrigin(0, 0) + final_image['r'].setOrigin(0, 0) + final_image['i'].setOrigin(0, 0) + + #################################################################### + + gal_ra = [] # .append(np.float64(galRa/3600.0)) + gal_dec = [] # .append(np.float64(galDec/3600.0)) + + phot_C1 = [] # .append(intscales['g']) + phot_C2 = [] # .append(intscales['r']) + phot_C3 = [] # .append(intscales['i']) + + posX_C1 = [] # .append(gx['g']) + posX_C2 = [] # .append(gx['r']) + posX_C3 = [] # .append(gx['i']) + + posY_C1 = [] # .append(gy['g']) + posY_C2 = [] # .append(gy['r']) + posY_C3 = [] # .append(gy['i']) + + galID = [] + haloID = [] + + gal_new_ra = [] + gal_new_dec = [] + + gal_ori_ra = [] + gal_ori_dec = [] + + gal_rescale = [] + gal_szs = [] + + gal_magu = [] + gal_magg = [] + gal_magr = [] + gal_magi = [] + gal_magz = [] + + ######################################## + + nlayccd = 0 + + # generate loc image + + if self.debug: + nk2 = 1 + else: + nk2 = 1 + + for k2 in (range(nk2)): + print('k2=', k2) + + # Lens_SED_IMG_0.025_230606 + + if self.lensing: + filename = self.information['dir_path'] + \ + 'MCI_inputData/galaxy_Input/Lens_SED_IMG_0.025_230606/Lens_img_cut_IMG_' + \ + str(k2+1)+'.fits' + + else: + filename = self.information['dir_path'] + \ + 'MCI_inputData/galaxy_Input/noLens_SED_IMG_0.025_230626/Lens_img_cut_IMG_' + \ + str(k2+1)+'.fits' + + self.log.info('galaxy_Input image path is: %s' % (filename)) + + if not os.path.exists(filename): + print('finish load all the input galaxy image fits files') + break + + srcs_cat = fits.open(filename) + + #### load galaxy SED fitsfile ### + + if self.lensing: + + filename = self.information['dir_path'] + \ + 'MCI_inputData/galaxy_Input/Lens_SED_IMG_0.025_230606/Lens_img_cut_SED_' + \ + str(k2+1)+'.fits' + + else: + + filename = self.information['dir_path'] + \ + 'MCI_inputData/galaxy_Input/noLens_SED_IMG_0.025_230626/Lens_img_cut_SED_' + \ + str(k2+1)+'.fits' + srcs_sed = fits.open(filename) + + self.log.info('galaxy_Input SED path is: %s' % (filename)) + + ###### do Tiance Effect ### + ra_list = [] + dec_list = [] + + for kkk in range(1, len(srcs_cat)): + t1 = srcs_cat[kkk].header['new_ra'] - \ + self.information['gal_ra'] + self.information['star_ra'] + ra_list.append(float(t1)) + + t2 = srcs_cat[kkk].header['new_dec'] - \ + self.information['gal_dec'] + self.information['star_dec'] + dec_list.append(float(t2)) + + if self.TianceEffect: + + pmra_list = [0.0 for i in range(len(ra_list))] + pmdec_list = [0.0 for i in range(len(ra_list))] + parallax_list = [0.0 for i in range(len(ra_list))] + rv_list = [0.0 for i in range(len(ra_list))] + + ################################################ + newRa, newDec = shao.onOrbitObsPosition(self.information['dir_path'], ra_list, dec_list, pmra_list, + pmdec_list, rv_list, parallax_list, len( + ra_list), + self.information['pos_x'], self.information['pos_y'], self.information['pos_z'], self.information['velocity_x'], self.information['velocity_y'], self.information['velocity_z'], "J2000", self.TianCe_day, self.TianCe_exp_start) + else: + newRa = ra_list + newDec = dec_list + + ############################################################################ + ################################################################### + for k1 in tqdm(range(len(newRa))): + + #print('k1=', k1) + + # *(srcs_cat[k1].header['new_ra'] -self.information['gal_ra'] +self.information['star_ra']) # ra of galaxies, arcsecond + galRa = newRa[k1]*3600 + # (srcs_cat[k1].header['new_dec']-self.information['gal_dec']+self.information['star_dec']) # dec of galaxies, arcsecond + galDec = newDec[k1]*3600 + + ########################################################################### + ############################################################################ + ### test code #### + + ################################################################# + fsx, fsy = cal_pos(center_ra, center_dec, + rotTelPos, rotSkyPos, galRa, galDec) + + if abs(fsx) > 235 or abs(fsy) > 235: + continue + + # # SED of j-th galaxy ,# unit of 10-17 erg/s/A/cm2 + # here is k1+1, not k1, k1 begins with 0 + gal_flux = srcs_sed[k1+1].data + + ################################ + ### rotate the lensed_images_g ### + if abs(theta.deg) > 0: + # here we choose reshape=False, the rotated image will + lensed_images_g = ndimage.rotate( + srcs_cat[k1+1].data, -theta.deg, order=1, reshape=True) + else: + lensed_images_g = srcs_cat[k1+1].data + + if lensed_images_g.sum() <= 0: + continue + + ################################################################ + for i in range(3): + ### + ch = channel[i] + galaxy_input[ch][int(nlayccd)-1, 0] = fsx # ra + galaxy_input[ch][int(nlayccd)-1, 1] = fsy # dec + # galaxy_input[ch][int(nlayccd)-1, 2] =mag[ch] # mag g + + ############# do field distottion ########## + + if self.distortion: + + for i in range(3): + + ch = channel[i] + + fpx, fpy = distortField(fsx, fsy, ch) + obj[ch][0] = fpx + obj[ch][1] = fpy + # obj[ch][2]=mag[ch] + + galaxy_output[ch][int(nlayccd)-1, 0] = fpx # ra + galaxy_output[ch][int(nlayccd)-1, 1] = fpy # dec + # galaxy_output[ch][int(nlayccd)-1, 2] =mag[ch] # mag g + + else: + + fpx = fsx + fpy = fsy + + for i in range(3): + + ch = channel[i] + + obj[ch][0] = fpx + obj[ch][1] = fpy + # obj[ch][2]=mag[ch] + + galaxy_output[ch][int(nlayccd)-1, 0] = fpx # ra + galaxy_output[ch][int(nlayccd)-1, 1] = fpy # dec + # galaxy_output[ch][int(nlayccd)-1, 2] =mag[ch] # mag g + + ####################################################################### + #self.log.info('Galaxy number=%i, Ra(in arcsec)=%f, Dec(in arcsec)=%f' % (nlayccd, fpx, fpy)) + #self.log.info('Galaxy number=%i, Mag_g=%f, Mag_r=%f,Mag_i=%f' % (nlayccd, mag['g'],mag['r'],mag['i'])) + + #print('Flux sum =', gal_flux.sum()) + + # cal_SED_photons(self, flux_arr): + intscales, PSF_eff_weight = self.cal_SED_photons(gal_flux) + + ls_imgsum = lensed_images_g.sum() + + if intscales['g'].max()/ls_imgsum < 0.1 and intscales['r'].max()/ls_imgsum < 0.1 and intscales['i'].max()/ls_imgsum < 0.1: + ######### + continue + ######## test code ########## + + if self.debug: + if nlayccd > 1: + ######### + break + + # 33 + + img = dict() + + lensed_images_g = lensed_images_g/ls_imgsum + img['g'] = lensed_images_g*intscales['g'] + img['r'] = lensed_images_g*intscales['r'] + img['i'] = lensed_images_g*intscales['i'] + + if img['g'].max() < 0.01 and img['r'].max() < 0.01 and img['i'].max() < 0.01: + ######### + continue + + nlayccd = nlayccd+1 + + ################################################ + + gx = dict() + gy = dict() + + for i in range(3): + + ch = channel[i] + + gy[ch] = obj[ch][1]/losscale+self.information['ysize'] / \ + 2-0.5 # row number on CCD image + gx[ch] = obj[ch][0]/losscale+self.information['xsize'] / \ + 2-0.5 # col number on CCD image + + ###################################################################### + + gal_ra.append(np.float64(galRa/3600.0)) + gal_dec.append(np.float64(galDec/3600.0)) + + phot_C1.append(intscales['g']) + phot_C2.append(intscales['r']) + phot_C3.append(intscales['i']) + + posX_C1.append(gx['g']) + posX_C2.append(gx['r']) + posX_C3.append(gx['i']) + + posY_C1.append(gy['g']) + posY_C2.append(gy['r']) + posY_C3.append(gy['i']) + + gal_magu.append(srcs_cat[k1+1].header['magu']) + gal_magg.append(srcs_cat[k1+1].header['magg']) + gal_magr.append(srcs_cat[k1+1].header['magr']) + gal_magi.append(srcs_cat[k1+1].header['magi']) + gal_magz.append(srcs_cat[k1+1].header['magz']) + + galID.append(srcs_cat[k1+1].header['GALAXID']) + haloID.append(srcs_cat[k1+1].header['HALOID']) + + gal_new_ra.append(srcs_cat[k1+1].header['NEW_RA']) + gal_new_dec.append(srcs_cat[k1+1].header['NEW_DEC']) + + gal_ori_ra.append(srcs_cat[k1+1].header['OLD_RA']) + gal_ori_dec.append(srcs_cat[k1+1].header['OLD_DEC']) + + gal_rescale.append(srcs_cat[k1+1].header['RESCALE']) + gal_szs.append(srcs_cat[k1+1].header['SZS']) + + #self.log.info('Galaxy number=%i, Mag_g=%f, Mag_r=%f,Mag_i=%f' % (nlayccd, mag['g'],mag['r'],mag['i'])) + + if nlayccd % 1000 == 0: + self.log.info('Galaxy number on CCD is = %i' % (nlayccd)) + + psf = dict() + + for i in range(3): + + ch = channel[i] + + psfmat = self.get_PSF(fsx, fsy, ch) + + temp = 0 + + for iwave in range(7): + temp = temp + \ + PSF_eff_weight[ch][iwave]*psfmat[:, :, iwave] + + temp = temp/temp.sum() + + # rotate the PSF data + if abs(theta.deg) > 0: + # here we choose reshape=False, the rotated image will + psf[ch] = ndimage.rotate( + temp, theta.deg, order=1, reshape=False) + else: + psf[ch] = temp + + conv = fftconvolve(img[ch], psf[ch], mode='full') + # continue + # suppress negative numbers + + conv[conv < 0.0] = 0.0 + + ################################################# + stamp_img = galsim.Image(conv, copy=True) + stamp_img.setOrigin(0, 0) + stamp_img.scale = 0.025 + + photons = galsim.PhotonArray.makeFromImage( + stamp_img, max_flux=1.0) + + cx0, cy0 = centroid(stamp_img.array) + ########### apply fat effect to galaxy image ##### + if self.appFatt: + + # apply treering and bright fatter and diffusion; + + SimpleTreeRing = galsim.SiliconSensor().simple_treerings( + amplitude=self.information['treering']) + + cx, cy = centroid(conv) + + # set subimge center as the pos + pos = galsim.PositionD(x=int(cx), y=int(cy)) + + siliconsensor = galsim.SiliconSensor( + strength=self.information['fatter'], diffusion_factor=self.information['difusion'], treering_func=SimpleTreeRing, treering_center=pos) + + image = galsim.Image(conv.shape[1], conv.shape[0]) + image.scale = 0.025 + image.setOrigin(0, 0) + + photons.x = photons.x/image.scale + photons.y = photons.y/image.scale + + siliconsensor.accumulate(photons, image) + + photons = galsim.PhotonArray.makeFromImage( + image, max_flux=1.0) + + cx0, cy0 = centroid(image.array) + ############################################################## + ############################################################## + ################################################################## + ################################################################## + + photons.x = photons.x-cx0*0.025 + photons.y = photons.y-cy0*0.025 + + photons.x = photons.x/0.05+gx[ch] + photons.y = photons.y/0.05+gy[ch] + + photons.addTo(final_image[ch]) + + #################################################### + ############################################ + ############################################ + ######################################################################## + + print('total los nummber on CCD is: ', nlayccd) + + self.log.info('total los nummber on CCD is = %i' % (nlayccd)) + + ################################################################# + tab = Table() + tab['gal_ra'] = np.asarray(gal_ra) + tab['gal_dec'] = np.asarray(gal_dec) + + tab['phot_C1'] = np.asarray(phot_C1) + tab['phot_C2'] = np.asarray(phot_C2) + tab['phot_C3'] = np.asarray(phot_C3) + + tab['posX_C1'] = np.asarray(posX_C1) + tab['posX_C2'] = np.asarray(posX_C2) + tab['posX_C3'] = np.asarray(posX_C3) + + tab['posY_C1'] = np.asarray(posY_C1) + tab['posY_C2'] = np.asarray(posY_C2) + tab['posY_C3'] = np.asarray(posY_C3) + + tab['gal_magu'] = np.asarray(gal_magu) + tab['gal_magg'] = np.asarray(gal_magg) + tab['gal_magr'] = np.asarray(gal_magr) + tab['gal_magi'] = np.asarray(gal_magi) + tab['gal_magz'] = np.asarray(gal_magz) + + tab['gal_ID'] = np.asarray(galID) + tab['halo_ID'] = np.asarray(haloID) + + tab['gal_new_ra'] = np.asarray(gal_new_ra) + tab['gal_new_dec'] = np.asarray(gal_new_dec) + + tab['gal_ori_ra'] = np.asarray(gal_ori_ra) + tab['gal_ori_dec'] = np.asarray(gal_ori_dec) + + tab['gal_rescale'] = np.asarray(gal_rescale) + tab['gal_szs'] = np.asarray(gal_szs) + + #################################################################### + + file = self.result_path+'/log_Data/gal_info' + \ + '_SN_'+str(self.information['simnumber'])+'.fits' + tab.write(file, overwrite=True) + + ####################################################################### + + fullimg['g'] = final_image['g'].array + fullimg['r'] = final_image['r'].array + fullimg['i'] = final_image['i'].array + + return fullimg + + +############################################################################### + +############################################################################### + + + def generatePRNU(self, ave=1.0, sigma=0.01): + """ + Creates a PRNU field image with given properties. + + :return: PRNU field image + :rtype: ndarray + """ + self.log.info('Generating a flat field...') + self.log.info( + 'The flat field has mean value of 1 and a given fluctuations, usually either 1 or 2 percent defined by sigma= %d...' % sigma) + + self.PRNU = dict() + + np.random.seed(self.filterP[self.filter_g]['seed']) + self.PRNU['g'] = np.random.normal(loc=ave, scale=sigma, size=( + self.information['ysize'], self.information['xsize'])) + + np.random.seed(self.filterP[self.filter_r]['seed']) + self.PRNU['r'] = np.random.normal(loc=ave, scale=sigma, size=( + self.information['ysize'], self.information['xsize'])) + + np.random.seed(self.filterP[self.filter_i]['seed']) + self.PRNU['i'] = np.random.normal(loc=ave, scale=sigma, size=( + self.information['ysize'], self.information['xsize'])) + + return self.PRNU +############################################################################### +############################################################################### + + def MakeFlatMatrix(self, img, seed): + #### + ysize, xsize = img.shape + np.random.seed(seed) + r1, r2, r3, r4 = np.random.random(4) + a1 = -0.5 + 0.2*r1 + a2 = -0.5 + 0.2*r2 + a3 = r3+5 + a4 = r4+5 + xmin, xmax, ymin, ymax = 0, xsize, 0, ysize + Flty, Fltx = np.mgrid[ymin:ymax, xmin:xmax] + np.random.seed(seed) + p1, p2, bg = np.random.poisson(1000, 3) + Fltz = 1e-6*(a1 * (Fltx-p1) ** 2 + a2 * (Flty-p2) + ** 2 - a3*Fltx - a4*Flty) + bg*20 + FlatMat = Fltz/np.mean(Fltz) + + return FlatMat +#################################################################################################### + + def applyFlat(self): + + # apply Flat field to image + self.FlatMat = dict() + ### + self.FlatMat['g'] = self.MakeFlatMatrix( + self.image_g, self.filterP[self.filter_g]['seed']) + self.image_g *= self.FlatMat['g'] + ### + self.FlatMat['r'] = self.MakeFlatMatrix( + self.image_r, self.filterP[self.filter_r]['seed']) + self.image_r *= self.FlatMat['r'] + ### + self.FlatMat['i'] = self.MakeFlatMatrix( + self.image_i, self.filterP[self.filter_i]['seed']) + self.image_g *= self.FlatMat['i'] + return + +###################################################################################################### + def applyPRNUeffect(self): + """ + Applies multiplicative flat field to emulate pixel-to-pixel non-uniformity. + + Because the pixel-to-pixel non-uniformity effect (i.e. multiplicative) flat fielding takes place + before CTI and other effects, the flat field file must be the same size as the pixels that see + the sky. + """ + # generate flatfiledfile, with + flatsigma = self.information['flatsigma'] + + self.generatePRNU(ave=1.0, sigma=flatsigma) + self.image_g *= self.PRNU['g'] + self.image_r *= self.PRNU['r'] + self.image_i *= self.PRNU['i'] + self.log.info('Applied flatfield to images.') + +################################################################################################## + +################################################################################################## + + def addCosmicRays(self): + """ + Add cosmic rays to the arrays based on a power-law intensity distribution for tracks. + Cosmic ray properties (such as location and angle) are chosen from random Uniform distribution. + For details, see the documentation for the cosmicrays class in the support package. + """ + self.readCosmicRayInformation() + + # to scale the number of cosmics with exposure time + self.cr['exptime'] = self.information['exptime'] + + # cosmic ray image + crImage = np.zeros( + (self.information['ysize'], self.information['xsize']), dtype=float) + + # cosmic ray instance + cosmics_g = cosmicrays.cosmicrays( + self.log, crImage, self.information['exptime'], crInfo=self.cr) + cosmics_r = cosmicrays.cosmicrays( + self.log, crImage, self.information['exptime'], crInfo=self.cr) + cosmics_i = cosmicrays.cosmicrays( + self.log, crImage, self.information['exptime'], crInfo=self.cr) + + CCD_cr_g = cosmics_g.addUpToFraction( + self.information['coveringfraction']*self.information['exptime']/300.0, limit=None) + CCD_cr_r = cosmics_r.addUpToFraction( + self.information['coveringfraction']*self.information['exptime']/300.0, limit=None) + CCD_cr_i = cosmics_i.addUpToFraction( + self.information['coveringfraction']*self.information['exptime']/300.0, limit=None) + + # paste the information + #self.image += CCD_cr + self.image_g += CCD_cr_g + self.image_r += CCD_cr_r + self.image_i += CCD_cr_i + + if self.save_cosmicrays and self.cosmicRays: + + self.log.info('Saved the cosmicRays fits...') + fits.writeto(self.result_path+'/ori_Sky/cosmicMap_C1_N_'+str( + self.information['simnumber'])+'.fits', np.int32(CCD_cr_g), overwrite=True) + fits.writeto(self.result_path+'/ori_Sky/cosmicMap_C2_N_'+str( + self.information['simnumber'])+'.fits', np.int32(CCD_cr_r), overwrite=True) + fits.writeto(self.result_path+'/ori_Sky/cosmicMap_C3_N_'+str( + self.information['simnumber'])+'.fits', np.int32(CCD_cr_i), overwrite=True) + + # count the covering factor + #area_cr = np.count_nonzero(self.cosmicMap) + area_cr_g = np.count_nonzero(CCD_cr_g) + area_cr_r = np.count_nonzero(CCD_cr_r) + area_cr_i = np.count_nonzero(CCD_cr_i) + + #self.log.info('The cosmic ray covering factor is %i pixels ' % area_cr) + self.log.info( + 'The cosmic ray in C1 channel covering factor is %i pixels ' % area_cr_g) + self.log.info( + 'The cosmic ray in C2 channel covering factor is %i pixels ' % area_cr_r) + self.log.info( + 'The cosmic ray in C3 channel covering factor is %i pixels ' % area_cr_i) + +######################################################################################## + + def ShutterEffectMat(self, img, t_exp=300, t_shutter=1.3, dist_bearing=4000, dt=0.001): + # Generate Shutter-Effect normalized image + # t_shutter: time of shutter movement + # dist_bearing: distance between two bearings of shutter leaves + # dt: delta_t of sampling + + SampleNumb = int(t_shutter/dt+1) + + DistHalf = dist_bearing/2 + + t = np.arange(SampleNumb)*dt + a_arr = 5.84*np.sin(2*math.pi/t_shutter*t) + v = np.zeros(SampleNumb) + theta = np.zeros(SampleNumb) + x = np.arange(SampleNumb)/(SampleNumb-1)*dist_bearing + s = np.zeros(SampleNumb) + s1 = np.zeros(SampleNumb) + s2 = np.zeros(SampleNumb) + brt = np.zeros(SampleNumb) + idx = np.arange(SampleNumb) + sidx = np.zeros(SampleNumb) + s1idx = np.zeros(SampleNumb) + s2idx = np.zeros(SampleNumb) + + v[0] = 0 + theta[0] = 0 + + for i in range(SampleNumb-1): + v[i+1] = v[i]+a_arr[i]*dt + theta[i+1] = theta[i]+v[i]*dt + s1[i] = DistHalf*np.cos(theta[i]) + s2[i] = dist_bearing-DistHalf*np.cos(theta[i]) + s1idx[i] = int(s1[i]/dist_bearing*(SampleNumb)) + s2idx[i] = int(s2[i]/dist_bearing*(SampleNumb)) + brt[(idx > s1idx[i]) & (idx < s2idx[i])] += dt + + if t_exp > t_shutter*2: + brt = brt*2+(t_exp-t_shutter*2) + else: + brt = brt*2 + + x = (x-dist_bearing/2)*100 + + intp = interpolate.splrep(x, brt, s=0) + + xmin = 0 # + xmax = img.shape[1] # row + ymin = 0 # + ymax = img.shape[0] # column + if xmin < np.min(x) or xmax > np.max(x): + raise LookupError("Out of focal-plane bounds in X-direction.") + if ymin < 0 or ymax > 25331: + raise LookupError("Out of focal-plane bounds in Y-direction.") + + sizex = xmax-xmin + sizey = ymax-ymin + xnewgrid = np.mgrid[xmin:(xmin+sizex)] + expeffect = interpolate.splev(xnewgrid, intp, der=0) + expeffect /= np.max(expeffect) + exparrnormal = np.tile(expeffect, (sizey, 1)) + + # Image *= exparrnormal + + shutter = (1-exparrnormal)*t_exp + shutter = shutter/shutter.max()*(t_shutter*2) + + return shutter + +######################################################################################## + def addshuttereffect(self): + ''' apply shutter effect to image''' + # + if self.information['exptime'] > 0: + self.shutter = self.ShutterEffectMat(self.image_g) + + self.shutterMat = 1-self.shutter/self.information['exptime'] + + self.image_g *= self.shutterMat + self.image_r *= self.shutterMat + self.image_i *= self.shutterMat + if self.information['simnumber'] == 1: + fits.writeto(self.result_path+'/cali_Data/' + + 'shutter.fits', self.shutter, overwrite=True) + + return + +######################################################################################## + def applyDarkCurrent(self): + """ + Apply dark current. Scales the dark with the exposure time. + + Additionally saves the image without noise to a FITS file. + """ + + # add dark + dark = self.information['exptime'] * self.information['dark'] + + self.image_g += dark + self.image_r += dark + self.image_i += dark + self.log.info('Added dark current = %f' % dark) + + ########################################################################### + ########################################################################### + def applyPoissonSkyNoise(self): + """ + Add Poisson sky noise to the image. + """ + skynoise_g = self.Sky_Noise['g']*np.ones_like(self.image_g) + + np.random.seed(11*self.information['simnumber']) + self.image_g = self.image_g + np.random.poisson(lam=skynoise_g) + self.log.info('Added Poisson sky noise on channel g') + + skynoise_r = self.Sky_Noise['r']*np.ones_like(self.image_r) + np.random.seed(111*self.information['simnumber']) + self.image_r = self.image_r + np.random.poisson(lam=skynoise_r) + self.log.info('Added Poisson sky noise on channel r') + + skynoise_i = self.Sky_Noise['i']*np.ones_like(self.image_i) + np.random.seed(1111*self.information['simnumber']) + self.image_i = self.image_i + np.random.poisson(lam=skynoise_i) + self.log.info('Added Poisson sky noise on channel i') + + return + +############################################################################### +############################################################################### + def applyCosmetics(self): + """ + Apply cosmetic defects described in the input file. + + #Number of hot and dead pixels from MSSL/Euclid/TR/12003 Issue 2 Draft b + + .. Warning:: This method does not work if the input file has exactly one line. + """ + ####################################################################### + cosmetics = np.loadtxt( + self.information['dir_path']+'MCI_inputData/data/Cosmetics_g.txt') + + x = np.round(cosmetics[:, 0]).astype(int) # row number + y = np.round(cosmetics[:, 1]).astype(int) # col number + value = np.round(cosmetics[:, 2]).astype(int) + + # cosmetics_g=np.zeros((4636,235526)) + + self.log.info('Adding cosmetic defects to C1 channel:') + for xc, yc, val in zip(x, y, value): + if 0 < xc < 4616 and 27 < (yc % 1499) < 1179: + + self.image_g[xc, yc] = val + # cosmetics_g[yc,xc]=val + self.log.info('x=%i, y=%i, value=%f' % (xc, yc, val)) + ###################################################################### + ####################################################################### + cosmetics = np.loadtxt( + self.information['dir_path']+'MCI_inputData/data/Cosmetics_r.txt') + + x = np.round(cosmetics[:, 0]).astype(int) # row number + y = np.round(cosmetics[:, 1]).astype(int) # col number + value = np.round(cosmetics[:, 2]).astype(int) + + self.log.info('Adding cosmetic defects to C2 channel:') + for xc, yc, val in zip(x, y, value): + if 0 < xc < 4616 and 27 < (yc % 1499) < 1179: + + self.image_r[xc, yc] = val + # cosmetics_r[yc,xc]=val + self.log.info('x=%i, y=%i, value=%f' % (xc, yc, val)) +############################################################################## + ####################################################################### + cosmetics = np.loadtxt( + self.information['dir_path']+'MCI_inputData/data/Cosmetics_i.txt') + + x = np.round(cosmetics[:, 0]).astype(int) # row number + y = np.round(cosmetics[:, 1]).astype(int) # col number + value = np.round(cosmetics[:, 2]).astype(int) + + # cosmetics_i=np.zeros((4636,235526)) + + self.log.info('Adding cosmetic defects to C3 channel:') + for xc, yc, val in zip(x, y, value): + if 0 < xc < 4616 and 27 < (yc % 1499) < 1179: + + self.image_i[xc, yc] = val + # cosmetics_i[yc,xc]=val + self.log.info('x=%i, y=%i, value=%f' % (xc, yc, val)) + + ####### save cosmetics image code #### + # if self.simnumber<2: + # #save cosmic ray image map + + # obsid=self.simnumber + # bluefile=self.result_path+'/calibration_Data/IFS_cosmetics'+'_'+str(obsid)+'_B.fits' + # redfile =self.result_path+'/calibration_Data/IFS_cosmetics'+'_'+str(obsid)+'_R.fits' + + # fits.writeto(bluefile, cosmetics_b, overwrite=True) + # fits.writeto(redfile, cosmetics_r, overwrite=True) + + # #output information to a FITS file + + # self.log.info('Save cosmetics images fits file to calibration_Data file ') +############################################################################################################################# + def applyRadiationDamage(self): + """ + Applies CDM03 radiation model to the image being constructed. + + .. seealso:: Class :`CDM03` + """ + # save image without CTI + self.noCTI_g = self.image_g.copy() + + self.log.debug('Starting to apply radiation damage model...') + # at this point we can give fake data... + cti = CDM03bidir(self.information, [], log=self.log) + # here we need the right input data + self.image_g = cti.applyRadiationDamage(self.image_g.copy( + ).transpose(), iquadrant=self.information['quadrant']).transpose() + self.log.info('Radiation damage added.') + + ################################################# + + self.noCTI_r = self.image_r.copy() + + self.log.debug('Starting to apply radiation damage model...') + # at this point we can give fake data... + cti = CDM03bidir(self.information, [], log=self.log) + # here we need the right input data + self.image_r = cti.applyRadiationDamage(self.image_r.copy( + ).transpose(), iquadrant=self.information['quadrant']).transpose() + self.log.info('Radiation damage added.') + + ################################################## + + self.noCTI_i = self.image_i.copy() + + self.log.debug('Starting to apply radiation damage model...') + # at this point we can give fake data... + cti = CDM03bidir(self.information, [], log=self.log) + # here we need the right input data + self.image_i = cti.applyRadiationDamage(self.image_i.copy( + ).transpose(), iquadrant=self.information['quadrant']).transpose() + self.log.info('Radiation damage added.') + + ########################################################################### + + def applyNonlinearity(self): + """ + Applies a CCD273 non-linearity model to the image being constructed. + + """ + + self.log.debug('Starting to apply non-linearity model...') + self.image_g = CCDnonLinearityModel(self.image_g.copy()) + self.log.info('Non-linearity effects included.') + + ######################################################################## + self.log.debug('Starting to apply non-linearity model...') + self.image_r = CCDnonLinearityModel(self.image_r.copy()) + self.log.info('Non-linearity effects included.') + + ######################################################################## + self.log.debug('Starting to apply non-linearity model...') + self.image_i = CCDnonLinearityModel(self.image_i.copy()) + + self.log.info('Non-linearity effects included.') + +############################################################################### + + def applyReadoutNoise(self): + """ + Applies readout noise to the image being constructed. + + The noise is drawn from a Normal (Gaussian) distribution with average=0.0 and std=readout noise. + """ + if self.overscans and self.information['xsize'] == 9216 and self.information['ysize'] == 9232: + xmin = dict() + xmax = dict() + ymin = dict() + ymax = dict() + + # calculate the non-over subimg of the final matrix of the image + for k in range(1, 17): + xmin[k] = 0 + xmax[k] = 4616+self.overscan + + ymin[k] = (1152+self.prescan+self.overscan)*(k-1) + ymax[k] = (1152+self.prescan+self.overscan)*(k) + ############### + np.random.seed(2*self.information['simnumber']) + noise_g = np.random.normal(loc=0.0, scale=self.information['g_rdnois'+str( + k)], size=(4616+self.overscan, 1152+self.overscan+self.prescan)) + self.image_g[xmin[k]:xmax[k], ymin[k]:ymax[k]] += noise_g + self.log.info('Bias of %i counts were added to the g band image' % + self.information['g_rdnois'+str(k)]) + + np.random.seed(22*self.information['simnumber']) + noise_r = np.random.normal(loc=0.0, scale=self.information['r_rdnois'+str( + k)], size=(4616+self.overscan, 1152+self.overscan+self.prescan)) + self.image_r[xmin[k]:xmax[k], ymin[k]:ymax[k]] += noise_r + self.log.info('Bias of %i counts were added to the r band image' % + self.information['r_rdnois'+str(k)]) + + np.random.seed(222*self.information['simnumber']) + noise_i = np.random.normal(loc=0.0, scale=self.information['i_rdnois'+str( + k)], size=(4616+self.overscan, 1152+self.overscan+self.prescan)) + self.image_i[xmin[k]:xmax[k], ymin[k]:ymax[k]] += noise_i + self.log.info('Bias of %i counts were added to the i band image' % + self.information['i_rdnois'+str(k)]) + + # end for + + else: + + np.random.seed(3*self.information['simnumber']) + noise = np.random.normal(loc=0.0, scale=4, size=self.image_g.shape) + self.log.info('Sum of readnoise in g channel = %f' % np.sum(noise)) + # add to the image + self.image_g += noise + + ####################################################### + np.random.seed(33*self.information['simnumber']) + noise = np.random.normal(loc=0.0, scale=4, size=self.image_r.shape) + self.log.info('Sum of readnoise in r channel= %f' % np.sum(noise)) + # add to the image + self.image_r += noise + + ######################################## + np.random.seed(333*self.information['simnumber']) + noise = np.random.normal(loc=0.0, scale=4, size=self.image_i.shape) + self.log.info('Sum of readnoise in i channel= %f' % np.sum(noise)) + # add to the image + self.image_i += noise + + return + +########################################################################################################## + + def electrons2ADU(self): + """ + Convert from electrons to ADUs using the value read from the configuration file. + """ + ##### + if self.overscans: + ###################### + xmin = dict() + xmax = dict() + ymin = dict() + ymax = dict() + + # calculate the non-over subimg of the final matrix of the image + for k in range(1, 17): + xmin[k] = 0 + xmax[k] = 4616+self.overscan + + ymin[k] = (1152+self.overscan+self.prescan)*(k-1) + ymax[k] = (1152+self.overscan+self.prescan)*(k) + ############### + + self.image_g[xmin[k]:xmax[k], ymin[k]:ymax[k] + ] /= self.information['g_gain'+str(k)] + + self.image_r[xmin[k]:xmax[k], ymin[k]:ymax[k] + ] /= self.information['r_gain'+str(k)] + + self.image_i[xmin[k]:xmax[k], ymin[k]:ymax[k] + ] /= self.information['i_gain'+str(k)] + + # end for + return + + #################################################################################### + + def applyBias(self): + """ + Adds a bias level to the image being constructed. + + The value of bias is read from the configure file and stored + in the information dictionary (key bias). + + """ + + xmin = dict() + xmax = dict() + ymin = dict() + ymax = dict() + + # calculate the non-over subimg of the final matrix of the image + for k in range(1, 17): + xmin[k] = 0 + xmax[k] = 4616+self.overscan + + ymin[k] = (1152+self.prescan+self.overscan)*(k-1) + ymax[k] = (1152+self.prescan+self.overscan)*(k) + ############### + self.image_g[xmin[k]:xmax[k], ymin[k]:ymax[k] + ] += self.information['g_detbia'+str(k)] + self.log.info('Bias of %i counts were added to the g band image' % + self.information['g_detbia'+str(k)]) + + self.image_r[xmin[k]:xmax[k], ymin[k]:ymax[k] + ] += self.information['r_detbia'+str(k)] + self.log.info('Bias of %i counts were added to the g band image' % + self.information['r_detbia'+str(k)]) + + self.image_i[xmin[k]:xmax[k], ymin[k]:ymax[k] + ] += self.information['i_detbia'+str(k)] + self.log.info('Bias of %i counts were added to the g band image' % + self.information['i_detbia'+str(k)]) + + # end for + return + +############################################################################## + def img_fits_save(self, data, filename): + #### + + ofd_g = fits.PrimaryHDU() + #### World coordinate system and related parameters ##### + hdu_g = fits.ImageHDU(data) + + hdu_g.header['WCSAXES'] = (np.int16(2), 'number of WCS axes') + hdu_g.header['CRPIX1'] = ( + round(float(self.information['CRPIX1']), 1), 'x-coordinate of reference pixel') + hdu_g.header['CRPIX2'] = ( + round(float(self.information['CRPIX2']), 1), 'y-coordinate of reference pixel') + + hdu_g.header['CRVAL1'] = ( + float(self.information['CRVAL1']), 'first axis value at reference pixel') + hdu_g.header['CRVAL2'] = ( + float(self.information['CRVAL2']), 'second axis value at reference pixel') + hdu_g.header['CTYPE1'] = ( + 'RA---TAN', 'the coordinate type for the first axis') + hdu_g.header['CTYPE2'] = ( + 'DEC--TAN', 'the coordinate type for the second axis') + hdu_g.header['CD1_1'] = ( + float(self.information['CD1_1']), 'partial of first axis coordinate w.r.t. x') + hdu_g.header['CD1_2'] = ( + float(self.information['CD1_2']), 'partial of first axis coordinate w.r.t. y') + hdu_g.header['CD2_1'] = ( + float(self.information['CD2_1']), 'partial of second axis coordinate w.r.t. x') + hdu_g.header['CD2_2'] = ( + float(self.information['CD2_2']), 'partial of second axis coordinate w.r.t. y') + + hdulist_g = fits.HDUList([ofd_g, hdu_g]) + + file_g = self.result_path+'/ori_Sky/'+filename + '.fits' + + hdulist_g.writeto(file_g, overwrite=True) + return + + +#################################################################################### + + + def applyBleeding(self, img, direction='not_horizon'): + """ + Apply bleeding along the CCD readout direction if the number of electrons in a pixel exceeds the full-well capacity. + + direction: this is the CCD readout direction, and default value is horizon direction + + :return: the bleeding image + """ + + data = img.copy() + + if data.max() < self.information['fullwellcapacity']: + return img + + else: + + if direction == 'horizon': + + # loop over each column, as bleeding is modelled column-wise + for i, column in enumerate(data): # select one solumnn + if column.max() <= self.information['fullwellcapacity']: + continue + sum = 0. + for j, value in enumerate(column): + # first round - from bottom to top (need to half the bleeding) + overload = value - self.information['fullwellcapacity'] + if overload > 0.: + overload /= 2. + #self.image[j, i] -= overload + data[i, j] -= overload + sum += overload + + elif sum > 0.: + if -overload > sum: + overload = -sum + #self.image[j, i] -= overload + data[i, j] -= overload + sum += overload + ###################################################### + for i, column in enumerate(data): + if column.max() <= self.information['fullwellcapacity']: + continue + sum = 0. + for j, value in enumerate(column[::-1]): + # second round - from top to bottom (bleeding was half'd already, so now full) + overload = value - self.information['fullwellcapacity'] + if overload > 0.: + #self.image[-j-1, i] -= overload + data[i, -j-1] -= overload + sum += overload + + elif sum > 0.: + if -overload > sum: + overload = -sum + #self.image[-j-1, i] -= overload + data[i, -j-1,] -= overload + sum += overload + + else: + + # loop over each column, as bleeding is modelled column-wise + for i, column in enumerate(data.T): + sum = 0. + for j, value in enumerate(column): + # first round - from bottom to top (need to half the bleeding) + overload = value - self.information['fullwellcapacity'] + if overload > 0.: + overload /= 2. + #self.image[j, i] -= overload + data[j, i] -= overload + + sum += overload + elif sum > 0.: + if -overload > sum: + overload = -sum + #self.image[j, i] -= overload + data[j, i] -= overload + sum += overload + ################################ + for i, column in enumerate(data.T): + sum = 0. + for j, value in enumerate(column[::-1]): + # second round - from top to bottom (bleeding was half'd already, so now full) + overload = value - self.information['fullwellcapacity'] + if overload > 0.: + #self.image[-j-1, i] -= overload + data[-j-1, i] -= overload + + sum += overload + elif sum > 0.: + if -overload > sum: + overload = -sum + #self.image[-j-1, i] -= overload + data[-j-1, i] -= overload + sum += overload + + #####print('Applying column bleeding finished.......') + return data + + ############################################################################ + + def discretise(self, max=2**16-1): + """ + Converts a floating point image array (self.image) to an integer array with max values + defined by the argument max. + + :param max: maximum value the the integer array may contain [default 65k] + :type max: float + + :return: None + """ + + # avoid negative numbers in case bias level was not added + self.image_g[self.image_g < 0.0] = 0. + # cut of the values larger than max + self.image_g[self.image_g > max] = max + + self.image_g = np.rint(self.image_g).astype(int) + self.log.info('Maximum and total values of the image are %i and %i, respectively' % (np.max(self.image_g), + np.sum(self.image_g))) + + # avoid negative numbers in case bias level was not added + self.image_r[self.image_r < 0.0] = 0. + # cut of the values larger than max + self.image_r[self.image_r > max] = max + + self.image_r = np.rint(self.image_r).astype(int) + self.log.info('Maximum and total values of the image are %i and %i, respectively' % (np.max(self.image_r), + np.sum(self.image_r))) + + ############################################################################### + + # avoid negative numbers in case bias level was not added + self.image_i[self.image_i < 0.0] = 0. + # cut of the values larger than max + self.image_i[self.image_i > max] = max + + self.image_i = np.rint(self.image_i).astype(int) + self.log.info('Maximum and total values of the image are %i and %i, respectively' % (np.max(self.image_i), + np.sum(self.image_i))) + + ################################################################################################## + def addScanEffect(self, img): + # read CCD image to matrix as defined data format + ''' + function: add prescan and overscan empty zone to the CCD physical image + + return: img + + ''' + + ### xsize is column, ysize is row + + # horizon direction , columnn + self.prescan = self.information['prescan'] + self.overscan = self.information['overscan'] # vertical direction, row + + row = self.information['ysize'] + col = self.information['xsize'] + if row != 9232 or col != 9216: + print('image size is not correct.......') + sys.exit(1) + #### + ########################################### + xmin = dict() + xmax = dict() + ymin = dict() + ymax = dict() + # physical CCD zone 1-16 + for n in range(1, 17): + + # physical zone index + if n < 9: + xmin[n] = 4616 # row number + xmax[n] = 4616*2 + + ymin[n] = (n-1)*1152 # col number + ymax[n] = ymin[n]+1152 + else: + xmin[n] = 0 + xmax[n] = 4616 + + ymin[n] = (16-n)*1152 + ymax[n] = ymin[n]+1152 + ########################################### + + temp = np.zeros( + (int(4616+self.overscan), int((1152+self.overscan+self.prescan)*16))) + + matxmin = dict() + matxmax = dict() + matymin = dict() + matymax = dict() + +# get the non-over subimg of the final matrix of the image + for k in range(1, 17): + + # big new image include prescan and overscan; + matxmin[k] = 0 # row number + matxmax[k] = 4616 # row number + + matymin[k] = (1152+self.overscan+self.prescan) * \ + (k-1)+self.prescan # col number start + matymax[k] = matymin[k]+1152 # col number end + + ################################ + + for k in range(1, 17): + if k <= 4: # zos1 zos2 zos3 zos4 + # get k-th subimg from input img + subimg = img[xmin[k]:xmax[k], ymin[k]:ymax[k]] + subimg = np.flipud(subimg) # down to up, up to down + + if k >= 5 and k <= 8: # zos5 zos6 zos7 zos8 + # get k-th subimg from input img + subimg = img[xmin[k]:xmax[k], ymin[k]:ymax[k]] + subimg = np.flipud(subimg) # down to up, up to down + subimg = np.fliplr(subimg) # left to right ,right to left + + if k >= 9 and k <= 12: # zos5 zos6 zos7 zos8 + # get k-th subimg from input img + subimg = img[xmin[k]:xmax[k], ymin[k]:ymax[k]] + subimg = np.fliplr(subimg) # left to right ,right to left + + if k >= 13 and k <= 16: # zos5 zos6 zos7 zos8 + # get k-th subimg from input img + subimg = img[xmin[k]:xmax[k], ymin[k]:ymax[k]] + + ###################### + + temp[matxmin[k]:matxmax[k], matymin[k]:matymax[k]] = subimg[:, :] + # end for + return temp + ########################################################################### + + ######################################################################################### + + def writeOutputs(self, obnum=1): + """ + Writes out a FITS file using PyFITS and converts the image array to 16bit unsigned integer as + appropriate for VIS. + + Updates header with the input values and flags used during simulation. + """ + + # Readout information + + HeaderTest = 'no' + ########################################################################## + if self.source == 'EXDF': + obsid = 20100000000+obnum + + if self.source == 'STAR': + obsid = 20200000000+obnum + + if self.source == 'PIPR': + obsid = 20300000000+obnum + + if self.source == 'TRNS': + obsid = 20400000000+obnum + + if self.source == 'COMB': + obsid = 20500000000+obnum + + if self.source == 'CALI': + obsid = 20600000000+obnum + + if self.source == 'BIAS': + obsid = 20700000000+obnum + + if self.source == 'DARK': + obsid = 20800000000+obnum + + if self.source == 'FLAT': + obsid = 20900000000+obnum + ############################# + + if self.source == 'OSBS': + obsid = 21000000000+obnum + + if self.source == 'OSDK': + obsid = 21100000000+obnum + + if self.source == 'OSFT': + obsid = 21200000000+obnum + + if self.source == 'FRNG': + obsid = 21300000000+obnum + + ####################################################################### + exp_starttime = self.dt.strftime("%Y%m%d%H%M%S") + + # exposure end time is t2 ; + t2 = self.dt+timedelta(seconds=self.information['exptime']) + exp_endtime = t2.strftime("%Y%m%d%H%M%S") + ### data read time is the exposure end time ### + # data read end time is t3, and this is the write file time; + t3 = self.dt+timedelta(seconds=self.information['exptime'])+timedelta( + seconds=self.information['readouttime']) + + if self.filterP[self.filter_g]['number'] > 9: + filternum = str(self.filterP[self.filter_g]['number']) + else: + filternum = '0'+str(self.filterP[self.filter_g]['number']) + + filename_g = 'CSST_MCI_C1_'+self.source+'_'+exp_starttime+'_'+exp_endtime+'_'+str(obsid)+'_' + \ + filternum+'_L0_V01' + + ####################################################################### + ####################################################################### + # create a new FITS file, using HDUList instance + ofd_g = fits.PrimaryHDU() + # add primary header + ofd_g = fits.PrimaryHDU() + ofd_g.header['SIMPLE'] = (True, 'conforms to FITS standard') + ofd_g.header['BITPIX'] = (0, 'array data type') + ofd_g.header['NAXIS'] = (0, 'number of array dimensions') + ofd_g.header['NEXTEND'] = (np.int16(1), 'number of array dimensions') + ####################################################################### + ########################### + temp = t3.utcnow() + data_time = temp.strftime("%Y-%m-%dT%H:%M:%S.%f") + ofd_g.header['DATE'] = ( + data_time[:21], 'written date (yyyy-mm-ddThh:mm:ss.s)') + ofd_g.header['FILENAME'] = (filename_g[:68], '') + ofd_g.header['FILETYPE'] = (self.source[:12], 'observation type') + ofd_g.header['TELESCOP'] = ('CSST', 'telescope name') + ofd_g.header['INSTRUME'] = ('MCI', 'instrument name') + ofd_g.header['RADECSYS'] = ( + 'ICRS', 'coordinate system of the object') + ofd_g.header['EQUINOX'] = (float(2000.0), '') + ofd_g.header['FITSSWV'] = ( + 'mci_sim_0.8.03', 'FITS creating software version') + ############################################################################## + ######### Object information ######################################### + ofd_g.header['OBJECT'] = ( + self.information['name_obj'][:30], 'object name') + ofd_g.header['TARGET'] = ( + self.information['target'][:15], 'target name (hhmmss.s+ddmmss)') + ####################################################################### + ofd_g.header['OBSID'] = (str(obsid), 'observation ID') + ofd_g.header['RA_OBJ'] = ( + float(self.information['ra_obj']), 'object RA (deg)') + ofd_g.header['DEC_OBJ'] = ( + float(self.information['dec_obj']), 'object Dec (deg)') + + ######## Telescope information ############### + + ofd_g.header['REFFRAME'] = ('CSSTGSC-1.0', 'guiding catalog version') + ofd_g.header['DATE-OBS'] = (data_time[:21], + 'observation date (yyyy-mm-ddThh:mm:ss.s)') + ofd_g.header['SATESWV'] = ( + 'softwave-1.0', 'satellite software version') + ofd_g.header['EXPSTART'] = (np.float64( + time2mjd(self.dt)), 'exposure start time (MJD)') + ofd_g.header['CABSTART'] = (np.float64( + time2jd(self.dt)), 'first cabin time after exposure start (MJD)') + ofd_g.header['SUNANGL0'] = (np.float32( + 0.0), 'angle between the Sun and opt axis at CABSTART') + ofd_g.header['MOONANG0'] = (np.float32( + 0.0), 'angle between the Moon and opt axis at CABSTART') + ofd_g.header['TEL_ALT0'] = (np.float64( + 0.0), 'angle between opt axis and Elimb at CABSTART') + ofd_g.header['POS_ANG0'] = (np.float64( + 0.0), 'angle between y axis and North Pole at CABSTART') + ofd_g.header['POSI0_X'] = (np.float64( + self.information['pos_x']), 'orbital position in X at CABSTART (km)') + ofd_g.header['POSI0_Y'] = (np.float64( + self.information['pos_y']), 'orbital position in Y at CABSTART (km)') + ofd_g.header['POSI0_Z'] = (np.float64( + self.information['pos_z']), 'orbital position in Z at CABSTART (km)') + ofd_g.header['VELO0_X'] = (np.float64( + self.information['velocity_x']), 'orbital velocity in X at CABSTART (km/s)') + ofd_g.header['VELO0_Y'] = (np.float64( + self.information['velocity_y']), 'orbital velocity in Y at CABSTART (km/s)') + ofd_g.header['VELO0_Z'] = (np.float64( + self.information['velocity_z']), 'orbital velocity in Z at CABSTART (km/s)') + ofd_g.header['Euler0_1'] = (np.float64( + 0.0), 'Euler angle 1 at CABSTART (deg)') + ofd_g.header['Euler0_2'] = (np.float64( + 0.0), 'Euler angle 2 at CABSTART (deg)') + ofd_g.header['Euler0_3'] = (np.float64( + 0.0), 'Euler angle 3 at CABSTART (deg)') + + ############################################################## + ##### + ofd_g.header['RA_PNT0'] = ( + float(self.information['ra_pnt0']), 'pointing RA at CABSTART (deg)') + ofd_g.header['DEC_PNT0'] = ( + float(self.information['dec_pnt0']), 'pointing Dec at CABSTART (deg)') + ofd_g.header['EXPEND'] = (float(time2mjd(t2)), + 'exposure end time (MJD)') + ofd_g.header['CABEND'] = (float(time2jd(t2)), + 'first cabin time after exposure end (MJD)') + ofd_g.header['SUNANGL1'] = (np.float32( + 0.0), 'angle between the Sun and opt axis at CABEND') + ofd_g.header['MOONANG1'] = (np.float32( + 0.0), 'angle between the Moon and opt axis at CABEND') + ofd_g.header['TEL_ALT1'] = (np.float64( + 0.0), 'angle between opt axis and Elimb at CABEND') + ofd_g.header['POS_ANG1'] = ( + float(0.0), 'angle between y axis and North Pole at CABEND') + + ############################################################## + + t2jd = time2jd(t2) + + if self.orbit_pars[-1, 0] < t2jd: # orbit parameters are not in currenct txt file + self.orbit_file_num = self.orbit_file_num+1 + fn = self.information['dir_path'] + \ + 'MCI_inputData/TianCe/orbit20160925/' + \ + str(self.orbit_file_num)+'.txt' + self.orbit_pars = np.loadtxt(fn) + self.orbit_exp_num = 0 + + for k in range(self.orbit_exp_num, len(self.orbit_pars), 1): + + if t2jd-self.orbit_pars[k, 0] < 0: + break + if k == 0: + deltaT = jd2time(self.orbit_pars[k, 0])-t2 + p1x = self.orbit_pars[k, 1]-(self.orbit_pars[k+1, 1] - + self.orbit_pars[k, 1])*deltaT.seconds/120 + p1y = self.orbit_pars[k, 2]-(self.orbit_pars[k+1, 2] - + self.orbit_pars[k, 2])*deltaT.seconds/120 + p1z = self.orbit_pars[k, 3]-(self.orbit_pars[k+1, 3] - + self.orbit_pars[k, 3])*deltaT.seconds/120 + + p1vx = self.orbit_pars[k, 4]-(self.orbit_pars[k+1, 4] - + self.orbit_pars[k, 4])*deltaT.seconds/120 + p1vx = self.orbit_pars[k, 5]-(self.orbit_pars[k+1, 5] - + self.orbit_pars[k, 5])*deltaT.seconds/120 + p1vx = self.orbit_pars[k, 6]-(self.orbit_pars[k+1, 6] - + self.orbit_pars[k, 6])*deltaT.seconds/120 + + else: + deltaT = jd2time(self.orbit_pars[k, 0])-t2 + p1x = self.orbit_pars[k-1, 1]+(self.orbit_pars[k, 1] - + self.orbit_pars[k-1, 1])*deltaT.seconds/120 + p1y = self.orbit_pars[k-1, 2]+(self.orbit_pars[k, 2] - + self.orbit_pars[k-1, 2])*deltaT.seconds/120 + p1z = self.orbit_pars[k-1, 3]+(self.orbit_pars[k, 3] - + self.orbit_pars[k-1, 3])*deltaT.seconds/120 + + p1vx = self.orbit_pars[k-1, 4]+(self.orbit_pars[k, 4] - + self.orbit_pars[k-1, 4])*deltaT.seconds/120 + p1vy = self.orbit_pars[k-1, 5]+(self.orbit_pars[k, 5] - + self.orbit_pars[k-1, 5])*deltaT.seconds/120 + p1vz = self.orbit_pars[k-1, 6]+(self.orbit_pars[k, 6] - + self.orbit_pars[k-1, 6])*deltaT.seconds/120 + + ofd_g.header['POSI1_X'] = ( + float(p1x), 'orbital position in X at CABEND (km)') + ofd_g.header['POSI1_Y'] = ( + float(p1y), 'orbital position in Y at CABEND (km)') + ofd_g.header['POSI1_Z'] = ( + float(p1z), 'orbital position in Z at CABEND (km)') + + ofd_g.header['VELO1_X'] = ( + float(p1vx), 'orbital velocity in X at CABEND (km/s)') + ofd_g.header['VELO1_Y'] = ( + float(p1vy), 'orbital velocity in Y at CABEND (km/s)') + ofd_g.header['VELO1_Z'] = ( + float(p1vz), 'orbital velocity in Z at CABEND (km/s)') + + ################### + + ofd_g.header['Euler1_1'] = (np.float64( + 0.0), 'Euler angle 1 at CABEND (deg)') + ofd_g.header['Euler1_2'] = (np.float64( + 0.0), 'Euler angle 2 at CABEND (deg)') + ofd_g.header['Euler1_3'] = (np.float64( + 0.0), 'Euler angle 3 at CABEND (deg)') + ofd_g.header['RA_PNT1'] = ( + float(ofd_g.header['RA_PNT0']), 'pointing RA at CABEND (deg)') + ofd_g.header['DEC_PNT1'] = ( + float(ofd_g.header['DEC_PNT0']), 'pointing Dec at CABEND (deg)') + ofd_g.header['EPOCH'] = (np.float32( + time2jd(t3.utcnow())), 'equinox of pointing RA and Dec') + ofd_g.header['EXPTIME'] = (np.float32( + self.information['exptime']), 'exposure time (s)') + + ofd_g.header['CHECKSUM'] = (0, '0') + ################################################ + ########## finish header for 0 layer #################### + + # new image HDU + + if HeaderTest == 'yes': + hdu_g = fits.ImageHDU() + else: + hdu_g = fits.ImageHDU(data=np.uint16(self.image_g)) + # + ### MCI header for image layer ### + # need review thie header + hdu_g.header['XTENSION'] = ('IMAGE', 'extension type') + hdu_g.header['BITPIX'] = (np.int16(16), 'array data type') + hdu_g.header['NAXIS'] = (np.int16(2), 'number of array dimensions') + + hdu_g.header['NAXIS1'] = (np.int32(4936), '') + hdu_g.header['NAXIS2'] = (np.int32(23984), '') + + hdu_g.header['PCOUNT'] = (np.int16(0), 'number of parameters') + hdu_g.header['GCOUNT'] = (np.int16(1), 'number of groups') + + hdu_g.header['BSCALE'] = (int(1), '') + hdu_g.header['BZERO'] = (int(32768), '') + + hdu_g.header['EXTNAME'] = ('SCI', '') + hdu_g.header['EXTVER'] = (np.int16(1), '') + + hdu_g.header['BUNIT'] = ('ADU', 'physical unit of array values') + + ###### file information ############### + + temp = t3.utcnow() + data_time = temp.strftime("%Y-%m-%dT%H:%M:%S.%f") + hdu_g.header['DATE'] = ( + data_time[:21], 'written date (yyyy-mm-ddThh:mm:ss.s)') + hdu_g.header['FILENAME'] = (filename_g[:68], '') + hdu_g.header['FILETYPE'] = (self.source[:12], 'observation type') + hdu_g.header['RADECSYS'] = ( + 'ICRS', 'coordinate system of the object') + hdu_g.header['EQUINOX'] = (float(2000.1), '') + hdu_g.header['FITSSWV'] = ('4.2.1', 'FITS creating software version') + + ######################################## + + #### instrument status and object information ###### + hdu_g.header['TELESCOP'] = ('CSST', 'telescope name') + hdu_g.header['INSTRUME'] = ('MCI', 'instrument name') + hdu_g.header['CHANNEL'] = ('C1', 'channel number') + hdu_g.header['FILTERNO'] = (filternum, 'filter number') + hdu_g.header['DIFFUSER'] = ( + bool(True), 'insert diffuser status for flat calibration') + hdu_g.header['FLAMP'] = (np.int16(0), 'status of flat lamp') + hdu_g.header['MCISTAT'] = ( + np.int16(0), 'MCI components status parameter') + + hdu_g.header['DATE-OBS'] = (data_time[:21], + 'observation date (yyyy-mm-ddThh:mm:ss.s)') + hdu_g.header['OBJECT'] = ('MCI_obj', 'object name') + + hdu_g.header['TARGET'] = ( + self.information['target'][:15], 'target name (hhmmss.s+ddmmss)') + ####################################################################### + hdu_g.header['OBSID'] = (str(obsid), 'observation ID') + hdu_g.header['RA_OBJ'] = ( + float(self.information['ra_obj']), 'object RA (deg)') + hdu_g.header['DEC_OBJ'] = ( + float(self.information['dec_obj']), 'object Dec (deg)') + + ##### detector and Filter information ##### + hdu_g.header['FILTER'] = (self.filter_g[:6], 'filter band') + ##### Detector information #### + hdu_g.header['DETSN'] = ( + 'E2V-CCD-290-0000000', 'detector serial number') + hdu_g.header['DETNAME'] = ('uv', 'detector name') + # need review here + + hdu_g.header['DETTEMP0'] = (np.float32( + 0), 'detector temperature at EXPSTART (K)') + hdu_g.header['DETTEMP1'] = (np.float32( + 0), 'detector temperature at EXPEND (K)') + hdu_g.header['DETTEMP2'] = (np.float32( + 0), 'detector temperature at READT1 (K)') + + hdu_g.header['DETSIZE'] = (str( + self.information['ysize'])+'*'+str(self.information['xsize']), 'detector size') + hdu_g.header['DATASECT'] = ( + str(self.image_g.shape[0])+'*'+str(self.image_g.shape[1]), 'data section') + + hdu_g.header['PIXSCAL1'] = ( + float(0.05), 'pixel scale for axis 1 (arcsec/pixel)') + hdu_g.header['PIXSCAL2'] = ( + float(0.05), 'pixel scale for axis 2 (arcsec/pixel)') + + hdu_g.header['PIXSIZE1'] = (int(10), 'pixel size for axis 1 (micron)') + hdu_g.header['PIXSIZE2'] = (int(10), 'pixel size for axis 2 (micron)') + + hdu_g.header['NCHANNEL'] = (np.int16(16), 'number of readout channels') + + hdu_g.header['PSCAN1'] = ( + np.int32(27), 'horizontal prescan width, per readout channel') + hdu_g.header['PSCAN2'] = ( + np.int32(0), 'vertical prescan width, per readout channel') + + hdu_g.header['OSCAN1'] = ( + np.int32(320), 'horizontal overscan width, per readout channel') + hdu_g.header['OSCAN2'] = ( + np.int32(320), 'vertical overscan width, per readout channel') + + hdu_g.header['BIN_X'] = (np.int16(1), 'bin number in X') + hdu_g.header['BIN_Y'] = (np.int16(1), 'bin number in Y') + + #### World coordinate system information ##### + + hdu_g.header['WCSAXES'] = (np.int16(2), 'number of WCS axes') + hdu_g.header['CRPIX1'] = ( + round(float(self.information['CRPIX1']), 1), 'x-coordinate of reference pixel') + hdu_g.header['CRPIX2'] = ( + round(float(self.information['CRPIX2']), 1), 'y-coordinate of reference pixel') + + hdu_g.header['CRVAL1'] = ( + float(self.information['CRVAL1']), 'value of reference pixel on axis 1') + hdu_g.header['CRVAL2'] = ( + float(self.information['CRVAL2']), 'value of reference pixel on axis 2') + + hdu_g.header['CTYPE1'] = ('RA---TAN', 'type of RA WCS projection') + hdu_g.header['CTYPE2'] = ('DEC--TAN', 'type of Dec WCS projection') + + hdu_g.header['CD1_1'] = ( + float(self.information['CD1_1']), 'transformation matrix element (deg/pix)') + hdu_g.header['CD1_2'] = ( + float(self.information['CD1_2']), 'transformation matrix element (deg/pix)') + hdu_g.header['CD2_1'] = ( + float(self.information['CD2_1']), 'transformation matrix element (deg/pix)') + hdu_g.header['CD2_2'] = ( + float(self.information['CD2_2']), 'transformation matrix element (deg/pix)') + + ###### Readout information ####### + hdu_g.header['GAINLVL'] = (np.float32( + 1.5), 'gain level (e-/ADU)') + hdu_g.header['GAIN01'] = (np.float32( + self.information['g_gain1']), 'CCD gain [channel 1] (e-/ADU)') + hdu_g.header['GAIN02'] = (np.float32( + self.information['g_gain2']), 'CCD gain [channel 2] (e-/ADU)') + hdu_g.header['GAIN03'] = (np.float32( + self.information['g_gain3']), 'CCD gain [channel 3] (e-/ADU)') + hdu_g.header['GAIN04'] = (np.float32( + self.information['g_gain4']), 'CCD gain [channel 4] (e-/ADU)') + hdu_g.header['GAIN05'] = (np.float32( + self.information['g_gain5']), 'CCD gain [channel 5] (e-/ADU)') + hdu_g.header['GAIN06'] = (np.float32( + self.information['g_gain6']), 'CCD gain [channel 6] (e-/ADU)') + hdu_g.header['GAIN07'] = (np.float32( + self.information['g_gain7']), 'CCD gain [channel 7] (e-/ADU)') + hdu_g.header['GAIN08'] = (np.float32( + self.information['g_gain8']), 'CCD gain [channel 8] (e-/ADU)') + hdu_g.header['GAIN09'] = (np.float32( + self.information['g_gain9']), 'CCD gain [channel 9] (e-/ADU)') + hdu_g.header['GAIN10'] = (np.float32( + self.information['g_gain10']), 'CCD gain [channel 10] (e-/ADU)') + hdu_g.header['GAIN11'] = (np.float32( + self.information['g_gain11']), 'CCD gain [channel 11] (e-/ADU)') + hdu_g.header['GAIN12'] = (np.float32( + self.information['g_gain12']), 'CCD gain [channel 12] (e-/ADU)') + hdu_g.header['GAIN13'] = (np.float32( + self.information['g_gain13']), 'CCD gain [channel 13] (e-/ADU)') + hdu_g.header['GAIN14'] = (np.float32( + self.information['g_gain14']), 'CCD gain [channel 14] (e-/ADU)') + hdu_g.header['GAIN15'] = (np.float32( + self.information['g_gain15']), 'CCD gain [channel 15] (e-/ADU)') + hdu_g.header['GAIN16'] = (np.float32( + self.information['g_gain16']), 'CCD gain [channel 16] (e-/ADU)') + ####### + hdu_g.header['RON01'] = (np.float32( + self.information['g_rdnois1']), 'readout noise [channel 01] (e-)') + hdu_g.header['RON02'] = (np.float32( + self.information['g_rdnois2']), 'readout noise [channel 02] (e-)') + hdu_g.header['RON03'] = (np.float32( + self.information['g_rdnois3']), 'readout noise [channel 03] (e-)') + hdu_g.header['RON04'] = (np.float32( + self.information['g_rdnois4']), 'readout noise [channel 04] (e-)') + hdu_g.header['RON05'] = (np.float32( + self.information['g_rdnois5']), 'readout noise [channel 05] (e-)') + hdu_g.header['RON06'] = (np.float32( + self.information['g_rdnois6']), 'readout noise [channel 06] (e-)') + hdu_g.header['RON07'] = (np.float32( + self.information['g_rdnois7']), 'readout noise [channel 07] (e-)') + hdu_g.header['RON08'] = (np.float32( + self.information['g_rdnois8']), 'readout noise [channel 08] (e-)') + hdu_g.header['RON09'] = (np.float32( + self.information['g_rdnois9']), 'readout noise [channel 09] (e-)') + hdu_g.header['RON10'] = (np.float32( + self.information['g_rdnois10']), 'readout noise [channel 10] (e-)') + hdu_g.header['RON11'] = (np.float32( + self.information['g_rdnois11']), 'readout noise [channel 11] (e-)') + hdu_g.header['RON12'] = (np.float32( + self.information['g_rdnois12']), 'readout noise [channel 12] (e-)') + hdu_g.header['RON13'] = (np.float32( + self.information['g_rdnois13']), 'readout noise [channel 13] (e-)') + hdu_g.header['RON14'] = (np.float32( + self.information['g_rdnois14']), 'readout noise [channel 14] (e-)') + hdu_g.header['RON15'] = (np.float32( + self.information['g_rdnois15']), 'readout noise [channel 15] (e-)') + hdu_g.header['RON16'] = (np.float32( + self.information['g_rdnois16']), 'readout noise [channel 16] (e-)') + + ####### + hdu_g.header['DETBIA01'] = (np.float32( + self.information['g_detbia1']), 'amplifier bias grey value [channel 1] (ADU)') + hdu_g.header['DETBIA02'] = (np.float32( + self.information['g_detbia2']), 'amplifier bias grey value [channel 2] (ADU)') + hdu_g.header['DETBIA03'] = (np.float32( + self.information['g_detbia3']), 'amplifier bias grey value [channel 3] (ADU)') + hdu_g.header['DETBIA04'] = (np.float32( + self.information['g_detbia4']), 'amplifier bias grey value [channel 4] (ADU)') + hdu_g.header['DETBIA05'] = (np.float32( + self.information['g_detbia5']), 'amplifier bias grey value [channel 5] (ADU)') + hdu_g.header['DETBIA06'] = (np.float32( + self.information['g_detbia6']), 'amplifier bias grey value [channel 6] (ADU)') + hdu_g.header['DETBIA07'] = (np.float32( + self.information['g_detbia7']), 'amplifier bias grey value [channel 7] (ADU)') + hdu_g.header['DETBIA08'] = (np.float32( + self.information['g_detbia8']), 'amplifier bias grey value [channel 8] (ADU)') + hdu_g.header['DETBIA09'] = (np.float32( + self.information['g_detbia9']), 'amplifier bias grey value [channel 9] (ADU)') + hdu_g.header['DETBIA10'] = (np.float32( + self.information['g_detbia10']), 'amplifier bias grey value [channel 10] (ADU)') + hdu_g.header['DETBIA11'] = (np.float32( + self.information['g_detbia11']), 'amplifier bias grey value [channel 11] (ADU)') + hdu_g.header['DETBIA12'] = (np.float32( + self.information['g_detbia12']), 'amplifier bias grey value [channel 12] (ADU)') + hdu_g.header['DETBIA13'] = (np.float32( + self.information['g_detbia13']), 'amplifier bias grey value [channel 13] (ADU)') + hdu_g.header['DETBIA14'] = (np.float32( + self.information['g_detbia14']), 'amplifier bias grey value [channel 14] (ADU)') + hdu_g.header['DETBIA15'] = (np.float32( + self.information['g_detbia15']), 'amplifier bias grey value [channel 15] (ADU)') + hdu_g.header['DETBIA16'] = (np.float32( + self.information['g_detbia16']), 'amplifier bias grey value [channel 16] (ADU)') + hdu_g.header['ROSPEED'] = (np.float32(100), 'readout speed (MHz)') + + ###################################### + #### exposure and shutter information ##### + hdu_g.header['EXPSTART'] = (np.float64( + time2mjd(self.dt)), 'exposure start time (MJD)') + hdu_g.header['EXPEND'] = (float(time2mjd(t2)), + 'exposure end time (MJD)') + hdu_g.header['EXPTIME'] = (np.float32( + self.information['exptime']), 'exposure time (s)') + hdu_g.header['DARKTIME'] = (np.float32( + self.information['exptime']), 'dark current time (s)') + hdu_g.header['SHTSTAT'] = (bool(True), 'shutter status') + + ###### satellite and its allitide information ############## + hdu_g.header['REFFRAME'] = ('CSSTGSC-1.0', 'guiding catalog version') + hdu_g.header['SATESWV'] = ( + 'softwave-1.0', 'satellite software version') + hdu_g.header['CABSTART'] = (np.float64( + time2jd(self.dt)), 'first cabin time after exposure start (MJD)') + hdu_g.header['SUNANGL0'] = (np.float32( + 0.0), 'angle between the Sun and opt axis at CABSTART') + hdu_g.header['MOONANG0'] = (np.float32( + 0.0), 'angle between the Moon and opt axis at CABSTART') + hdu_g.header['TEL_ALT0'] = (np.float64( + 0.0), 'angle between opt axis and Elimb at CABSTART') + hdu_g.header['POS_ANG0'] = (np.float64( + 0.0), 'angle between y axis and North Pole at CABSTART') + hdu_g.header['POSI0_X'] = (np.float64( + self.information['pos_x']), 'orbital position in X at CABSTART (km)') + hdu_g.header['POSI0_Y'] = (np.float64( + self.information['pos_y']), 'orbital position in Y at CABSTART (km)') + hdu_g.header['POSI0_Z'] = (np.float64( + self.information['pos_z']), 'orbital position in Z at CABSTART (km)') + hdu_g.header['VELO0_X'] = (np.float64( + self.information['velocity_x']), 'orbital velocity in X at CABSTART (km/s)') + hdu_g.header['VELO0_Y'] = (np.float64( + self.information['velocity_y']), 'orbital velocity in Y at CABSTART (km/s)') + hdu_g.header['VELO0_Z'] = (np.float64( + self.information['velocity_z']), 'orbital velocity in Z at CABSTART (km/s)') + hdu_g.header['Euler0_1'] = (np.float64( + 0.0), 'Euler angle 1 at CABSTART (deg)') + hdu_g.header['Euler0_2'] = (np.float64( + 0.0), 'Euler angle 2 at CABSTART (deg)') + hdu_g.header['Euler0_3'] = (np.float64( + 0.0), 'Euler angle 3 at CABSTART (deg)') + hdu_g.header['RA_PNT0'] = ( + float(self.information['ra_pnt0']), 'pointing RA at CABSTART (deg)') + hdu_g.header['DEC_PNT0'] = ( + float(self.information['dec_pnt0']), 'pointing Dec at CABSTART (deg)') + hdu_g.header['CABEND'] = (float(time2jd(t2)), + 'first cabin time after exposure end (MJD)') + hdu_g.header['SUNANGL1'] = (np.float32( + 0.0), 'angle between the Sun and opt axis at CABEND') + hdu_g.header['MOONANG1'] = (np.float32( + 0.0), 'angle between the Moon and opt axis at CABEND') + hdu_g.header['TEL_ALT1'] = (np.float64( + 0.0), 'angle between opt axis and Elimb at CABEND') + hdu_g.header['POS_ANG1'] = ( + float(0.0), 'angle between y axis and North Pole at CABEND') + hdu_g.header['POSI1_X'] = ( + float(p1x), 'orbital position in X at CABEND (km)') + hdu_g.header['POSI1_Y'] = ( + float(p1y), 'orbital position in Y at CABEND (km)') + hdu_g.header['POSI1_Z'] = ( + float(p1z), 'orbital position in Z at CABEND (km)') + hdu_g.header['VELO1_X'] = ( + float(p1vx), 'orbital velocity in X at CABEND (km/s)') + hdu_g.header['VELO1_Y'] = ( + float(p1vy), 'orbital velocity in Y at CABEND (km/s)') + hdu_g.header['VELO1_Z'] = ( + float(p1vz), 'orbital velocity in Z at CABEND (km/s)') + hdu_g.header['Euler1_1'] = (np.float64( + 0.0), 'Euler angle 1 at CABEND (deg)') + hdu_g.header['Euler1_2'] = (np.float64( + 0.0), 'Euler angle 2 at CABEND (deg)') + hdu_g.header['Euler1_3'] = (np.float64( + 0.0), 'Euler angle 3 at CABEND (deg)') + hdu_g.header['RA_PNT1'] = ( + float(hdu_g.header['RA_PNT0']), 'pointing RA at CABEND (deg)') + hdu_g.header['DEC_PNT1'] = ( + float(hdu_g.header['DEC_PNT0']), 'pointing Dec at CABEND (deg)') + hdu_g.header['EPOCH'] = (np.float32( + time2jd(t3.utcnow())), 'equinox of pointing RA and Dec') + hdu_g.header['CHECKSUM'] = (0, '0') + ############################################################## + + # finish MCI header for image layer in G channel + # 333 + + hdulist_g = fits.HDUList([ofd_g, hdu_g]) + #### + ofd_g.header.add_comment( + '========================================================================', after='FITSSWV') + ofd_g.header.add_comment('OBJECT INFORMATION', after='FITSSWV') + ofd_g.header.add_comment( + '========================================================================', after='FITSSWV') + + ofd_g.header.add_comment( + '========================================================================', after='DEC_OBJ') + ofd_g.header.add_comment('TELESCOPE INFORMATION', after='DEC_OBJ') + ofd_g.header.add_comment( + '========================================================================', after='DEC_OBJ') + + ofd_g.header.add_comment( + '========================================================================', after='EXPTIME') + ofd_g.header.add_comment('VERIFICATION INFORMATION', after='EXPTIME') + ofd_g.header.add_comment( + '========================================================================', after='EXPTIME') + + ####################################################################### + hdu_g.header.add_comment( + '========================================================================', after='BUNIT') + hdu_g.header.add_comment('FILE INFORMATION', after='BUNIT') + hdu_g.header.add_comment( + '========================================================================', after='BUNIT') + + hdu_g.header.add_comment( + '========================================================================', after='FITSSWV') + hdu_g.header.add_comment( + 'INSTRUMENT STATUS AND OBJECT INFORMATION', after='FITSSWV') + hdu_g.header.add_comment( + '========================================================================', after='FITSSWV') + + hdu_g.header.add_comment( + '========================================================================', after='DEC_OBJ') + hdu_g.header.add_comment( + 'DETECTOR AND FILTER INFORMATION', after='DEC_OBJ') + hdu_g.header.add_comment( + '========================================================================', after='DEC_OBJ') + + hdu_g.header.add_comment( + '========================================================================', after='BIN_Y') + hdu_g.header.add_comment( + 'WORLD COORDINATE SYSTEM INFORMATION', after='BIN_Y') + hdu_g.header.add_comment( + '========================================================================', after='BIN_Y') + + hdu_g.header.add_comment( + '========================================================================', after='CD2_2') + hdu_g.header.add_comment('READOUT INFORMATION', after='CD2_2') + hdu_g.header.add_comment( + '========================================================================', after='CD2_2') + + hdu_g.header.add_comment( + '========================================================================', after='ROSPEED') + hdu_g.header.add_comment( + 'EXPOSURE AND SHUTTER INFORMATION', after='ROSPEED') + hdu_g.header.add_comment( + '========================================================================', after='ROSPEED') + + hdu_g.header.add_comment( + '========================================================================', after='SHTSTAT') + hdu_g.header.add_comment( + 'SATELLITE AND ITS ATTITUDE INFORMATION', after='SHTSTAT') + hdu_g.header.add_comment( + '========================================================================', after='SHTSTAT') + + hdu_g.header.add_comment( + '========================================================================', after='EPOCH') + hdu_g.header.add_comment('VERIFICATION INFORMATION', after='EPOCH') + hdu_g.header.add_comment( + '========================================================================', after='EPOCH') + + ##################################### + + if self.source == 'DARK' or self.source == 'FLAT' or self.source == 'BIAS': + file_g = self.result_path+'/cali_Data/'+filename_g+'.fits' + else: + file_g = self.result_path+'/sky_Data/'+filename_g + '.fits' + + hdulist_g.writeto(file_g, checksum=True) + + ########################################################################### + # r band + + # create a new FITS file, using HDUList instance + ofd_r = fits.PrimaryHDU() + + if self.filterP[self.filter_r]['number'] > 9: + filternum = str(self.filterP[self.filter_r]['number']) + else: + filternum = '0'+str(self.filterP[self.filter_r]['number']) + + filename_r = 'CSST_MCI_C2_'+self.source+'_'+exp_starttime+'_'+exp_endtime+'_'+str(obsid)+'_' + \ + filternum+'_L0_V01' + + ####################################################################### + ####################################################################### + # create a new FITS file, using HDUList instance + ofd_r = fits.PrimaryHDU() + # add primary header + + ########## + ofd_r.header['SIMPLE'] = (True, 'conforms to FITS standard') + ofd_r.header['BITPIX'] = (0, 'array data type') + ofd_r.header['NAXIS'] = (0, 'number of array dimensions') + ofd_r.header['NEXTEND'] = (np.int16(1), 'number of array dimensions') + + temp = t3.utcnow() + data_time = temp.strftime("%Y-%m-%dT%H:%M:%S.%f") + ofd_r.header['DATE'] = ( + data_time[:21], 'written date (yyyy-mm-ddThh:mm:ss.s)') + ofd_r.header['FILENAME'] = (filename_r[:68], '') + ofd_r.header['FILETYPE'] = (self.source[:12], 'observation type') + ofd_r.header['TELESCOP'] = ('CSST', 'telescope name') + ofd_r.header['INSTRUME'] = ('MCI', 'instrument name') + ofd_r.header['RADECSYS'] = ( + 'ICRS', 'coordinate system of the object') + ofd_r.header['EQUINOX'] = (float(2000.0), '') + ofd_r.header['FITSSWV'] = ( + 'mci_sim_0.8.03', 'FITS creating software version') + ######### Object information ######################################### + ofd_r.header['OBJECT'] = ( + self.information['name_obj'][:30], 'object name') + ofd_r.header['TARGET'] = ( + self.information['target'][:15], 'target name (hhmmss.s+ddmmss)') + ####################################################################### + ofd_r.header['OBSID'] = (str(obsid), 'observation ID') + ofd_r.header['RA_OBJ'] = ( + float(self.information['ra_obj']), 'object RA (deg)') + ofd_r.header['DEC_OBJ'] = ( + float(self.information['dec_obj']), 'object Dec (deg)') + + ######## Telescope information ############### + + ofd_r.header['REFFRAME'] = ('CSSTGSC-1.0', 'guiding catalog version') + ofd_r.header['DATE-OBS'] = (data_time[:21], + 'observation date (yyyy-mm-ddThh:mm:ss.s)') + ofd_r.header['SATESWV'] = ( + 'softwave-1.0', 'satellite software version') + ofd_r.header['EXPSTART'] = (np.float64( + time2mjd(self.dt)), 'exposure start time (MJD)') + ofd_r.header['CABSTART'] = (np.float64( + time2jd(self.dt)), 'first cabin time after exposure start (MJD)') + ofd_r.header['SUNANGL0'] = (np.float32( + 0.0), 'angle between the Sun and opt axis at CABSTART') + ofd_r.header['MOONANG0'] = (np.float32( + 0.0), 'angle between the Moon and opt axis at CABSTART') + ofd_r.header['TEL_ALT0'] = (np.float64( + 0.0), 'angle between opt axis and Elimb at CABSTART') + ofd_r.header['POS_ANG0'] = (np.float64( + 0.0), 'angle between y axis and North Pole at CABSTART') + ofd_r.header['POSI0_X'] = (np.float64( + self.information['pos_x']), 'orbital position in X at CABSTART (km)') + ofd_r.header['POSI0_Y'] = (np.float64( + self.information['pos_y']), 'orbital position in Y at CABSTART (km)') + ofd_r.header['POSI0_Z'] = (np.float64( + self.information['pos_z']), 'orbital position in Z at CABSTART (km)') + ofd_r.header['VELO0_X'] = (np.float64( + self.information['velocity_x']), 'orbital velocity in X at CABSTART (km/s)') + ofd_r.header['VELO0_Y'] = (np.float64( + self.information['velocity_y']), 'orbital velocity in Y at CABSTART (km/s)') + ofd_r.header['VELO0_Z'] = (np.float64( + self.information['velocity_z']), 'orbital velocity in Z at CABSTART (km/s)') + ofd_r.header['Euler0_1'] = (np.float64( + 0.0), 'Euler angle 1 at CABSTART (deg)') + ofd_r.header['Euler0_2'] = (np.float64( + 0.0), 'Euler angle 2 at CABSTART (deg)') + ofd_r.header['Euler0_3'] = (np.float64( + 0.0), 'Euler angle 3 at CABSTART (deg)') + + ############################################################## + ##### + ofd_r.header['RA_PNT0'] = ( + float(self.information['ra_pnt0']), 'pointing RA at CABSTART (deg)') + ofd_r.header['DEC_PNT0'] = ( + float(self.information['dec_pnt0']), 'pointing Dec at CABSTART (deg)') + ofd_r.header['EXPEND'] = (float(time2mjd(t2)), + 'exposure end time (MJD)') + ofd_r.header['CABEND'] = (float(time2jd(t2)), + 'first cabin time after exposure end (MJD)') + ofd_r.header['SUNANGL1'] = (np.float32( + 0.0), 'angle between the Sun and opt axis at CABEND') + ofd_r.header['MOONANG1'] = (np.float32( + 0.0), 'angle between the Moon and opt axis at CABEND') + ofd_r.header['TEL_ALT1'] = (np.float64( + 0.0), 'angle between opt axis and Elimb at CABEND') + ofd_r.header['POS_ANG1'] = ( + float(0.0), 'angle between y axis and North Pole at CABEND') + + ############################################################## + + ofd_r.header['POSI1_X'] = ( + float(p1x), 'orbital position in X at CABEND (km)') + ofd_r.header['POSI1_Y'] = ( + float(p1y), 'orbital position in Y at CABEND (km)') + ofd_r.header['POSI1_Z'] = ( + float(p1z), 'orbital position in Z at CABEND (km)') + + ofd_r.header['VELO1_X'] = ( + float(p1vx), 'orbital velocity in X at CABEND (km/s)') + ofd_r.header['VELO1_Y'] = ( + float(p1vy), 'orbital velocity in Y at CABEND (km/s)') + ofd_r.header['VELO1_Z'] = ( + float(p1vz), 'orbital velocity in Z at CABEND (km/s)') + + ofd_r.header['Euler1_1'] = (np.float64( + 0.0), 'Euler angle 1 at CABEND (deg)') + ofd_r.header['Euler1_2'] = (np.float64( + 0.0), 'Euler angle 2 at CABEND (deg)') + ofd_r.header['Euler1_3'] = (np.float64( + 0.0), 'Euler angle 3 at CABEND (deg)') + ofd_r.header['RA_PNT1'] = ( + float(ofd_r.header['RA_PNT0']), 'pointing RA at CABEND (deg)') + ofd_r.header['DEC_PNT1'] = ( + float(ofd_r.header['DEC_PNT0']), 'pointing Dec at CABEND (deg)') + ofd_r.header['EPOCH'] = (np.float32( + time2jd(t3.utcnow())), 'equinox of pointing RA and Dec') + ofd_r.header['EXPTIME'] = (np.float32( + self.information['exptime']), 'exposure time (s)') + + ofd_r.header['CHECKSUM'] = (0, '0') + ################################################ + ########## finish header for 0 layer #################### + ########## finish header for 0 layer #################### + # finish header for 0 layer + + if HeaderTest == 'yes': + hdu_r = fits.ImageHDU() + else: + hdu_r = fits.ImageHDU(data=np.uint16(self.image_r)) + ### MCI header for image layer ### + # need review thie header + hdu_r.header['XTENSION'] = ('IMAGE', 'extension type') + hdu_r.header['BITPIX'] = (np.int16(16), 'array data type') + hdu_r.header['NAXIS'] = (np.int16(2), 'number of array dimensions') + + hdu_r.header['NAXIS1'] = (np.int32(4936), '') + hdu_r.header['NAXIS2'] = (np.int32(23984), '') + + hdu_r.header['PCOUNT'] = (np.int16(0), 'number of parameters') + hdu_r.header['GCOUNT'] = (np.int16(1), 'number of groups') + + hdu_r.header['BSCALE'] = (int(1), '') + hdu_r.header['BZERO'] = (int(32768), '') + + hdu_r.header['EXTNAME'] = ('SCI', '') + hdu_r.header['EXTVER'] = (np.int16(1), '') + + hdu_r.header['BUNIT'] = ('ADU', 'physical unit of array values') + + ###### file information ############### + + temp = t3.utcnow() + data_time = temp.strftime("%Y-%m-%dT%H:%M:%S.%f") + hdu_r.header['DATE'] = ( + data_time[:21], 'written date (yyyy-mm-ddThh:mm:ss.s)') + hdu_r.header['FILENAME'] = (filename_r[:68], '') + hdu_r.header['FILETYPE'] = (self.source[:12], 'observation type') + hdu_r.header['RADECSYS'] = ( + 'ICRS', 'coordinate system of the object') + hdu_r.header['EQUINOX'] = (float(2000.1), '') + hdu_r.header['FITSSWV'] = ('4.2.1', 'FITS creating software version') + + ######################################## + + #### instrument status and object information ###### + hdu_r.header['TELESCOP'] = ('CSST', 'telescope name') + hdu_r.header['INSTRUME'] = ('MCI', 'instrument name') + hdu_r.header['CHANNEL'] = ('C2', 'channel number') + hdu_r.header['FILTERNO'] = (filternum, 'filter number') + hdu_r.header['DIFFUSER'] = ( + bool(True), 'insert diffuser status for flat calibration') + hdu_r.header['FLAMP'] = (np.int16(0), 'status of flat lamp') + hdu_r.header['MCISTAT'] = ( + np.int16(0), 'MCI components status parameter') + + hdu_r.header['DATE-OBS'] = (data_time[:21], + 'observation date (yyyy-mm-ddThh:mm:ss.s)') + hdu_r.header['OBJECT'] = ('MCI_obj', 'object name') + + hdu_r.header['TARGET'] = ( + self.information['target'][:15], 'target name (hhmmss.s+ddmmss)') + ####################################################################### + hdu_r.header['OBSID'] = (str(obsid), 'observation ID') + hdu_r.header['RA_OBJ'] = ( + float(self.information['ra_obj']), 'object RA (deg)') + hdu_r.header['DEC_OBJ'] = ( + float(self.information['dec_obj']), 'object Dec (deg)') + + ##### detector and Filter information ##### + hdu_r.header['FILTER'] = (self.filter_r[:6], 'filter band') + ##### Detector information #### + hdu_r.header['DETSN'] = ( + 'E2V-CCD-290-0000000', 'detector serial number') + hdu_r.header['DETNAME'] = ('blue-opt', 'detector name') + # need review here + + hdu_r.header['DETTEMP0'] = (np.float32( + 0), 'detector temperature at EXPSTART (K)') + hdu_r.header['DETTEMP1'] = (np.float32( + 0), 'detector temperature at EXPEND (K)') + hdu_r.header['DETTEMP2'] = (np.float32( + 0), 'detector temperature at READT1 (K)') + + hdu_r.header['DETSIZE'] = (str( + self.information['ysize'])+'*'+str(self.information['xsize']), 'detector size') + hdu_r.header['DATASECT'] = ( + str(self.image_g.shape[0])+'*'+str(self.image_g.shape[1]), 'data section') + + hdu_r.header['PIXSCAL1'] = ( + float(0.05), 'pixel scale for axis 1 (arcsec/pixel)') + hdu_r.header['PIXSCAL2'] = ( + float(0.05), 'pixel scale for axis 2 (arcsec/pixel)') + + hdu_r.header['PIXSIZE1'] = (int(10), 'pixel size for axis 1 (micron)') + hdu_r.header['PIXSIZE2'] = (int(10), 'pixel size for axis 2 (micron)') + + hdu_r.header['NCHANNEL'] = (np.int16(16), 'number of readout channels') + + hdu_r.header['PSCAN1'] = ( + np.int32(27), 'horizontal prescan width, per readout channel') + hdu_r.header['PSCAN2'] = ( + np.int32(0), 'vertical prescan width, per readout channel') + + hdu_r.header['OSCAN1'] = ( + np.int32(320), 'horizontal overscan width, per readout channel') + hdu_r.header['OSCAN2'] = ( + np.int32(320), 'vertical overscan width, per readout channel') + + hdu_r.header['BIN_X'] = (np.int16(1), 'bin number in X') + hdu_r.header['BIN_Y'] = (np.int16(1), 'bin number in Y') + + #### World coordinate system information ##### + + hdu_r.header['WCSAXES'] = (np.int16(2), 'number of WCS axes') + hdu_r.header['CRPIX1'] = ( + round(float(self.information['CRPIX1']), 1), 'x-coordinate of reference pixel') + hdu_r.header['CRPIX2'] = ( + round(float(self.information['CRPIX2']), 1), 'y-coordinate of reference pixel') + + hdu_r.header['CRVAL1'] = ( + float(self.information['CRVAL1']), 'value of reference pixel on axis 1') + hdu_r.header['CRVAL2'] = ( + float(self.information['CRVAL2']), 'value of reference pixel on axis 2') + + hdu_r.header['CTYPE1'] = ('RA---TAN', 'type of RA WCS projection') + hdu_r.header['CTYPE2'] = ('DEC--TAN', 'type of Dec WCS projection') + + hdu_r.header['CD1_1'] = ( + float(self.information['CD1_1']), 'transformation matrix element (deg/pix)') + hdu_r.header['CD1_2'] = ( + float(self.information['CD1_2']), 'transformation matrix element (deg/pix)') + hdu_r.header['CD2_1'] = ( + float(self.information['CD2_1']), 'transformation matrix element (deg/pix)') + hdu_r.header['CD2_2'] = ( + float(self.information['CD2_2']), 'transformation matrix element (deg/pix)') + + ###### Readout information ####### + hdu_r.header['GAINLVL'] = (np.float32( + 1.5), 'gain level (e-/ADU)') + hdu_r.header['GAIN01'] = (np.float32( + self.information['r_gain1']), 'CCD gain [channel 1] (e-/ADU)') + hdu_r.header['GAIN02'] = (np.float32( + self.information['r_gain2']), 'CCD gain [channel 2] (e-/ADU)') + hdu_r.header['GAIN03'] = (np.float32( + self.information['r_gain3']), 'CCD gain [channel 3] (e-/ADU)') + hdu_r.header['GAIN04'] = (np.float32( + self.information['r_gain4']), 'CCD gain [channel 4] (e-/ADU)') + hdu_r.header['GAIN05'] = (np.float32( + self.information['r_gain5']), 'CCD gain [channel 5] (e-/ADU)') + hdu_r.header['GAIN06'] = (np.float32( + self.information['r_gain6']), 'CCD gain [channel 6] (e-/ADU)') + hdu_r.header['GAIN07'] = (np.float32( + self.information['r_gain7']), 'CCD gain [channel 7] (e-/ADU)') + hdu_r.header['GAIN08'] = (np.float32( + self.information['r_gain8']), 'CCD gain [channel 8] (e-/ADU)') + hdu_r.header['GAIN09'] = (np.float32( + self.information['r_gain9']), 'CCD gain [channel 9] (e-/ADU)') + hdu_r.header['GAIN10'] = (np.float32( + self.information['r_gain10']), 'CCD gain [channel 10] (e-/ADU)') + hdu_r.header['GAIN11'] = (np.float32( + self.information['r_gain11']), 'CCD gain [channel 11] (e-/ADU)') + hdu_r.header['GAIN12'] = (np.float32( + self.information['r_gain12']), 'CCD gain [channel 12] (e-/ADU)') + hdu_r.header['GAIN13'] = (np.float32( + self.information['r_gain13']), 'CCD gain [channel 13] (e-/ADU)') + hdu_r.header['GAIN14'] = (np.float32( + self.information['r_gain14']), 'CCD gain [channel 14] (e-/ADU)') + hdu_r.header['GAIN15'] = (np.float32( + self.information['r_gain15']), 'CCD gain [channel 15] (e-/ADU)') + hdu_r.header['GAIN16'] = (np.float32( + self.information['r_gain16']), 'CCD gain [channel 16] (e-/ADU)') + ####### + hdu_r.header['RON01'] = (np.float32( + self.information['r_rdnois1']), 'readout noise [channel 01] (e-)') + hdu_r.header['RON02'] = (np.float32( + self.information['r_rdnois2']), 'readout noise [channel 02] (e-)') + hdu_r.header['RON03'] = (np.float32( + self.information['r_rdnois3']), 'readout noise [channel 03] (e-)') + hdu_r.header['RON04'] = (np.float32( + self.information['r_rdnois4']), 'readout noise [channel 04] (e-)') + hdu_r.header['RON05'] = (np.float32( + self.information['r_rdnois5']), 'readout noise [channel 05] (e-)') + hdu_r.header['RON06'] = (np.float32( + self.information['r_rdnois6']), 'readout noise [channel 06] (e-)') + hdu_r.header['RON07'] = (np.float32( + self.information['r_rdnois7']), 'readout noise [channel 07] (e-)') + hdu_r.header['RON08'] = (np.float32( + self.information['r_rdnois8']), 'readout noise [channel 08] (e-)') + hdu_r.header['RON09'] = (np.float32( + self.information['r_rdnois9']), 'readout noise [channel 09] (e-)') + hdu_r.header['RON10'] = (np.float32( + self.information['r_rdnois10']), 'readout noise [channel 10] (e-)') + hdu_r.header['RON11'] = (np.float32( + self.information['r_rdnois11']), 'readout noise [channel 11] (e-)') + hdu_r.header['RON12'] = (np.float32( + self.information['r_rdnois12']), 'readout noise [channel 12] (e-)') + hdu_r.header['RON13'] = (np.float32( + self.information['r_rdnois13']), 'readout noise [channel 13] (e-)') + hdu_r.header['RON14'] = (np.float32( + self.information['r_rdnois14']), 'readout noise [channel 14] (e-)') + hdu_r.header['RON15'] = (np.float32( + self.information['r_rdnois15']), 'readout noise [channel 15] (e-)') + hdu_r.header['RON16'] = (np.float32( + self.information['r_rdnois16']), 'readout noise [channel 16] (e-)') + + ####### + hdu_r.header['DETBIA01'] = (np.float32( + self.information['r_detbia1']), 'amplifier bias grey value [channel 1] (ADU)') + hdu_r.header['DETBIA02'] = (np.float32( + self.information['r_detbia2']), 'amplifier bias grey value [channel 2] (ADU)') + hdu_r.header['DETBIA03'] = (np.float32( + self.information['r_detbia3']), 'amplifier bias grey value [channel 3] (ADU)') + hdu_r.header['DETBIA04'] = (np.float32( + self.information['r_detbia4']), 'amplifier bias grey value [channel 4] (ADU)') + hdu_r.header['DETBIA05'] = (np.float32( + self.information['r_detbia5']), 'amplifier bias grey value [channel 5] (ADU)') + hdu_r.header['DETBIA06'] = (np.float32( + self.information['r_detbia6']), 'amplifier bias grey value [channel 6] (ADU)') + hdu_r.header['DETBIA07'] = (np.float32( + self.information['r_detbia7']), 'amplifier bias grey value [channel 7] (ADU)') + hdu_r.header['DETBIA08'] = (np.float32( + self.information['r_detbia8']), 'amplifier bias grey value [channel 8] (ADU)') + hdu_r.header['DETBIA09'] = (np.float32( + self.information['r_detbia9']), 'amplifier bias grey value [channel 9] (ADU)') + hdu_r.header['DETBIA10'] = (np.float32( + self.information['r_detbia10']), 'amplifier bias grey value [channel 10] (ADU)') + hdu_r.header['DETBIA11'] = (np.float32( + self.information['r_detbia11']), 'amplifier bias grey value [channel 11] (ADU)') + hdu_r.header['DETBIA12'] = (np.float32( + self.information['r_detbia12']), 'amplifier bias grey value [channel 12] (ADU)') + hdu_r.header['DETBIA13'] = (np.float32( + self.information['r_detbia13']), 'amplifier bias grey value [channel 13] (ADU)') + hdu_r.header['DETBIA14'] = (np.float32( + self.information['r_detbia14']), 'amplifier bias grey value [channel 14] (ADU)') + hdu_r.header['DETBIA15'] = (np.float32( + self.information['r_detbia15']), 'amplifier bias grey value [channel 15] (ADU)') + hdu_r.header['DETBIA16'] = (np.float32( + self.information['r_detbia16']), 'amplifier bias grey value [channel 16] (ADU)') + hdu_r.header['ROSPEED'] = (np.float32(100), 'readout speed (MHz)') + + ###################################### + #### exposure and shutter information ##### + hdu_r.header['EXPSTART'] = (np.float64( + time2mjd(self.dt)), 'exposure start time (MJD)') + hdu_r.header['EXPEND'] = (float(time2mjd(t2)), + 'exposure end time (MJD)') + hdu_r.header['EXPTIME'] = (np.float32( + self.information['exptime']), 'exposure time (s)') + hdu_r.header['DARKTIME'] = (np.float32( + self.information['exptime']), 'dark current time (s)') + hdu_r.header['SHTSTAT'] = (bool(True), 'shutter status') + + ###### satellite and its allitide information ############## + hdu_r.header['REFFRAME'] = ('CSSTGSC-1.0', 'guiding catalog version') + hdu_r.header['SATESWV'] = ( + 'softwave-1.0', 'satellite software version') + hdu_r.header['CABSTART'] = (np.float64( + time2jd(self.dt)), 'first cabin time after exposure start (MJD)') + hdu_r.header['SUNANGL0'] = (np.float32( + 0.0), 'angle between the Sun and opt axis at CABSTART') + hdu_r.header['MOONANG0'] = (np.float32( + 0.0), 'angle between the Moon and opt axis at CABSTART') + hdu_r.header['TEL_ALT0'] = (np.float64( + 0.0), 'angle between opt axis and Elimb at CABSTART') + hdu_r.header['POS_ANG0'] = (np.float64( + 0.0), 'angle between y axis and North Pole at CABSTART') + hdu_r.header['POSI0_X'] = (np.float64( + self.information['pos_x']), 'orbital position in X at CABSTART (km)') + hdu_r.header['POSI0_Y'] = (np.float64( + self.information['pos_y']), 'orbital position in Y at CABSTART (km)') + hdu_r.header['POSI0_Z'] = (np.float64( + self.information['pos_z']), 'orbital position in Z at CABSTART (km)') + hdu_r.header['VELO0_X'] = (np.float64( + self.information['velocity_x']), 'orbital velocity in X at CABSTART (km/s)') + hdu_r.header['VELO0_Y'] = (np.float64( + self.information['velocity_y']), 'orbital velocity in Y at CABSTART (km/s)') + hdu_r.header['VELO0_Z'] = (np.float64( + self.information['velocity_z']), 'orbital velocity in Z at CABSTART (km/s)') + hdu_r.header['Euler0_1'] = (np.float64( + 0.0), 'Euler angle 1 at CABSTART (deg)') + hdu_r.header['Euler0_2'] = (np.float64( + 0.0), 'Euler angle 2 at CABSTART (deg)') + hdu_r.header['Euler0_3'] = (np.float64( + 0.0), 'Euler angle 3 at CABSTART (deg)') + hdu_r.header['RA_PNT0'] = ( + float(self.information['ra_pnt0']), 'pointing RA at CABSTART (deg)') + hdu_r.header['DEC_PNT0'] = ( + float(self.information['dec_pnt0']), 'pointing Dec at CABSTART (deg)') + hdu_r.header['CABEND'] = (float(time2jd(t2)), + 'first cabin time after exposure end (MJD)') + hdu_r.header['SUNANGL1'] = (np.float32( + 0.0), 'angle between the Sun and opt axis at CABEND') + hdu_r.header['MOONANG1'] = (np.float32( + 0.0), 'angle between the Moon and opt axis at CABEND') + hdu_r.header['TEL_ALT1'] = (np.float64( + 0.0), 'angle between opt axis and Elimb at CABEND') + hdu_r.header['POS_ANG1'] = ( + float(0.0), 'angle between y axis and North Pole at CABEND') + hdu_r.header['POSI1_X'] = ( + float(p1x), 'orbital position in X at CABEND (km)') + hdu_r.header['POSI1_Y'] = ( + float(p1y), 'orbital position in Y at CABEND (km)') + hdu_r.header['POSI1_Z'] = ( + float(p1z), 'orbital position in Z at CABEND (km)') + hdu_r.header['VELO1_X'] = ( + float(p1vx), 'orbital velocity in X at CABEND (km/s)') + hdu_r.header['VELO1_Y'] = ( + float(p1vy), 'orbital velocity in Y at CABEND (km/s)') + hdu_r.header['VELO1_Z'] = ( + float(p1vz), 'orbital velocity in Z at CABEND (km/s)') + hdu_r.header['Euler1_1'] = (np.float64( + 0.0), 'Euler angle 1 at CABEND (deg)') + hdu_r.header['Euler1_2'] = (np.float64( + 0.0), 'Euler angle 2 at CABEND (deg)') + hdu_r.header['Euler1_3'] = (np.float64( + 0.0), 'Euler angle 3 at CABEND (deg)') + hdu_r.header['RA_PNT1'] = ( + float(hdu_r.header['RA_PNT0']), 'pointing RA at CABEND (deg)') + hdu_r.header['DEC_PNT1'] = ( + float(hdu_r.header['DEC_PNT0']), 'pointing Dec at CABEND (deg)') + hdu_r.header['EPOCH'] = (np.float32( + time2jd(t3.utcnow())), 'equinox of pointing RA and Dec') + hdu_r.header['CHECKSUM'] = (0, '0') + ############################################################## + # finish MCI header for image layer in G channel + # 333 + + hdulist_r = fits.HDUList([ofd_r, hdu_r]) + #### + ofd_r.header.add_comment( + '========================================================================', after='FITSSWV') + ofd_r.header.add_comment('OBJECT INFORMATION', after='FITSSWV') + ofd_r.header.add_comment( + '========================================================================', after='FITSSWV') + + ofd_r.header.add_comment( + '========================================================================', after='DEC_OBJ') + ofd_r.header.add_comment('TELESCOPE INFORMATION', after='DEC_OBJ') + ofd_r.header.add_comment( + '========================================================================', after='DEC_OBJ') + + ofd_r.header.add_comment( + '========================================================================', after='EXPTIME') + ofd_r.header.add_comment('VERIFICATION INFORMATION', after='EXPTIME') + ofd_r.header.add_comment( + '========================================================================', after='EXPTIME') + + ####################################################################### + hdu_r.header.add_comment( + '========================================================================', after='BUNIT') + hdu_r.header.add_comment('FILE INFORMATION', after='BUNIT') + hdu_r.header.add_comment( + '========================================================================', after='BUNIT') + + hdu_r.header.add_comment( + '========================================================================', after='FITSSWV') + hdu_r.header.add_comment( + 'INSTRUMENT STATUS AND OBJECT INFORMATION', after='FITSSWV') + hdu_r.header.add_comment( + '========================================================================', after='FITSSWV') + + hdu_r.header.add_comment( + '========================================================================', after='DEC_OBJ') + hdu_r.header.add_comment( + 'DETECTOR AND FILTER INFORMATION', after='DEC_OBJ') + hdu_r.header.add_comment( + '========================================================================', after='DEC_OBJ') + + hdu_r.header.add_comment( + '========================================================================', after='BIN_Y') + hdu_r.header.add_comment( + 'WORLD COORDINATE SYSTEM INFORMATION', after='BIN_Y') + hdu_r.header.add_comment( + '========================================================================', after='BIN_Y') + + hdu_r.header.add_comment( + '========================================================================', after='CD2_2') + hdu_r.header.add_comment('READOUT INFORMATION', after='CD2_2') + hdu_r.header.add_comment( + '========================================================================', after='CD2_2') + + hdu_r.header.add_comment( + '========================================================================', after='ROSPEED') + hdu_r.header.add_comment( + 'EXPOSURE AND SHUTTER INFORMATION', after='ROSPEED') + hdu_r.header.add_comment( + '========================================================================', after='ROSPEED') + + hdu_r.header.add_comment( + '========================================================================', after='SHTSTAT') + hdu_r.header.add_comment( + 'SATELLITE AND ITS ATTITUDE INFORMATION', after='SHTSTAT') + hdu_r.header.add_comment( + '========================================================================', after='SHTSTAT') + + hdu_r.header.add_comment( + '========================================================================', after='EPOCH') + hdu_r.header.add_comment('VERIFICATION INFORMATION', after='EPOCH') + hdu_r.header.add_comment( + '========================================================================', after='EPOCH') + + ##################################### + # finish MCI header for image layer in R channel + + if self.source == 'DARK' or self.source == 'FLAT' or self.source == 'BIAS': + file_r = self.result_path+'/cali_Data/'+filename_r+'.fits' + else: + file_r = self.result_path+'/sky_Data/'+filename_r + '.fits' + + hdulist_r.writeto(file_r, checksum=True) + + print('MCI_r.fits is created ') + + ################################################################################# + # i band + + # create a new FITS file, using HDUList instance + ofd_i = fits.PrimaryHDU() + + if self.filterP[self.filter_i]['number'] > 9: + filternum = str(self.filterP[self.filter_i]['number']) + else: + filternum = '0'+str(self.filterP[self.filter_i]['number']) + + filename_i = 'CSST_MCI_C3_'+self.source+'_'+exp_starttime+'_'+exp_endtime+'_'+str(obsid)+'_' + \ + filternum+'_L0_V01' + + # add primary header + + ofd_i.header['SIMPLE'] = (True, 'conforms to FITS standard') + ofd_i.header['BITPIX'] = (0, 'array data type') + ofd_i.header['NAXIS'] = (0, 'number of array dimensions') + ofd_i.header['NEXTEND'] = (np.int16(1), 'number of array dimensions') + ############################### + temp = t3.utcnow() + data_time = temp.strftime("%Y-%m-%dT%H:%M:%S.%f") + ofd_i.header['DATE'] = ( + data_time[:21], 'written date (yyyy-mm-ddThh:mm:ss.s)') + ofd_i.header['FILENAME'] = (filename_i[:68], '') + ofd_i.header['FILETYPE'] = (self.source[:12], 'observation type') + ofd_i.header['TELESCOP'] = ('CSST', 'telescope name') + ofd_i.header['INSTRUME'] = ('MCI', 'instrument name') + ofd_i.header['RADECSYS'] = ( + 'ICRS', 'coordinate system of the object') + ofd_i.header['EQUINOX'] = (float(2000.0), '') + ofd_i.header['FITSSWV'] = ( + 'mci_sim_0.8.03', 'FITS creating software version') + ############################################################################## + ######### Object information ######################################### + ofd_i.header['OBJECT'] = ( + self.information['name_obj'][:30], 'object name') + ofd_i.header['TARGET'] = ( + self.information['target'][:15], 'target name (hhmmss.s+ddmmss)') + ####################################################################### + ofd_i.header['OBSID'] = (str(obsid), 'observation ID') + ofd_i.header['RA_OBJ'] = ( + float(self.information['ra_obj']), 'object RA (deg)') + ofd_i.header['DEC_OBJ'] = ( + float(self.information['dec_obj']), 'object Dec (deg)') + + ######## Telescope information ############### + ofd_i.header['REFFRAME'] = ('CSSTGSC-1.0', 'guiding catalog version') + ofd_i.header['DATE-OBS'] = (data_time[:21], + 'observation date (yyyy-mm-ddThh:mm:ss.s)') + ofd_i.header['SATESWV'] = ( + 'softwave-1.0', 'satellite software version') + ofd_i.header['EXPSTART'] = (np.float64( + time2mjd(self.dt)), 'exposure start time (MJD)') + ofd_i.header['CABSTART'] = (np.float64( + time2jd(self.dt)), 'first cabin time after exposure start (MJD)') + ofd_i.header['SUNANGL0'] = (np.float32( + 0.0), 'angle between the Sun and opt axis at CABSTART') + ofd_i.header['MOONANG0'] = (np.float32( + 0.0), 'angle between the Moon and opt axis at CABSTART') + ofd_i.header['TEL_ALT0'] = (np.float64( + 0.0), 'angle between opt axis and Elimb at CABSTART') + ofd_i.header['POS_ANG0'] = (np.float64( + 0.0), 'angle between y axis and North Pole at CABSTART') + ofd_i.header['POSI0_X'] = (np.float64( + self.information['pos_x']), 'orbital position in X at CABSTART (km)') + ofd_i.header['POSI0_Y'] = (np.float64( + self.information['pos_y']), 'orbital position in Y at CABSTART (km)') + ofd_i.header['POSI0_Z'] = (np.float64( + self.information['pos_z']), 'orbital position in Z at CABSTART (km)') + ofd_i.header['VELO0_X'] = (np.float64( + self.information['velocity_x']), 'orbital velocity in X at CABSTART (km/s)') + ofd_i.header['VELO0_Y'] = (np.float64( + self.information['velocity_y']), 'orbital velocity in Y at CABSTART (km/s)') + ofd_i.header['VELO0_Z'] = (np.float64( + self.information['velocity_z']), 'orbital velocity in Z at CABSTART (km/s)') + ofd_i.header['Euler0_1'] = (np.float64( + 0.0), 'Euler angle 1 at CABSTART (deg)') + ofd_i.header['Euler0_2'] = (np.float64( + 0.0), 'Euler angle 2 at CABSTART (deg)') + ofd_i.header['Euler0_3'] = (np.float64( + 0.0), 'Euler angle 3 at CABSTART (deg)') + + ofd_i.header['RA_PNT0'] = ( + float(self.information['ra_pnt0']), 'pointing RA at CABSTART (deg)') + ofd_i.header['DEC_PNT0'] = ( + float(self.information['dec_pnt0']), 'pointing Dec at CABSTART (deg)') + ofd_i.header['EXPEND'] = (float(time2mjd(t2)), + 'exposure end time (MJD)') + ofd_i.header['CABEND'] = (float(time2jd(t2)), + 'first cabin time after exposure end (MJD)') + ofd_i.header['SUNANGL1'] = (np.float32( + 0.0), 'angle between the Sun and opt axis at CABEND') + ofd_i.header['MOONANG1'] = (np.float32( + 0.0), 'angle between the Moon and opt axis at CABEND') + ofd_i.header['TEL_ALT1'] = (np.float64( + 0.0), 'angle between opt axis and Elimb at CABEND') + ofd_i.header['POS_ANG1'] = ( + float(0.0), 'angle between y axis and North Pole at CABEND') + + ############################################################## + + ofd_i.header['POSI1_X'] = ( + float(p1x), 'orbital position in X at CABEND (km)') + ofd_i.header['POSI1_Y'] = ( + float(p1y), 'orbital position in Y at CABEND (km)') + ofd_i.header['POSI1_Z'] = ( + float(p1z), 'orbital position in Z at CABEND (km)') + + ofd_i.header['VELO1_X'] = ( + float(p1vx), 'orbital velocity in X at CABEND (km/s)') + ofd_i.header['VELO1_Y'] = ( + float(p1vy), 'orbital velocity in Y at CABEND (km/s)') + ofd_i.header['VELO1_Z'] = ( + float(p1vz), 'orbital velocity in Z at CABEND (km/s)') + + ################### + + ofd_i.header['Euler1_1'] = (np.float64( + 0.0), 'Euler angle 1 at CABEND (deg)') + ofd_i.header['Euler1_2'] = (np.float64( + 0.0), 'Euler angle 2 at CABEND (deg)') + ofd_i.header['Euler1_3'] = (np.float64( + 0.0), 'Euler angle 3 at CABEND (deg)') + ofd_i.header['RA_PNT1'] = ( + float(ofd_i.header['RA_PNT0']), 'pointing RA at CABEND (deg)') + ofd_i.header['DEC_PNT1'] = ( + float(ofd_i.header['DEC_PNT0']), 'pointing Dec at CABEND (deg)') + ofd_i.header['EPOCH'] = (np.float32( + time2jd(t3.utcnow())), 'equinox of pointing RA and Dec') + ofd_i.header['EXPTIME'] = (np.float32( + self.information['exptime']), 'exposure time (s)') + + ofd_i.header['CHECKSUM'] = (0, '0') + ################################################ + ########## finish header for 0 layer #################### + + # finish header for 0 layer + + # new image HDU + if HeaderTest == 'yes': + hdu_i = fits.ImageHDU() + else: + hdu_i = fits.ImageHDU(data=np.uint16(self.image_i)) + + # need review thie header + ### MCI header for image layer ### + # need review thie header + hdu_i.header['XTENSION'] = ('IMAGE', 'extension type') + hdu_i.header['BITPIX'] = (np.int16(16), 'array data type') + hdu_i.header['NAXIS'] = (np.int16(2), 'number of array dimensions') + + hdu_i.header['NAXIS1'] = (np.int32(4936), '') + hdu_i.header['NAXIS2'] = (np.int32(23984), '') + + hdu_i.header['PCOUNT'] = (np.int16(0), 'number of parameters') + hdu_i.header['GCOUNT'] = (np.int16(1), 'number of groups') + + hdu_i.header['BSCALE'] = (int(1), '') + hdu_i.header['BZERO'] = (int(32768), '') + + hdu_i.header['EXTNAME'] = ('SCI', '') + hdu_i.header['EXTVER'] = (np.int16(1), '') + + hdu_i.header['BUNIT'] = ('ADU', 'physical unit of array values') + + ###### file information ############### + + temp = t3.utcnow() + data_time = temp.strftime("%Y-%m-%dT%H:%M:%S.%f") + hdu_i.header['DATE'] = ( + data_time[:21], 'written date (yyyy-mm-ddThh:mm:ss.s)') + hdu_i.header['FILENAME'] = (filename_i[:68], '') + hdu_i.header['FILETYPE'] = (self.source[:12], 'observation type') + hdu_i.header['RADECSYS'] = ( + 'ICRS', 'coordinate system of the object') + hdu_i.header['EQUINOX'] = (float(2000.1), '') + hdu_i.header['FITSSWV'] = ('4.2.1', 'FITS creating software version') + + ######################################## + + #### instrument status and object information ###### + hdu_i.header['TELESCOP'] = ('CSST', 'telescope name') + hdu_i.header['INSTRUME'] = ('MCI', 'instrument name') + hdu_i.header['CHANNEL'] = ('C3', 'channel number') + hdu_i.header['FILTERNO'] = (filternum, 'filter number') + hdu_i.header['DIFFUSER'] = ( + bool(True), 'insert diffuser status for flat calibration') + hdu_i.header['FLAMP'] = (np.int16(0), 'status of flat lamp') + hdu_i.header['MCISTAT'] = ( + np.int16(0), 'MCI components status parameter') + + hdu_i.header['DATE-OBS'] = (data_time[:21], + 'observation date (yyyy-mm-ddThh:mm:ss.s)') + hdu_i.header['OBJECT'] = ('MCI_obj', 'object name') + + hdu_i.header['TARGET'] = ( + self.information['target'][:15], 'target name (hhmmss.s+ddmmss)') + ####################################################################### + hdu_i.header['OBSID'] = (str(obsid), 'observation ID') + hdu_i.header['RA_OBJ'] = ( + float(self.information['ra_obj']), 'object RA (deg)') + hdu_i.header['DEC_OBJ'] = ( + float(self.information['dec_obj']), 'object Dec (deg)') + + ##### detector and Filter information ##### + hdu_i.header['FILTER'] = (self.filter_i[:6], 'filter band') + ##### Detector information #### + hdu_i.header['DETSN'] = ( + 'E2V-CCD-290-0000000', 'detector serial number') + hdu_i.header['DETNAME'] = ('opt-red', 'detector name') + # need review here + + hdu_i.header['DETTEMP0'] = (np.float32( + 0), 'detector temperature at EXPSTART (K)') + hdu_i.header['DETTEMP1'] = (np.float32( + 0), 'detector temperature at EXPEND (K)') + hdu_i.header['DETTEMP2'] = (np.float32( + 0), 'detector temperature at READT1 (K)') + + hdu_i.header['DETSIZE'] = (str( + self.information['ysize'])+'*'+str(self.information['xsize']), 'detector size') + hdu_i.header['DATASECT'] = ( + str(self.image_g.shape[0])+'*'+str(self.image_g.shape[1]), 'data section') + + hdu_i.header['PIXSCAL1'] = ( + float(0.05), 'pixel scale for axis 1 (arcsec/pixel)') + hdu_i.header['PIXSCAL2'] = ( + float(0.05), 'pixel scale for axis 2 (arcsec/pixel)') + + hdu_i.header['PIXSIZE1'] = (int(10), 'pixel size for axis 1 (micron)') + hdu_i.header['PIXSIZE2'] = (int(10), 'pixel size for axis 2 (micron)') + + hdu_i.header['NCHANNEL'] = (np.int16(16), 'number of readout channels') + + hdu_i.header['PSCAN1'] = ( + np.int32(27), 'horizontal prescan width, per readout channel') + hdu_i.header['PSCAN2'] = ( + np.int32(0), 'vertical prescan width, per readout channel') + + hdu_i.header['OSCAN1'] = ( + np.int32(320), 'horizontal overscan width, per readout channel') + hdu_i.header['OSCAN2'] = ( + np.int32(320), 'vertical overscan width, per readout channel') + + hdu_i.header['BIN_X'] = (np.int16(1), 'bin number in X') + hdu_i.header['BIN_Y'] = (np.int16(1), 'bin number in Y') + + #### World coordinate system information ##### + + hdu_i.header['WCSAXES'] = (np.int16(2), 'number of WCS axes') + hdu_i.header['CRPIX1'] = ( + round(float(self.information['CRPIX1']), 1), 'x-coordinate of reference pixel') + hdu_i.header['CRPIX2'] = ( + round(float(self.information['CRPIX2']), 1), 'y-coordinate of reference pixel') + + hdu_i.header['CRVAL1'] = ( + float(self.information['CRVAL1']), 'value of reference pixel on axis 1') + hdu_i.header['CRVAL2'] = ( + float(self.information['CRVAL2']), 'value of reference pixel on axis 2') + + hdu_i.header['CTYPE1'] = ('RA---TAN', 'type of RA WCS projection') + hdu_i.header['CTYPE2'] = ('DEC--TAN', 'type of Dec WCS projection') + + hdu_i.header['CD1_1'] = ( + float(self.information['CD1_1']), 'transformation matrix element (deg/pix)') + hdu_i.header['CD1_2'] = ( + float(self.information['CD1_2']), 'transformation matrix element (deg/pix)') + hdu_i.header['CD2_1'] = ( + float(self.information['CD2_1']), 'transformation matrix element (deg/pix)') + hdu_i.header['CD2_2'] = ( + float(self.information['CD2_2']), 'transformation matrix element (deg/pix)') + + ###### Readout information ####### + hdu_i.header['GAINLVL'] = (np.float32( + 1.5), 'gain level (e-/ADU)') + hdu_i.header['GAIN01'] = (np.float32( + self.information['i_gain1']), 'CCD gain [channel 1] (e-/ADU)') + hdu_i.header['GAIN02'] = (np.float32( + self.information['i_gain2']), 'CCD gain [channel 2] (e-/ADU)') + hdu_i.header['GAIN03'] = (np.float32( + self.information['i_gain3']), 'CCD gain [channel 3] (e-/ADU)') + hdu_i.header['GAIN04'] = (np.float32( + self.information['i_gain4']), 'CCD gain [channel 4] (e-/ADU)') + hdu_i.header['GAIN05'] = (np.float32( + self.information['i_gain5']), 'CCD gain [channel 5] (e-/ADU)') + hdu_i.header['GAIN06'] = (np.float32( + self.information['i_gain6']), 'CCD gain [channel 6] (e-/ADU)') + hdu_i.header['GAIN07'] = (np.float32( + self.information['i_gain7']), 'CCD gain [channel 7] (e-/ADU)') + hdu_i.header['GAIN08'] = (np.float32( + self.information['i_gain8']), 'CCD gain [channel 8] (e-/ADU)') + hdu_i.header['GAIN09'] = (np.float32( + self.information['i_gain9']), 'CCD gain [channel 9] (e-/ADU)') + hdu_i.header['GAIN10'] = (np.float32( + self.information['i_gain10']), 'CCD gain [channel 10] (e-/ADU)') + hdu_i.header['GAIN11'] = (np.float32( + self.information['i_gain11']), 'CCD gain [channel 11] (e-/ADU)') + hdu_i.header['GAIN12'] = (np.float32( + self.information['i_gain12']), 'CCD gain [channel 12] (e-/ADU)') + hdu_i.header['GAIN13'] = (np.float32( + self.information['i_gain13']), 'CCD gain [channel 13] (e-/ADU)') + hdu_i.header['GAIN14'] = (np.float32( + self.information['i_gain14']), 'CCD gain [channel 14] (e-/ADU)') + hdu_i.header['GAIN15'] = (np.float32( + self.information['i_gain15']), 'CCD gain [channel 15] (e-/ADU)') + hdu_i.header['GAIN16'] = (np.float32( + self.information['i_gain16']), 'CCD gain [channel 16] (e-/ADU)') + ####### + hdu_i.header['RON01'] = (np.float32( + self.information['i_rdnois1']), 'readout noise [channel 01] (e-)') + hdu_i.header['RON02'] = (np.float32( + self.information['i_rdnois2']), 'readout noise [channel 02] (e-)') + hdu_i.header['RON03'] = (np.float32( + self.information['i_rdnois3']), 'readout noise [channel 03] (e-)') + hdu_i.header['RON04'] = (np.float32( + self.information['i_rdnois4']), 'readout noise [channel 04] (e-)') + hdu_i.header['RON05'] = (np.float32( + self.information['i_rdnois5']), 'readout noise [channel 05] (e-)') + hdu_i.header['RON06'] = (np.float32( + self.information['i_rdnois6']), 'readout noise [channel 06] (e-)') + hdu_i.header['RON07'] = (np.float32( + self.information['i_rdnois7']), 'readout noise [channel 07] (e-)') + hdu_i.header['RON08'] = (np.float32( + self.information['i_rdnois8']), 'readout noise [channel 08] (e-)') + hdu_i.header['RON09'] = (np.float32( + self.information['i_rdnois9']), 'readout noise [channel 09] (e-)') + hdu_i.header['RON10'] = (np.float32( + self.information['i_rdnois10']), 'readout noise [channel 10] (e-)') + hdu_i.header['RON11'] = (np.float32( + self.information['i_rdnois11']), 'readout noise [channel 11] (e-)') + hdu_i.header['RON12'] = (np.float32( + self.information['i_rdnois12']), 'readout noise [channel 12] (e-)') + hdu_i.header['RON13'] = (np.float32( + self.information['i_rdnois13']), 'readout noise [channel 13] (e-)') + hdu_i.header['RON14'] = (np.float32( + self.information['i_rdnois14']), 'readout noise [channel 14] (e-)') + hdu_i.header['RON15'] = (np.float32( + self.information['i_rdnois15']), 'readout noise [channel 15] (e-)') + hdu_i.header['RON16'] = (np.float32( + self.information['i_rdnois16']), 'readout noise [channel 16] (e-)') + + ####### + hdu_i.header['DETBIA01'] = (np.float32( + self.information['i_detbia1']), 'amplifier bias grey value [channel 1] (ADU)') + hdu_i.header['DETBIA02'] = (np.float32( + self.information['i_detbia2']), 'amplifier bias grey value [channel 2] (ADU)') + hdu_i.header['DETBIA03'] = (np.float32( + self.information['i_detbia3']), 'amplifier bias grey value [channel 3] (ADU)') + hdu_i.header['DETBIA04'] = (np.float32( + self.information['i_detbia4']), 'amplifier bias grey value [channel 4] (ADU)') + hdu_i.header['DETBIA05'] = (np.float32( + self.information['i_detbia5']), 'amplifier bias grey value [channel 5] (ADU)') + hdu_i.header['DETBIA06'] = (np.float32( + self.information['i_detbia6']), 'amplifier bias grey value [channel 6] (ADU)') + hdu_i.header['DETBIA07'] = (np.float32( + self.information['i_detbia7']), 'amplifier bias grey value [channel 7] (ADU)') + hdu_i.header['DETBIA08'] = (np.float32( + self.information['i_detbia8']), 'amplifier bias grey value [channel 8] (ADU)') + hdu_i.header['DETBIA09'] = (np.float32( + self.information['i_detbia9']), 'amplifier bias grey value [channel 9] (ADU)') + hdu_i.header['DETBIA10'] = (np.float32( + self.information['i_detbia10']), 'amplifier bias grey value [channel 10] (ADU)') + hdu_i.header['DETBIA11'] = (np.float32( + self.information['i_detbia11']), 'amplifier bias grey value [channel 11] (ADU)') + hdu_i.header['DETBIA12'] = (np.float32( + self.information['i_detbia12']), 'amplifier bias grey value [channel 12] (ADU)') + hdu_i.header['DETBIA13'] = (np.float32( + self.information['i_detbia13']), 'amplifier bias grey value [channel 13] (ADU)') + hdu_i.header['DETBIA14'] = (np.float32( + self.information['i_detbia14']), 'amplifier bias grey value [channel 14] (ADU)') + hdu_i.header['DETBIA15'] = (np.float32( + self.information['i_detbia15']), 'amplifier bias grey value [channel 15] (ADU)') + hdu_i.header['DETBIA16'] = (np.float32( + self.information['i_detbia16']), 'amplifier bias grey value [channel 16] (ADU)') + hdu_i.header['ROSPEED'] = (np.float32(100), 'readout speed (MHz)') + + ###################################### + #### exposure and shutter information ##### + hdu_i.header['EXPSTART'] = (np.float64( + time2mjd(self.dt)), 'exposure start time (MJD)') + hdu_i.header['EXPEND'] = (float(time2mjd(t2)), + 'exposure end time (MJD)') + hdu_i.header['EXPTIME'] = (np.float32( + self.information['exptime']), 'exposure time (s)') + hdu_i.header['DARKTIME'] = (np.float32( + self.information['exptime']), 'dark current time (s)') + hdu_i.header['SHTSTAT'] = (bool(True), 'shutter status') + + ###### satellite and its allitide information ############## + hdu_i.header['REFFRAME'] = ('CSSTGSC-1.0', 'guiding catalog version') + hdu_i.header['SATESWV'] = ( + 'softwave-1.0', 'satellite software version') + hdu_i.header['CABSTART'] = (np.float64( + time2jd(self.dt)), 'first cabin time after exposure start (MJD)') + hdu_i.header['SUNANGL0'] = (np.float32( + 0.0), 'angle between the Sun and opt axis at CABSTART') + hdu_i.header['MOONANG0'] = (np.float32( + 0.0), 'angle between the Moon and opt axis at CABSTART') + hdu_i.header['TEL_ALT0'] = (np.float64( + 0.0), 'angle between opt axis and Elimb at CABSTART') + hdu_i.header['POS_ANG0'] = (np.float64( + 0.0), 'angle between y axis and North Pole at CABSTART') + hdu_i.header['POSI0_X'] = (np.float64( + self.information['pos_x']), 'orbital position in X at CABSTART (km)') + hdu_i.header['POSI0_Y'] = (np.float64( + self.information['pos_y']), 'orbital position in Y at CABSTART (km)') + hdu_i.header['POSI0_Z'] = (np.float64( + self.information['pos_z']), 'orbital position in Z at CABSTART (km)') + hdu_i.header['VELO0_X'] = (np.float64( + self.information['velocity_x']), 'orbital velocity in X at CABSTART (km/s)') + hdu_i.header['VELO0_Y'] = (np.float64( + self.information['velocity_y']), 'orbital velocity in Y at CABSTART (km/s)') + hdu_i.header['VELO0_Z'] = (np.float64( + self.information['velocity_z']), 'orbital velocity in Z at CABSTART (km/s)') + hdu_i.header['Euler0_1'] = (np.float64( + 0.0), 'Euler angle 1 at CABSTART (deg)') + hdu_i.header['Euler0_2'] = (np.float64( + 0.0), 'Euler angle 2 at CABSTART (deg)') + hdu_i.header['Euler0_3'] = (np.float64( + 0.0), 'Euler angle 3 at CABSTART (deg)') + hdu_i.header['RA_PNT0'] = ( + float(self.information['ra_pnt0']), 'pointing RA at CABSTART (deg)') + hdu_i.header['DEC_PNT0'] = ( + float(self.information['dec_pnt0']), 'pointing Dec at CABSTART (deg)') + hdu_i.header['CABEND'] = (float(time2jd(t2)), + 'first cabin time after exposure end (MJD)') + hdu_i.header['SUNANGL1'] = (np.float32( + 0.0), 'angle between the Sun and opt axis at CABEND') + hdu_i.header['MOONANG1'] = (np.float32( + 0.0), 'angle between the Moon and opt axis at CABEND') + hdu_i.header['TEL_ALT1'] = (np.float64( + 0.0), 'angle between opt axis and Elimb at CABEND') + hdu_i.header['POS_ANG1'] = ( + float(0.0), 'angle between y axis and North Pole at CABEND') + hdu_i.header['POSI1_X'] = ( + float(p1x), 'orbital position in X at CABEND (km)') + hdu_i.header['POSI1_Y'] = ( + float(p1y), 'orbital position in Y at CABEND (km)') + hdu_i.header['POSI1_Z'] = ( + float(p1z), 'orbital position in Z at CABEND (km)') + hdu_i.header['VELO1_X'] = ( + float(p1vx), 'orbital velocity in X at CABEND (km/s)') + hdu_i.header['VELO1_Y'] = ( + float(p1vy), 'orbital velocity in Y at CABEND (km/s)') + hdu_i.header['VELO1_Z'] = ( + float(p1vz), 'orbital velocity in Z at CABEND (km/s)') + hdu_i.header['Euler1_1'] = (np.float64( + 0.0), 'Euler angle 1 at CABEND (deg)') + hdu_i.header['Euler1_2'] = (np.float64( + 0.0), 'Euler angle 2 at CABEND (deg)') + hdu_i.header['Euler1_3'] = (np.float64( + 0.0), 'Euler angle 3 at CABEND (deg)') + hdu_i.header['RA_PNT1'] = ( + float(hdu_i.header['RA_PNT0']), 'pointing RA at CABEND (deg)') + hdu_i.header['DEC_PNT1'] = ( + float(hdu_i.header['DEC_PNT0']), 'pointing Dec at CABEND (deg)') + hdu_i.header['EPOCH'] = (np.float32( + time2jd(t3.utcnow())), 'equinox of pointing RA and Dec') + hdu_i.header['CHECKSUM'] = (0, '0') + + ############################################################## + # finish MCI header for image layer in G channel + # 333 + + hdulist_i = fits.HDUList([ofd_i, hdu_i]) + #### + ofd_i.header.add_comment( + '========================================================================', after='FITSSWV') + ofd_i.header.add_comment('OBJECT INFORMATION', after='FITSSWV') + ofd_i.header.add_comment( + '========================================================================', after='FITSSWV') + + ofd_i.header.add_comment( + '========================================================================', after='DEC_OBJ') + ofd_i.header.add_comment('TELESCOPE INFORMATION', after='DEC_OBJ') + ofd_i.header.add_comment( + '========================================================================', after='DEC_OBJ') + + ofd_i.header.add_comment( + '========================================================================', after='EXPTIME') + ofd_i.header.add_comment('VERIFICATION INFORMATION', after='EXPTIME') + ofd_i.header.add_comment( + '========================================================================', after='EXPTIME') + + ####################################################################### + hdu_i.header.add_comment( + '========================================================================', after='BUNIT') + hdu_i.header.add_comment('FILE INFORMATION', after='BUNIT') + hdu_i.header.add_comment( + '========================================================================', after='BUNIT') + + hdu_i.header.add_comment( + '========================================================================', after='FITSSWV') + hdu_i.header.add_comment( + 'INSTRUMENT STATUS AND OBJECT INFORMATION', after='FITSSWV') + hdu_i.header.add_comment( + '========================================================================', after='FITSSWV') + + hdu_i.header.add_comment( + '========================================================================', after='DEC_OBJ') + hdu_i.header.add_comment( + 'DETECTOR AND FILTER INFORMATION', after='DEC_OBJ') + hdu_i.header.add_comment( + '========================================================================', after='DEC_OBJ') + + hdu_i.header.add_comment( + '========================================================================', after='BIN_Y') + hdu_i.header.add_comment( + 'WORLD COORDINATE SYSTEM INFORMATION', after='BIN_Y') + hdu_i.header.add_comment( + '========================================================================', after='BIN_Y') + + hdu_i.header.add_comment( + '========================================================================', after='CD2_2') + hdu_i.header.add_comment('READOUT INFORMATION', after='CD2_2') + hdu_i.header.add_comment( + '========================================================================', after='CD2_2') + + hdu_i.header.add_comment( + '========================================================================', after='ROSPEED') + hdu_i.header.add_comment( + 'EXPOSURE AND SHUTTER INFORMATION', after='ROSPEED') + hdu_i.header.add_comment( + '========================================================================', after='ROSPEED') + + hdu_i.header.add_comment( + '========================================================================', after='SHTSTAT') + hdu_i.header.add_comment( + 'SATELLITE AND ITS ATTITUDE INFORMATION', after='SHTSTAT') + hdu_i.header.add_comment( + '========================================================================', after='SHTSTAT') + + hdu_i.header.add_comment( + '========================================================================', after='EPOCH') + hdu_i.header.add_comment('VERIFICATION INFORMATION', after='EPOCH') + hdu_i.header.add_comment( + '========================================================================', after='EPOCH') + + # hdu_i.header.add_comment('========================================================================',before='CHECKSUM') + # hdu_i.header.add_comment('VERIFICATION INFORMATION',before='CHECKSUM') + # hdu_i.header.add_comment('========================================================================',before='CHECKSUM') + + ##################################### + # finish MCI header for image layer in i channel + # finish MCI header for image layer + if self.source == 'DARK' or self.source == 'FLAT' or self.source == 'BIAS': + file_i = self.result_path+'/cali_Data/'+filename_i+'.fits' + else: + file_i = self.result_path+'/sky_Data/'+filename_i + '.fits' + + hdulist_i.writeto(file_i, checksum=True) + +################################################################################################# + + def simulate(self, simnumber): + """ + Create a single simulated image of a quadrant defined by the configuration file. + Will do all steps defined in the config file sequentially. + + :return: None + """ + + ############################################################################################## + ############################################################################################## + sourcelist = ['EXDF', 'STAR', 'PIPR', 'TRNS', 'COMB', 'CALI'] + + sim_list = ['FLAT', 'BIAS', 'DARK', 'EXDF', + 'STAR', 'PIPR', 'TRNS', 'COMB', 'CALI'] + + if self.source not in sim_list: + print('Sourcein is wrong') + sys.exit(1) + + self._loadGhostModel() + self.load_filter_PSF() + self.information['simnumber'] = simnumber + + self.debug = self.information['debug'] + ########################### + flag = 0 + + self.dt_num = int((self.information['simnumber']-1)*( + self.information['exptime']+self.information['readouttime']+50)/120) + now_dt = datetime.utcnow() + now_jd = time2jd(now_dt) + + for k in range(1, 50, 1): + + fn = self.information['dir_path'] + \ + 'MCI_inputData/TianCe/orbit20160925/'+str(k)+'.txt' + d = np.loadtxt(fn) + + for kk in range(len(d[:, 0])): + if now_jd-d[kk, 0] <= 0: + flag = 1 + break + if flag == 1: + print('find time index...... ') + break + + ################################################################ + + if kk + self.dt_num < len(d[:, 0]): + + exptime_start_jd = d[kk+self.dt_num, 0] + + self.orbit_pars = d + self.orbit_file_num = k + self.orbit_exp_num = kk + + self.information['pos_x'] = float( + self.orbit_pars[self.orbit_exp_num, 1]) + self.information['pos_y'] = float( + self.orbit_pars[self.orbit_exp_num, 2]) + self.information['pos_z'] = float( + self.orbit_pars[self.orbit_exp_num, 3]) + self.information['velocity_x'] = float( + self.orbit_pars[self.orbit_exp_num, 4]) + self.information['velocity_y'] = float( + self.orbit_pars[self.orbit_exp_num, 5]) + self.information['velocity_z'] = float( + self.orbit_pars[self.orbit_exp_num, 6]) + + else: + + fn = self.information['dir_path'] + \ + 'MCI_inputData/TianCe/orbit20160925/'+str(k+1)+'.txt' + d = np.loadtxt(fn) + + self.orbit_pars = d + self.orbit_file_num = k+1 + self.orbit_exp_num = self.dt_num + + self.information['pos_x'] = float( + self.orbit_pars[self.orbit_exp_num, 1]) + self.information['pos_y'] = float( + self.orbit_pars[self.orbit_exp_num, 2]) + self.information['pos_z'] = float( + self.orbit_pars[self.orbit_exp_num, 3]) + self.information['velocity_x'] = float( + self.orbit_pars[self.orbit_exp_num, 4]) + self.information['velocity_y'] = float( + self.orbit_pars[self.orbit_exp_num, 5]) + self.information['velocity_z'] = float( + self.orbit_pars[self.orbit_exp_num, 6]) + ########################################## + + ###self.dt=julian.from_jd(exptime_start_jd, fmt='jd') + self.dt = jd2time(exptime_start_jd) + + # str(self.dt.year)+'-'+str(self.dt.month)+'-'+str(self.dt.day) + self.TianCe_day = str(self.dt.year)+self.dt.strftime("-%m-%d") + + self.TianCe_exp_start = dt2hmd(self.dt) + + ############################################################# + + if self.source in sourcelist: + + # self.earthshine_theta=30.0 # in degree + + ################################### + if self.sky_shift_rot: + + if self.information['simnumber'] > 1: + + np.random.seed(self.information['simnumber']) + ud = np.random.random() # Choose a random pointing ra in degree + dis_ra = 2*(ud-1)*10/3600.0 + + np.random.seed(10*self.information['simnumber']) + ud = np.random.random() # Choose a random pointing dec in degree + dis_dec = 2*(ud-1)*10/3600.0 + + np.random.seed(100*self.information['simnumber']) + ud = np.random.random() # Choose a random telescope rotation in degree + # the image will rotate in clockwise rotaion between -10-10 degree + rotTelPos = 2*(ud-1)*5.0*galsim.degrees + rotSkyPos = 0.0*galsim.degrees + + ################################ + # ### test code + # rotTelPos = 90.0*galsim.degrees + # dis_dec= -200*0.05/3600 + # dis_ra = 100*0.05/3600 + + else: + + dis_ra = 0 # + dis_dec = 0 # T + + # the image will rotate in clockwise rotaion between -10-10 degree + rotTelPos = 0.0*galsim.degrees + rotSkyPos = 0.0*galsim.degrees + + else: + + dis_ra = 0 # + dis_dec = 0 # T + + # the image will rotate in clockwise rotaion between -10-10 degree + rotTelPos = 0.0*galsim.degrees + rotSkyPos = 0.0*galsim.degrees + + ##################################################### + + # telescope pointing Ra displacement + self.information['T_disRa'] = dis_ra + # telescope pointing Dec displacement + self.information['T_disDec'] = dis_dec + + # telescope rotation in degree + self.information['rotTelPos'] = rotTelPos + # sky rotation in degree + self.information['rotSkyPos'] = rotSkyPos + + self.information['name_obj'] = 'MCI_obj' + + theta = self.information['rotTelPos'] - \ + self.information['rotSkyPos'] + + self.information['rotangle'] = theta.deg + #print('rotangle',theta.deg ) + + self.log.info('dis_ra(in pixel)=%f, dis_dec(in pixel)=%f, sky_rot(in deg)=%f' % ( + dis_ra*3600/0.05, dis_dec*3600/0.05, theta.deg)) + ####################################################################################################### + + ############################################################################# + ######### simulate star images from star SED data ################## + + if self.sim_star: + ########################### + + starimg = self.cal_StarSED_img() + + if 'starimg' in dir(): + self.log.info('Generate artifical stars finished...') + # fits.writeto(self.result_path+'/ori_Sky/starimg_C1_'+str(simnumber)+'.fits',starimg['g'], overwrite=True) + # fits.writeto(self.result_path+'/ori_Sky/starimg_C2_'+str(simnumber)+'.fits',starimg['r'], overwrite=True) + # fits.writeto(self.result_path+'/ori_Sky/starimg_C3_'+str(simnumber)+'.fits',starimg['i'], overwrite=True) + + self.img_fits_save( + starimg['g'], 'star_ori_C1_'+str(simnumber)) + self.img_fits_save( + starimg['r'], 'star_ori_C2_'+str(simnumber)) + self.img_fits_save( + starimg['i'], 'star_ori_C3_'+str(simnumber)) + + print('star image finished') + + ########################################################################## + # # ######################################################################## + # ###### generate galaxy SED image without Lensing effect ########## + + if self.sim_galaxy and self.source == 'EXDF': + + losimg = self.cal_Lensing_galaxy_img() + + if 'losimg' in dir(): + self.log.info('Generate los finished...') + print('losimage finished') + # fits.writeto(self.result_path+'/ori_Sky/galimg_C1_'+str(simnumber)+'.fits',losimg['g'], overwrite=True) + # fits.writeto(self.result_path+'/ori_Sky/galimg_C2_'+str(simnumber)+'.fits',losimg['r'], overwrite=True) + # fits.writeto(self.result_path+'/ori_Sky/galimg_C3_'+str(simnumber)+'.fits',losimg['i'], overwrite=True) + + self.img_fits_save( + losimg['g'], 'gal_ori_C1_'+str(simnumber)) + self.img_fits_save( + losimg['r'], 'gal_ori_C2_'+str(simnumber)) + self.img_fits_save( + losimg['i'], 'gal_ori_C3_'+str(simnumber)) + + ################################################################## + + # stack the image + + if 'losimg' in dir(): + self.image_g += losimg['g'] + self.image_r += losimg['r'] + self.image_i += losimg['i'] + + if 'starimg' in dir(): + self.image_g += starimg['g'] + self.image_r += starimg['r'] + self.image_i += starimg['i'] + + ########################################################################################################### + + fits.writeto(self.result_path+'/ori_Sky/original_C1_' + + str(simnumber)+'.fits', self.image_g, overwrite=True) + fits.writeto(self.result_path+'/ori_Sky/original_C2_' + + str(simnumber)+'.fits', self.image_r, overwrite=True) + fits.writeto(self.result_path+'/ori_Sky/original_C3_' + + str(simnumber)+'.fits', self.image_i, overwrite=True) + + # self.img_fits_save(self.image_g, 'full_ori_C1_'+str(simnumber)) + # self.img_fits_save(self.image_g, 'full_ori_C2_'+str(simnumber)) + # self.img_fits_save(self.image_g, 'full_ori_C3_'+str(simnumber)) + + ######################################################################################################### + + # end if source =='SKY' + + elif self.source == 'FLAT': + + theta = 0 + # np.ceil(self.information['simnumber']/2.0) + self.information['exptime'] = 100*self.information['simnumber'] + self.information['name_obj'] = 'FLAT' + + self.information['rotangle'] = 0 + + self.information['CRVAL1'] = 0 + self.information['CRVAL2'] = 0 + + self.information['tel_Pra'] = 0 + self.information['tel_Pdec'] = 0 + + self.information['CRPIX1'] = self.information['tel_Pra'] * \ + self.information['pixel_size'] + self.information['CRPIX2'] = self.information['tel_Pdec'] * \ + self.information['pixel_size'] + + self.information['CD1_1'] = - \ + np.cos(theta)*self.information['pixel_size'] + self.information['CD1_2'] = np.sin( + theta)*self.information['pixel_size'] + + self.information['CD2_1'] = np.sin( + theta)*self.information['pixel_size'] + self.information['CD2_2'] = np.cos( + theta)*self.information['pixel_size'] + + self.information['ra_obj'] = 0 + self.information['dec_obj'] = 0 + self.information['ra_pnt0'] = 0 + self.information['dec_pnt0'] = 0 + self.information['target'] = 'FLAT' + + self.image_g += 5000*self.information['simnumber'] + self.image_r += 5000*self.information['simnumber'] + self.image_i += 5000*self.information['simnumber'] + + elif self.source == 'DARK': + + theta = 0 + + self.information['exptime'] = 2*3600*self.information['simnumber'] + self.information['name_obj'] = 'DARK' + + self.information['rotangle'] = 0 + + self.information['CRVAL1'] = 0 + self.information['CRVAL2'] = 0 + + self.information['tel_Pra'] = 0 + self.information['tel_Pdec'] = 0 + self.information['CRPIX1'] = self.information['tel_Pra'] * \ + self.information['pixel_size'] + self.information['CRPIX2'] = self.information['tel_Pdec'] * \ + self.information['pixel_size'] + + self.information['CD1_1'] = - \ + np.cos(theta)*self.information['pixel_size'] + self.information['CD1_2'] = np.sin( + theta)*self.information['pixel_size'] + + self.information['CD2_1'] = np.sin( + theta)*self.information['pixel_size'] + self.information['CD2_2'] = np.cos( + theta)*self.information['pixel_size'] + + self.information['ra_obj'] = 0 + self.information['dec_obj'] = 0 + self.information['ra_pnt0'] = 0 + self.information['dec_pnt0'] = 0 + self.information['target'] = 'DARK' + + self.image_g = np.zeros( + (self.information['ysize'], self.information['xsize']), dtype=float) + self.image_r = np.zeros( + (self.information['ysize'], self.information['xsize']), dtype=float) + self.image_i = np.zeros( + (self.information['ysize'], self.information['xsize']), dtype=float) + + elif self.source == 'BIAS': + + theta = 0 + + self.information['exptime'] = 0.0 + + self.information['name_obj'] = 'BIAS' + + self.information['rotangle'] = 0 + + self.information['CRVAL1'] = 0 + self.information['CRVAL2'] = 0 + + self.information['tel_Pra'] = 0 + self.information['tel_Pdec'] = 0 + + self.information['CRPIX1'] = self.information['tel_Pra'] * \ + self.information['pixel_size'] + self.information['CRPIX2'] = self.information['tel_Pdec'] * \ + self.information['pixel_size'] + + self.information['CD1_1'] = - \ + np.cos(theta)*self.information['pixel_size'] + self.information['CD1_2'] = np.sin( + theta)*self.information['pixel_size'] + + self.information['CD2_1'] = np.sin( + theta)*self.information['pixel_size'] + self.information['CD2_2'] = np.cos( + theta)*self.information['pixel_size'] + + self.image_g = np.zeros( + (self.information['ysize'], self.information['xsize']), dtype=float) + self.image_r = np.zeros( + (self.information['ysize'], self.information['xsize']), dtype=float) + self.image_i = np.zeros( + (self.information['ysize'], self.information['xsize']), dtype=float) + self.information['target'] = 'BIAS' + self.information['ra_obj'] = 0 + self.information['dec_obj'] = 0 + self.information['ra_pnt0'] = 0 + self.information['dec_pnt0'] = 0 + + else: + print('Souce is not correct and programe will return') + sys.exit(1) + + ######################################################################################################################## + ######################################################################################################################## + + # Apply flat-field large scale structure for one chip + if self.flatfieldM: + self.applyFlat() + print('applyFlat') + ######################################################## + if self.PRNUeffect: + self.applyPRNUeffect() + print('applyPRNUeffect') + + #################################################################### + if self.source in sourcelist: + if self.cosmicRays: + + self.addCosmicRays() + print('addCosmicRays finisth') + + ################################################## + if self.skyback: + if self.source in sourcelist: + self.applyPoissonSkyNoise() + print('apply Poisson Sky Noise') + ################################################# + + if self.darknoise: + self.applyDarkCurrent() + print('applyDarkCurrent finisth') + + ################################################# + + if self.shutterEffect and self.information['exptime'] > 2.6: + self.addshuttereffect() + print('add shutter effect ') + + ################################################# + + if self.bleeding: + + self.image_g = self.applyBleeding(self.image_g.copy()) + self.image_r = self.applyBleeding(self.image_r.copy()) + self.image_i = self.applyBleeding(self.image_i.copy()) + print('apply bleeding effect finisth') + + ################################################ + + if self.nonlinearity: + self.applyNonlinearity() + print('applyNonlinearity') + + # ################################################ + # read CCD image to output matrix as defined + + if self.overscans and self.information['xsize'] == 9216 and self.information['ysize'] == 9232: + + self.image_g = self.addScanEffect(self.image_g.copy()) + self.image_r = self.addScanEffect(self.image_r.copy()) + self.image_i = self.addScanEffect(self.image_i.copy()) + print('apply overscan zone 1-16') + + ################################################################# + if self.source in sourcelist: + if self.radiationDamage: + + self.applyRadiationDamage() + print('applyRadiationDamage()') + ###################################################################### + ###### apply readoutNoise ###### + + if self.readoutNoise: + self.applyReadoutNoise() + print('apply readout noise') + + ###################################################################### + if self.information['xsize'] == 9216 and self.information['ysize'] == 9232 and self.overscans: + self.electrons2ADU() # apply gain + print('apply 1-16 zone by different gain ') + + else: + self.image_g = self.image_g/1.5 + self.image_r = self.image_r/1.5 + self.image_i = self.image_i/1.5 + print('apply gain finish') + ###################################################################### + # size of the output image array, xsize is column, ysize is row, xsize = 9216,ysize = 9232 + + if self.information['xsize'] == 9216 and self.information['ysize'] == 9232 and self.overscans: + self.applyBias() + print('apply bias finish') + else: + self.image_g = self.image_g+500.0 + self.image_r = self.image_r+500.0 + self.image_i = self.image_i+500.0 + + ################################################################### + if self.cosmetics and self.information['xsize'] == 9216 and self.information['ysize'] == 9232: + + self.applyCosmetics() + print('apply Cosmetics...') + + ################################################## + if self.intscale: + self.discretise() + ################################################## + self.writeOutputs(simnumber) + + ####################################################### + self.log.info('Using the following input values:') + + for key, value in self.information.items(): + self.log.info('%s = %s' % (key, value)) + self.log.info('Using the following booleans:') + + for key, value in self.booleans.items(): + self.log.info('%s = %s' % (key, value)) + ####################################################### + + self.log.info('Finished the ith_Exposure = %i' % (simnumber)) + + print('The iLoop= % d simlaiton finished. ' % simnumber) + + +################################################################################################ + +def runMCIsim(sourcein, configfile, dir_path, result_path, debug, iLoop): + + print('Path Test：dir_path', dir_path) + + sim = dict() + sim[iLoop] = MCIsimulator(configfile) + + sim[iLoop].configure(iLoop, sourcein, dir_path, + result_path) # load the configfile; + + sim[iLoop].information['sourcein'] = sourcein + sim[iLoop].information['debug'] = debug + sim[iLoop].information['result_path'] = result_path + + sim[iLoop].simulate(iLoop) + + return + +################################################################################################