diff --git a/csst_mci_sim/csst_mci_sim.py b/csst_mci_sim/csst_mci_sim.py deleted file mode 100644 index 4759f1bf89ab51c4112fcb4e5442ef4a3c65392d..0000000000000000000000000000000000000000 --- a/csst_mci_sim/csst_mci_sim.py +++ /dev/null @@ -1,5064 +0,0 @@ -#!/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 - -------- - -""" -import galsim -import pandas as pd -#from scipy.integrate import simps -from datetime import datetime, timedelta -from astropy.time import Time -from astropy.coordinates import get_sun -import numpy as np -from tqdm import tqdm -from astropy.table import Table -from astropy.io import fits -from astropy import units as u -import os, sys, math -import configparser as ConfigParser -#from matplotlib import pyplot as plt - -from scipy import ndimage -###################### -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 -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 - -########################### 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 -""" -import numpy as np - - -#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 -''' - -import scipy.spatial as spatial - -###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 - - wave = np.linspace(2500,10000,8501) # set the wavelenth for MCI from 250nm to 1100nm - 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，望远镜聚光面积�? - fluxlam=1e-17*(flux) # convert original unit to unit of erg/s/A/cm^2 - # wave in nm ;; - ephoton=planckh*cc/wave # single photon energy in erg/photon; cc与lambda单位需要一致； - - 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) - - iwave=self.filter_psf[ch]['psf_iwave'][:] ### seven wavelength for the seleced filter - - 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 - - wave = np.linspace(2500,11000, 8501) # set the wavelenth for MCI from 250nm to 1100nm - - 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, - Nphoton_flux =skynoise_flux/ephoton*telarea*0.05*0.05*self.information['exptime'] # each pixel will have the photons in unit of photons/s/A - - #### 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 ##################### - starcat=self.information['starcat'] #'selection_MACSJ0744.9 3927_20230517_concat.csv' - - - - 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: - psf[ch]=ndimage.rotate(temp, theta.deg, order=1, reshape=False) # here we choose reshape=False, the rotated image will - 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= fsy/pixelscale+self.information['ysize']/2-0.5 ### row number on CCD image - col= fsx/pixelscale+self.information['xsize']/2-0.5 ### col number on CCD image - - #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): - psf[ch]=ndimage.rotate(temp, theta.deg, order=1, reshape=True) # here we choose reshape=False, the rotated image will - 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) - - pos = galsim.PositionD(x=int(cx), y=int(cy)) ### set subimge center as the pos - - 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 - spec_erg2 = spec_erg / 4.25452e10 # erg/s/cm^2/A/arcsec^2 - - self.zodiacal_wave=wave_A # in A - - self.zodiacal_flux=spec_erg2 - - return wave_A, spec_erg2 - ################################################################################## - - ########################################################################## - def cal_Lensing_galaxy_img(self): - ##################### - - - self.dsx_arc = self.information['pixel_size'] # arcsec, pixel scale size; - self.bsz_arc = max(self.information['ysize'] ,self.information['xsize'])*self.dsx_arc # arcsec, field size; - 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) - - galRa = newRa[k1]*3600##*(srcs_cat[k1].header['new_ra'] -self.information['gal_ra'] +self.information['star_ra']) # ra of galaxies, arcsecond - galDec = newDec[k1]*3600####(srcs_cat[k1].header['new_dec']-self.information['gal_dec']+self.information['star_dec']) # dec of galaxies, arcsecond - - ########################################################################### - ############################################################################ - ### 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 - gal_flux=srcs_sed[k1+1].data ## here is k1+1, not k1, k1 begins with 0 - - - - ################################ - ### rotate the lensed_images_g ### - if abs(theta.deg)>0: - lensed_images_g=ndimage.rotate(srcs_cat[k1+1].data, -theta.deg, order=1, reshape=True) # here we choose reshape=False, the rotated image will - 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: - psf[ch]=ndimage.rotate(temp, theta.deg, order=1, reshape=False) # here we choose reshape=False, the rotated image will - 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) - - pos = galsim.PositionD(x=int(cx), y=int(cy)) ### set subimge center as the pos - - 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() - - self.cr['exptime'] = self.information['exptime'] #to scale the number of cosmics with exposure time - - #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]) & (idxt_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 xminnp.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() 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 - - self.prescan =self.information['prescan'] ## horizon direction , columnn - 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 - subimg=img[xmin[k]:xmax[k],ymin[k]:ymax[k]] # get k-th subimg from input img - subimg=np.flipud(subimg) ## down to up, up to down - - - - if k>=5 and k<=8: ### zos5 zos6 zos7 zos8 - subimg=img[xmin[k]:xmax[k],ymin[k]:ymax[k]] # get k-th subimg from input img - 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 - subimg=img[xmin[k]:xmax[k],ymin[k]:ymax[k]] # get k-th subimg from input img - subimg=np.fliplr(subimg) ## left to right ,right to left - - - - if k>=13 and k<=16: ### zos5 zos6 zos7 zos8 - subimg=img[xmin[k]:xmax[k],ymin[k]:ymax[k]] # get k-th subimg from input img - - - ###################### - - 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]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) - - self.TianCe_day=str(self.dt.year)+self.dt.strftime("-%m-%d") ####str(self.dt.year)+'-'+str(self.dt.month)+'-'+str(self.dt.day) - - 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 - rotTelPos = 2*(ud-1)*5.0*galsim.degrees # the image will rotate in clockwise rotaion between -10-10 degree - 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 - - rotTelPos = 0.0*galsim.degrees # the image will rotate in clockwise rotaion between -10-10 degree - rotSkyPos = 0.0*galsim.degrees - - - else: - - dis_ra = 0 # - dis_dec = 0 # T - - rotTelPos = 0.0*galsim.degrees # the image will rotate in clockwise rotaion between -10-10 degree - rotSkyPos = 0.0*galsim.degrees - - ##################################################### - - self.information['T_disRa'] =dis_ra # telescope pointing Ra displacement - self.information['T_disDec'] =dis_dec # telescope pointing Dec displacement - - self.information['rotTelPos']=rotTelPos # telescope rotation in degree - self.information['rotSkyPos']=rotSkyPos # sky rotation in degree - - 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 - self.information['exptime']=100*self.information['simnumber']#np.ceil(self.information['simnumber']/2.0) - 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 - -################################################################################################