diff --git a/.idea/inspectionProfiles/profiles_settings.xml b/.idea/inspectionProfiles/profiles_settings.xml new file mode 100644 index 0000000000000000000000000000000000000000..105ce2da2d6447d11dfe32bfb846c3d5b199fc99 --- /dev/null +++ b/.idea/inspectionProfiles/profiles_settings.xml @@ -0,0 +1,6 @@ + + + + \ No newline at end of file diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000000000000000000000000000000000000..a9b5e744c0d05f2592bc82552323c2044a5bfbce --- /dev/null +++ b/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2022 CSST-L1 + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/README.md b/README.md index f9a5328c3c39539219a075366937c977d253077e..aed373c9603307c3bddd511cc45de92fabdf830a 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -# Ifs +# csst_ifs_common @@ -15,14 +15,14 @@ Already a pro? Just edit this README.md and make it your own. Want to make it ea ``` cd existing_repo -git remote add origin https://csst-tb.bao.ac.cn/code/shaosim/ifs.git +git remote add origin https://csst-tb.bao.ac.cn/code/csst-l1/ifs/csst_ifs_common.git git branch -M main git push -uf origin main ``` ## Integrate with your tools -- [ ] [Set up project integrations](https://csst-tb.bao.ac.cn/code/shaosim/ifs/-/settings/integrations) +- [ ] [Set up project integrations](http://10.3.10.28/code/csst-l1/ifs/csst_ifs_common/-/settings/integrations) ## Collaborate with your team diff --git a/csst_ifs_sim/CTI/CTI.py b/csst_ifs_sim/CTI/CTI.py new file mode 100644 index 0000000000000000000000000000000000000000..48c87cb8aee4207a0b4d764e32a0cea44a3d6987 --- /dev/null +++ b/csst_ifs_sim/CTI/CTI.py @@ -0,0 +1,495 @@ +""" +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 + +# try: +# import cdm03bidir +# #import cdm03bidirTest as cdm03bidir #for testing purposes only +# except ImportError: +# print('import CTI module') +# #print ('No CDM03bidir module available, please compile it: f2py -c -m cdm03bidir cdm03bidir.f90') + + + +#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=1.) + #update with inputs + self.params.update(self.values) + + #read in trap information + trapdata = np.loadtxt(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['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 radiateFullCCD(self): + """ + This routine allows the whole CCD to be run through a radiation damage mode. + The routine takes into account the fact that the amplifiers are in the corners + of the CCD. The routine assumes that the CCD is using four amplifiers. + + There is an excess of .copy() calls, which should probably be cleaned up. However, + given that I had problem with the Fortran code, I have kept the calls. If memory + becomes an issue then this should be cleaned. + + :return: radiation damaged image + :rtype: ndarray + """ + ydim, xdim = self.data.shape + out = np.zeros((xdim, ydim)) + + #transpose the data, because Python has different convention than Fortran + data = self.data.transpose().copy() + + for quad in self.values['quads']: + if self.logger: + self.log.info('Adding CTI to Q%i' % quad) + + if quad == 0: + d = data[0:self.values['xsize'], 0:self.values['ysize']].copy() + tmp = self.applyRadiationDamage(d, iquadrant=quad).copy() + out[0:self.values['xsize'], 0:self.values['ysize']] = tmp + elif quad == 1: + d = data[self.values['xsize']:, :self.values['ysize']].copy() + tmp = self.applyRadiationDamage(d, iquadrant=quad).copy() + out[self.values['xsize']:, :self.values['ysize']] = tmp + elif quad == 2: + d = data[:self.values['xsize'], self.values['ysize']:].copy() + tmp = self.applyRadiationDamage(d, iquadrant=quad).copy() + out[:self.values['xsize'], self.values['ysize']:] = tmp + elif quad == 3: + d = data[self.values['xsize']:, self.values['ysize']:].copy() + tmp = self.applyRadiationDamage(d, iquadrant=quad).copy() + out[self.values['xsize']:, self.values['ysize']:] = tmp + else: + print( 'ERROR -- too many quadrants!!' ) + self.log.error('Too many quadrants! This method allows only four quadrants.') + + return out.transpose() + + + def radiateFullCCD2(self): + """ + This routine allows the whole CCD to be run through a radiation damage mode. + The routine takes into account the fact that the amplifiers are in the corners + of the CCD. The routine assumes that the CCD is using four amplifiers. + + There is an excess of .copy() calls, which should probably be cleaned up. However, + given that I had problem with the Fortran code, I have kept the calls. If memory + becomes an issue then this should be cleaned. + + :return: radiation damaged image + :rtype: ndarray + """ + ydim, xdim = self.data.shape + out = np.empty((ydim, xdim)) + + #transpose the data, because Python has different convention than Fortran + data = self.data.copy() + + for quad in self.values['quads']: + if self.logger: + self.log.info('Adding CTI to Q%i' % quad) + + if quad == 0: + d = data[:self.values['ysize'], :self.values['xsize']].copy() + tmp = self.applyRadiationDamage(d, iquadrant=quad).copy() + out[:self.values['ysize'], :self.values['xsize']] = tmp + elif quad == 1: + d = data[:self.values['ysize'], self.values['xsize']:].copy() + tmp = self.applyRadiationDamage(d, iquadrant=quad).copy() + out[:self.values['ysize'], self.values['xsize']:] = tmp + elif quad == 2: + d = data[self.values['ysize']:, :self.values['xsize']].copy() + tmp = self.applyRadiationDamage(d, iquadrant=quad).copy() + out[self.values['ysize']:, :self.values['xsize']] = tmp + elif quad == 3: + d = data[self.values['ysize']:, self.values['xsize']:].copy() + tmp = self.applyRadiationDamage(d, iquadrant=quad).copy() + out[self.values['ysize']:, self.values['xsize']:] = tmp + else: + print( 'ERROR -- too many quadrants!!') + self.log.error('Too many quadrants! This method allows only four quadrants.') + + return out + + + 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) + +################################################################################# +###modify + import sys + #sys.path.append('../so') + from ifs_so import cdm03bidir + # from ifs_so.cdm03.cpython-38-x86_64-linux-gnu import cdm03bidir + # import cdm03bidir + 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) + +################################################################################################################# + + + +class CDM03(): + """ + Class to run CDM03 CTI model, class Fortran routine to perform the actual CDM03 calculations. + + :param data: input data to be radiated + :type data: ndarray + :param input: input parameters + :type input: dictionary + :param log: instance to Python logging + :type log: logging instance + """ + def __init__(self, input, data, log=None): + """ + Class constructor. + + :param data: input data to be radiated + :type data: ndarray + :param input: input parameters + :type input: dictionary + :param log: instance to Python logging + :type log: logging instance + """ + try: + import cdm03 + except ImportError: + print( 'No CDM03 module available, please compile it: f2py -c -m cdm03 cdm03.f90') + + self.data = data + self.values = dict(quads=(0,1,2,3), xsize=2048, ysize=2066, dob=0.0, rdose=8.0e9) + self.values.update(input) + self.log = log + self._setupLogger() + + + def _setupLogger(self): + """ + Set up the logger. + """ + self.logger = True + if self.log is None: + self.logger = False + + + def radiateFullCCD(self): + """ + This routine allows the whole CCD to be run through a radiation damage mode. + The routine takes into account the fact that the amplifiers are in the corners + of the CCD. The routine assumes that the CCD is using four amplifiers. + + There is an excess of .copy() calls, which should probably be cleaned up. However, + given that I had problem with the Fortran code, I have kept the calls. If memory + becomes an issue then this should be cleaned. + + :return: radiation damaged image + :rtype: ndarray + """ + ydim, xdim = self.data.shape + out = np.zeros((xdim, ydim)) + + #transpose the data, because Python has different convention than Fortran + data = self.data.transpose().copy() + + for quad in self.values['quads']: + if self.logger: + self.log.info('Adding CTI to Q%i' % quad) + + if quad == 0: + d = data[0:self.values['xsize'], 0:self.values['ysize']].copy() + tmp = self.applyRadiationDamage(d, iquadrant=quad).copy() + out[0:self.values['xsize'], 0:self.values['ysize']] = tmp + elif quad == 1: + d = data[self.values['xsize']:, :self.values['ysize']].copy() + tmp = self.applyRadiationDamage(d, iquadrant=quad).copy() + out[self.values['xsize']:, :self.values['ysize']] = tmp + elif quad == 2: + d = data[:self.values['xsize'], self.values['ysize']:].copy() + tmp = self.applyRadiationDamage(d, iquadrant=quad).copy() + out[:self.values['xsize'], self.values['ysize']:] = tmp + elif quad == 3: + d = data[self.values['xsize']:, self.values['ysize']:].copy() + tmp = self.applyRadiationDamage(d, iquadrant=quad).copy() + out[self.values['xsize']:, self.values['ysize']:] = tmp + else: + print ('ERROR -- too many quadrants!!') + self.log.error('Too many quadrants! This method allows only four quadrants.') + + return out.transpose() + + + def radiateFullCCD2(self): + """ + This routine allows the whole CCD to be run through a radiation damage mode. + The routine takes into account the fact that the amplifiers are in the corners + of the CCD. The routine assumes that the CCD is using four amplifiers. + + There is an excess of .copy() calls, which should probably be cleaned up. However, + given that I had problem with the Fortran code, I have kept the calls. If memory + becomes an issue then this should be cleaned. + + :return: radiation damaged image + :rtype: ndarray + """ + ydim, xdim = self.data.shape + out = np.empty((ydim, xdim)) + + #transpose the data, because Python has different convention than Fortran + data = self.data.copy() + + for quad in self.values['quads']: + if self.logger: + self.log.info('Adding CTI to Q%i' % quad) + + if quad == 0: + d = data[:self.values['ysize'], :self.values['xsize']].copy() + tmp = self.applyRadiationDamage(d, iquadrant=quad).copy() + out[:self.values['ysize'], :self.values['xsize']] = tmp + elif quad == 1: + d = data[:self.values['ysize'], self.values['xsize']:].copy() + tmp = self.applyRadiationDamage(d, iquadrant=quad).copy() + out[:self.values['ysize'], self.values['xsize']:] = tmp + elif quad == 2: + d = data[self.values['ysize']:, :self.values['xsize']].copy() + tmp = self.applyRadiationDamage(d, iquadrant=quad).copy() + out[self.values['ysize']:, :self.values['xsize']] = tmp + elif quad == 3: + d = data[self.values['ysize']:, self.values['xsize']:].copy() + tmp = self.applyRadiationDamage(d, iquadrant=quad).copy() + out[self.values['ysize']:, self.values['xsize']:] = tmp + else: + print ('ERROR -- too many quadrants!!') + self.log.error('Too many quadrants! This method allows only four quadrants.') + + return out + + + 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 + """ + #read in trap information + trapdata = np.loadtxt(self.values['trapfile']) + nt = trapdata[:, 0] + sigma = trapdata[:, 1] + taur = trapdata[:, 2] + + iflip = iquadrant / 2 + jflip = iquadrant % 2 + + if self.logger: + self.log.info('nt=' + str(nt)) + self.log.info('sigma= ' + str(sigma)) + self.log.info('taur= ' + str(taur)) + 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) + + + # #call Fortran routine + # CTIed = cdm03.cdm03(np.asfortranarray(data), + # iflip, jflip, + # self.values['dob'], self.values['rdose'], + # nt, sigma, taur) + ###modify + import sys + sys.path.append('../CTI') + import cdm03 + + ################################################################################# + + CTIed = cdm03.cdm03(np.asfortranarray(data), + jflip, iflip, + self.values['dob'], self.values['rdose'], + nt,sigma,taur ) + + + return np.asanyarray(CTIed) + + +################################################################################################################# + diff --git a/csst_ifs_sim/CTI/__init__.py b/csst_ifs_sim/CTI/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/csst_ifs_sim/__init__.py b/csst_ifs_sim/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/csst_ifs_sim/csst_ifs_sim.py b/csst_ifs_sim/csst_ifs_sim.py new file mode 100644 index 0000000000000000000000000000000000000000..1dbd918a2cdaf4cbe59781b2edfd439be1e427a8 --- /dev/null +++ b/csst_ifs_sim/csst_ifs_sim.py @@ -0,0 +1,2624 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" + +@author: yan +""" + +""" +The CSST IFS Image Simulator +============================================= + +This file contains an image simulator for the CSST IFS. + +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, oversampling factor, 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). + + #. Load the wavefront aberration data used to calculate PSF with defined wavelength and field of view. + + #. Apply calibration unit flux to mimic flat field exposures [optional]. + #. 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 and background light from Zodiacal light [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 from a Gaussian distribution [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. + +:version: 2023.04 + + +Future Work +----------- + +.. todo:: + + #. + #. + #. ...... + +Contact Information: zhaojunyan@shao.ac.cn + +------- + +""" + +from scipy.integrate import simps +import pandas as pd +from datetime import datetime, timedelta +import julian +from astropy.coordinates import get_sun +from astropy.time import Time +from astropy import units as u + +from astropy.coordinates import SkyCoord, EarthLocation + +from scipy import interpolate + +import numpy as np + +import scipy.io as sio + +from astropy.io import fits + +from scipy.signal import fftconvolve + +import galsim + +from scipy import ndimage + +import cmath + +import os, sys + +import configparser as ConfigParser + +from optparse import OptionParser + +sys.path.append('../') + +from CTI import CTI +from support import logger as lg +from support import cosmicrays +from support import IFSinstrumentModel + +############################################################################### + +def str2time(strTime): + if len(strTime)>20:# + msec=int(float('0.'+strTime[20:])*1000000) # + str2=strTime[0:19]+' '+str(msec) + return datetime.strptime(str2,'%Y %m %d %H %M %S %f') + +#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(ra, dec, rou=False): + '''convert deg to ra's HMS or dec's DHS''' + RA, DEC, rs, ds = '000000', '000000', '', '' + if dec: + if str(dec)[0] == '-': + ds, dec = '-', abs(dec) + deg = int(dec) + decM = abs(int((dec-deg)*60)) + + if rou: + decS = round((abs((dec-deg)*60)-decM)*60,1) + else: + decS = int((abs((dec-deg)*60)-decM)*60) + + if deg ==0: + deg ="00" + elif deg <10: + deg = "0%s"%deg + + if decM ==0: + decM ="00" + elif decM <10: + decM = "0%s"%decM + + if decS ==0: + decS ="00" + elif decS <10: + decS = "0%s"%decS + + DEC = '{0}{1}{2}{3}'.format(ds, deg, decM, decS) + + if ra: + if str(ra)[0] == '-': + rs, ra = '-', abs(ra) + raH = int(ra/15) + raM = int(((ra/15)-raH)*60) + + if rou: + raS = round(((((ra/15)-raH)*60)-raM)*60,1) + else: + raS = int(((((ra/15)-raH)*60)-raM)*60) + + if raH ==0: + raH = "00" + + elif raH <10: + raH = "0%s"%raH + + if raM ==0: + raM = "00" + elif raM <10: + raM = "0%s"%raM + + if raS ==0: + raS = "00" + + elif raS <10: + raS = "0%s"%raS + + RA = '{0}{1}{2}{3}'.format(rs, raH, raM, raS) + + if ds=='-': + return RA+DEC + else: + return RA+'+'+DEC +############################################################################### + +########################################################### +def beta_angle( x_sat, y_sat, z_sat, vx_sat, vy_sat, vz_sat, ra_obj, dec_obj): + + # get the vector for next motion + x_sat2 = x_sat + vx_sat + y_sat2 = y_sat + vy_sat + z_sat2 = z_sat + vz_sat + + vector1 = np.stack((x_sat, y_sat, z_sat), axis=0) + vector2 = np.stack((x_sat2, y_sat2, z_sat2), axis=0) + vector_normal = np.cross(vector1, vector2) + + location_normal = EarthLocation.from_geocentric(vector_normal[0], + vector_normal[1], + vector_normal[2], + unit=u.km) + radec_normal = SkyCoord(ra=location_normal.lon, + dec=location_normal.lat, + frame='gcrs') + lb_normal = radec_normal.transform_to('geocentrictrueecliptic') + + + radec_sun = SkyCoord(ra=ra_obj*u.degree, dec=dec_obj*u.degree, frame='gcrs') + lb_sun = radec_sun.transform_to('geocentrictrueecliptic') + + # get the angle between normal and solar + angle = 90 - lb_normal.separation(lb_sun).degree + + return angle + +########################################################## +########################################################## +def LSR_velocity(ra,dec,velocity,Obstime): + + local = EarthLocation.from_geodetic(lat=22.349368*u.deg, lon=113.584068*u.deg, height=10*u.m) + ### convert ra and dec to + source = SkyCoord(ra*u.deg, dec*u.deg,frame='icrs', unit=(u.hourangle,u.deg) ) + l = source.galactic.l.deg + b = source.galactic.b.deg + c = SkyCoord(l=l*u.degree, b=b*u.degree, frame='galactic') + c_icrs = c.transform_to('icrs') + barycorr = c_icrs.radial_velocity_correction(obstime=Time(Obstime), location=local) + velocity = velocity + barycorr.value/1000 + #print(barycorr.value/1000) + l = l * np.pi / 180 + b = b * np.pi / 180 + return velocity + 9 * np.cos(l) * np.cos(b) + 12 * np.sin(l) * np.cos(b) + 7 * np.sin(b) + +############################################################################### +def rotation_yan(x0,y0,x1,y1,angle): + #% point A (x0,y0) + #% point B(x1,y1) + #% roate angle ,point B rotate with point A + alpha=angle/180*3.1415926 #% in radians + x2=(x1-x0)*np.cos(alpha)-(y1-y0)*np.sin(alpha)+x0 + y2=(x1-x0)*np.sin(alpha)+(y1-y0)*np.cos(alpha)+y0 + return x2, y2 +###################################################################### +def centroid(data): + h,w = np.shape(data) + x = np.arange(0,w) + y = np.arange(0,h) + x1=np.ones((1,h)) + y1=np.ones((w,1)) + cx=(np.dot(np.dot(x1,data),y)) /(np.dot(np.dot(x1,data),y1)) + cy=(np.dot(np.dot(x ,data),y1)) /(np.dot(np.dot(x1,data),y1)) + return np.float64(cx),np.float64(cy) +#################################################################### +def SpecCube2photon(data,wave): + # calcutle photons from original spec cube data, + # data: input data, two-dimentional, unit : 1e-17 erg/s/A/cm^2 + # wave: the relative wavefront, unit in nm; + planckh= 6.62620*10**-27 # % erg s; + cc=2.99792458*10**17 # in nm/s + # fov2=0.01 # in arcsec^2 + # telarea=3.1415926*100*100 # in cm^2， + fluxlam=1e-17*data # convert original unit to unit of erg/s/A/cm^2 + lam=wave # wave in nm ;; + ephoton=planckh*cc/lam # in erg/photon, + Nphoton =fluxlam/ephoton # in unit of photons/cm2/s/A + return Nphoton + +###################################################################### + +###################################################################### +def zero_pad(image, shape, position='corner'): + """ + Extends image to a certain size with zeros + Parameters + ---------- + image: real 2d `np.ndarray` + Input image + shape: tuple of int + Desired output shape of the image + position : str, optional + The position of the input image in the output one: + * 'corner' + top-left corner (default) + * 'center' + centered + Returns + ------- + padded_img: real `np.ndarray` + The zero-padded image + """ + shape = np.asarray(shape, dtype=int) + imshape = np.asarray(image.shape, dtype=int) + + if np.alltrue(imshape == shape): + return image + + if np.any(shape <= 0): + raise ValueError("ZERO_PAD: null or negative shape given") + + dshape = shape - imshape + if np.any(dshape < 0): + raise ValueError("ZERO_PAD: target size smaller than source one") + + pad_img = np.zeros(shape, dtype=image.dtype) + + idx, idy = np.indices(imshape) + + if position == 'center': + if np.any(dshape % 2 != 0): + raise ValueError("ZERO_PAD: source and target shapes " + "have different parity.") + offx, offy = dshape // 2 + else: + offx, offy = (0, 0) + + pad_img[idx + offx, idy + offy] = image + + return pad_img + +def anySampledPSF(wavefront, pupil , Q , sizeout): + ''' + written on 2020.12.04 by Yan @Shao + % wavefront sampled in the united circle; + % pupil function samped as wavefront; + % sample ratio Q=lambda/D/pixel; + % lambda is wavelength; + % D is diameter; + % pixel is in arcsec; + % or sample ratio Q=lambda*f/D/pixelsize + % f is focal length; + % pixelsize is the actural size of the detector; + % make sure all the varia have the same unit; + % the returned PSF has sum value of 1 + ''' + + m,n = np.shape(wavefront) + + phase=2.0*np.pi*wavefront #% the phase of the input wavefront; + + ####generalized pupil function of channel1; + if Q>=1: + pk1=zero_pad(pupil*np.exp(cmath.sqrt(-1)*(phase)),( int(np.floor(m*Q)),int(np.floor(m*Q))), position='corner') + psf=abs(np.fft.fftshift(np.fft.fft2(pk1)))**2 + #psf=psf/psf.sum() + else: #%%% case: Q<1 + # %fft method + if Q>=0.5: #% in this case Q<1 and Q>=0.5 + pk1=zero_pad(pupil*np.exp(cmath.sqrt(-1)*(phase)), (int(2*np.floor(m*Q)),int(2*np.floor(m*Q))), position='corner') + mypsf=np.fft.fft2(pk1); + t=mypsf[0::2,0::2]; + psf=abs(np.fft.fftshift(t))**2 + #psf=mypsf3/mypsf3.sum(); + + else: + if Q>=0.25: + pk1=zero_pad(pupil*np.exp(cmath.sqrt(-1)*(phase)),( int(4*np.floor(m*Q)),int(4*np.floor(m*Q))), position='corner') + mypsf=np.fft.fft2(pk1); + t=mypsf[0::4,0::4]; + psf=abs(np.fft.fftshift(t))**2 + #psf=mypsf3/mypsf3.sum(); + else: + print('---- Q<0.25 , stop running') + sys.exit(0) + ######################################################################################## + mm,nn=np.shape(psf) + if np.mod(sizeout,2)==0: + Nx=sizeout/2-1 + else: + Nx=(sizeout+1)/2-1 + + + if max(mm,nn)3*sigma: + + SpmOut=1.0/20000 + else: + + column=np.where(d==min(d)) + SpmOut=Qt[column[0]]*np.exp(-(lam-wave[column[0]])**2/2/sigma**2) + + return SpmOut + +##################################################################################### +#################################################################################################################### + +def processArgs(printHelp=False): + """ + Processes command line arguments. + """ + parser = OptionParser() + + parser.add_option('-c', '--configfile', dest='configfile', + help="Name of the configuration file", metavar="string") + parser.add_option('-s', '--section', dest='section', + help="Name of the section of the config file [SCIENCE]", metavar="string") + parser.add_option('-q', '--quadrant', dest='quadrant', help='CCD quadrant to simulate [0, 1, 2, 3]', + metavar='int') + parser.add_option('-x', '--xCCD', dest='xCCD', help='CCD number in X-direction within the FPA matrix', + metavar='int') + parser.add_option('-y', '--yCCD', dest='yCCD', help='CCD number in Y-direction within the FPA matrix', + metavar='int') + parser.add_option('-d', '--debug', dest='debug', action='store_true', + help='Debugging mode on') + parser.add_option('-t', '--test', dest='test', action='store_true', + help='Run unittest') + parser.add_option('-f', '--fixed', dest='fixed', help='Use a fixed seed for the random number generators', + metavar='int') + if printHelp: + parser.print_help() + else: + return parser.parse_args() +############################################################################################## +class IFSsimulator(): + """ + CSST IFS Simulator + + The image that is being build is in:: + + self.image + self.image_b blue channel + self.image_r red channel + + :param opts: OptionParser instance + :type opts: OptionParser instance + """ + + def __init__(self, opts): + """ + Class Constructor. + + :param opts: OptionParser instance + :type opts: OptionParser instance + + """ + + #################################### + + self.configfile = opts.configfile + + if opts.section is None: + ####self.section = 'DEFAULT' + self.section = 'TEST' #####simulation section; + + else: + self.section = opts.section + + #load instrument model, these values are also stored in the FITS header + self.information = IFSinstrumentModel.IFSinformation() + + #update settings with defaults + self.information.update(dict(quadrant=int(0), + ccdx=int(0), + ccdy=int(0), + psfoversampling=1.0, + ccdxgap=1.643, + ccdygap=8.116, + fullwellcapacity=90000, + readouttime=4, + rdose=8.0e9, + coveringfraction=0.5 )) + + ###################################################################### + self.configure() #print the configfile name and path; + + self.information.update(dict( + + cosmicraylengths =self.information['indata_path']+'/cdf_cr_length.dat', + cosmicraydistance=self.information['indata_path']+'/cdf_cr_total.dat', + parallelTrapfile =self.information['indata_path']+'/cdm_euclid_parallel.dat', + serialTrapfile =self.information['indata_path']+'/cdm_euclid_serial.dat', + cosmeticsfile_b =self.information['indata_path']+'/Cosmetics_b.txt', + cosmeticsfile_r =self.information['indata_path']+'/Cosmetics_r.txt' )) + + # public system information + self.pixelscale=0.1 # the pixel size is 0.1 arcsec + self.pixelsize=15 # the pixel physical size is 15 micron + self.Fnum=15.469875 # the f number= focal_length/D; + self.telD=2 # tht telescope size is 2 meter, diamter size; + self.telarea=3.1415926*(self.telD/2)**2*100*100 # telescope square in cm^2 + self.fov2=0.01 # the pixel square + self.planckh= 6.62620*1e-27 # % erg s; + self.cc=2.99792458*1e17 # in nm/s + self.light_FWHM=0.175 + self.HgArsigma=self.light_FWHM/2.35; ## sigma value of the light source; + + #### load system optical and CCD efficiency data + matfn0=self.information['indata_path']+'/TotalQ200923.mat' + da0=sio.loadmat(matfn0) + self.optical_blue_Q=da0['opticalQ_1'] # optical efficiency of blue channel + self.optical_red_Q =da0['opticalQ_2'] # optical efficiency of red channel + self.CCD_Qe=da0['ccdQE'] # red channel + + ############################################################# load all useful data; + ############### load wavefront data; + matfn2=self.information['indata_path']+'/opd_638nm.mat' + da2=sio.loadmat(matfn2) + opd0=da2['opd'] # opd unit is in meter + self.opd0=opd0*1e9 # convert opd to nm + self.pupil=abs(opd0)>0.0 + #################### + + da=fits.open(self.information['indata_path']+'/opd1.fits') + self.opd1=da[0].data + + da=fits.open(self.information['indata_path']+'/opd2.fits') + self.opd2=da[0].data + + da=fits.open(self.information['indata_path']+'/opd3.fits') + self.opd3=da[0].data + + + ################################################################################ + ## slice information; the slice is put in vertial direction; + # the silce long size is 64 pixels + slice_blue=dict() + slice_red =dict() + + randRedpos=np.array([0.52835362, 1.1105926 , -0.81667794, 0.88860884, -2.78092636, + -0.15810022, -1.56726852, -0.71601855, -1.31768647, 1.73107581, + 0.4933876 , 2.83673377, 0.22226286, -0.02634998, 0.35539383, + -0.91989574, -2.44856212, 0.91020484, -3.03097852, -1.11638071, + 1.28360669, -0.12521128, -0.3907698 , 0.70183575, 1.00578099, + 1.67339662, 0.18067182, -0.56303075, 0.40959616, 1.45518379, + -0.93194744, 0.41492972]) + + randBluepos=np.array([ 0.97449008, -0.21371406, -1.62513338, -3.06938604, 1.72615283, + 0.73333374, 0.80923919, -0.9418576 , -0.16806578, -1.04416631, + 2.20155068, -0.0900156 , 0.07597706, 0.76208373, 0.29426592, + -0.89434703, 0.34017826, 1.16936499, 0.10977829, -1.31179539, + -0.50859787, -1.01891651, -0.95791581, -1.53018041, 0.88358827, + 0.07837641, -0.86157157, -1.18070438, 0.53970599, 1.4913733 , + 2.10938775, 1.23213412]) + + ##################### + for i in range(32): + if i==0: + slice_blue['py']=np.zeros(32) + slice_blue['px']=np.zeros(32) + + slice_red['py']=np.zeros(32) + slice_red['px']=np.zeros(32) + + if i<16: + slice_blue['py'][i]=50+randBluepos[i]*4 + slice_blue['px'][i]=3.55/0.015*i+166.0 + + + slice_red['py'][i]=50+randBluepos[i]*4 + slice_red['px'][i]=3.55/0.015*i+1190.0 + + else: + slice_blue['py'][i]=50+250+randBluepos[i]*4 + slice_blue['px'][i]=3.55/0.015*(i-16)+166.0+118 + + + slice_red['py'][i]=50+250+randBluepos[i]*4 + slice_red['px'][i]=3.55/0.015*(i-16)+1190.0+118 + ####### + + self.slice_blue=slice_blue + self.slice_red=slice_red + + ############################################################################### + maskSlice=dict() + maskSlit=dict() + sizeout=100 + for k in range(32): + maskSlice[str(k)]=np.zeros((sizeout,sizeout)) + maskSlice[str(k)][2*k+18:2*k+20, int(sizeout/2)-32:int(sizeout/2)+32]=1 + + maskSlit[str(k)]=np.zeros((sizeout,sizeout)) + maskSlit[str(k)][2*k+18:2*k+20, int(sizeout/2)-37:int(sizeout/2)+37]=1 + + self.maskSlice=maskSlice + self.maskSlit =maskSlit + ################################################################################################ + + + + +############################################################################3 + def readConfigs(self): + """ + Reads the config file information using configParser and sets up a logger. + """ + self.config = ConfigParser.RawConfigParser() + + self.config.read_file(open(self.configfile)) + + +############################################################################### + + + def processConfigs(self): + """ + Processes configuration information and save the information to a dictionary self.information. + + The configuration file may look as follows:: + + [TEST] + + + 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.IFSinstrumentModel.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) + + 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.skyback = self.config.getboolean(self.section, 'skyback') + + + #these don't need to be in the config file + + try: + self.nonlinearity = self.config.getboolean(self.section, 'nonlinearity') + except: + self.nonlinearity = False + + try: + self.flatfieldM = self.config.getboolean(self.section, 'flatfieldM') + except: + self.flatfieldM = False + + try: + self.readoutNoise = self.config.getboolean(self.section, 'readoutNoise') + except: + self.readoutNoise = True + + #################################################################### + self.booleans = dict(nonlinearity=self.nonlinearity, + flatfieldM=self.flatfieldM, + cosmicRays=self.cosmicRays, + darknoise=self.darknoise, + cosmetics=self.cosmetics, + radiationDamage=self.radiationDamage, + bleeding=self.bleeding, + skyback =self.skyback) + ##################################################################### + + now=datetime.now() + + result_day=now.strftime("%Y-%m-%d") + + self.result_path=self.information['result_path']+'/'+result_day + + if os.path.isdir(self.result_path)==False: + os.mkdir(self.result_path) + os.mkdir(self.result_path+'/log_file') + os.mkdir(self.result_path+'/sky_Data') + + now=datetime.now() + + data_time=now.strftime("%Y-%m-%d-%H-%M-%S") + + self.log = lg.setUpLogger(self.result_path+'/log_file/IFS_'+'_'+data_time+'.log') + + self.log.info('STARTING A NEW SIMULATION') + + self.log.info(self.information) + + 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['cosmicraylengths']) + crDists = np.loadtxt(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]) + +############################################################################################## + + + def configure(self): + """ + Configures the simulator with input information and creates and empty array to which the final image will + be build on. + """ + self.readConfigs() + + self.processConfigs() + + self.log.info('Read in the configuration files and created an empty array') + +################################################################################################################# + + 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 addLampFlux(self): + """ + Include flux from the calibration source. + """ + + + self.image_b += fits.getdata(self.information['flatflux']) + + self.image_r += fits.getdata(self.information['flatflux']) + + self.log.info('Flux from the calibration unit included (%s)' % self.information['flatflux']) + + + def applyflatfield(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. + """ + + ### + flat_b=self.MakeFlatMatrix(self.image_b, 100) + + flat_r=self.MakeFlatMatrix(self.image_r, 200) + + self.image_b *= flat_b + + self.image_r *= flat_r + + self.log.info('Applied flatfield to images.') + + return + + ######################################################## + +############################################################################### + 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_b = np.zeros((2048, 4096), dtype=np.float64) + + crImage_r = np.zeros((3072, 6144), dtype=np.float64) + + #cosmic ray instance + cosmics_b = cosmicrays.cosmicrays(self.log, crImage_b, crInfo=self.cr) + cosmics_r = cosmicrays.cosmicrays(self.log, crImage_r, crInfo=self.cr) + + + #add cosmic rays up to the covering fraction + #CCD_cr = cosmics.addUpToFraction(self.information['coveringFraction'], limit=None) + + CCD_cr_b = cosmics_b.addUpToFraction(self.information['coveringfraction'], limit=None) + CCD_cr_r = cosmics_r.addUpToFraction(self.information['coveringfraction'], limit=None) + + + #paste the information + self.image_b += CCD_cr_b + self.image_r += CCD_cr_r + + + #count the covering factor + area_cr_b = np.count_nonzero(CCD_cr_b) + area_cr_r = np.count_nonzero(CCD_cr_r) + + #self.log.info('The cosmic ray covering factor is %i pixels ' % area_cr) + self.log.info('The cosmic ray in blue channel covering factor is %i pixels ' % area_cr_b) + self.log.info('The cosmic ray in red channel covering factor is %i pixels ' % area_cr_r) + + ######################################################### + +######################################################################### + +######################################################################### + + + def applyDarkCurrent(self): + """ + Apply dark current. Scales the dark with the exposure time. + + Additionally saves the image without noise to a FITS file. + """ + + self.log.info('Added dark current to bule and red channel' ) + + ########## blue zone 1 + self.image_b[0:1024,0:2048] += self.information['exptime'] * self.information['dark1_b'] + + ########## zone 4 ################# + self.image_b[1024:2048,0:2048] += self.information['exptime'] * self.information['dark4_b'] + + ########## zone 2 ################### + self.image_b[0:1024,2048:4096] += self.information['exptime'] * self.information['dark2_b'] + + ########## zone 3 + self.image_b[1024:2048,2048:4096]+= self.information['exptime'] * self.information['dark3_b'] + + + + ########## red zone 1 + self.image_r[0:1536, 0:3072] += self.information['exptime'] * self.information['dark1_r'] + ########## zone 4 ################# + + self.image_r[1536:3712,0:3072] += self.information['exptime'] * self.information['dark4_r'] + + ########## zone 2 ################### + self.image_r[0:1536,3072:6144] += self.information['exptime'] * self.information['dark2_r'] + + ########## zone 3 + self.image_r[1536:3072,3072:6144] += self.information['exptime'] * self.information['dark3_r'] + + +#######################################################################################################3 + + def applyCosmicBackground(self): + """ + Apply dark the cosmic background. Scales the background with the exposure time. + + Additionally saves the image without noise to a FITS file. + """ + + #add background + bcgr = self.information['exptime'] * self.information['cosmic_bkgd'] + #self.image += bcgr + self.image_b += bcgr + self.image_r += bcgr + + self.log.info('Added cosmic background = %f' % bcgr) + + if self.cosmicRays: + #self.imagenoCR += bcgr + self.imagenoCR_b += bcgr + self.imagenoCR_r += bcgr + + ########################################################################################## + + + +############################################################################## + def applyPoissonNoise(self): + """ + Add Poisson noise to the image. + """ + + rounded = np.rint(self.image_b) ### round to + residual = self.image_b.copy() - rounded #ugly workaround for multiple rounding operations... + rounded[rounded < 0.0] = 0.0 + + np.random.seed() + self.image_b = np.random.poisson(rounded).astype(np.float64) + self.log.info('Added Poisson noise on channel blue') + self.image_b += residual + + rounded = np.rint(self.image_r) ### round to + residual = self.image_r.copy() - rounded #ugly workaround for multiple rounding operations... + rounded[rounded < 0.0] = 0.0 + + np.random.seed() + self.image_r = np.random.poisson(rounded).astype(np.float64) + self.log.info('Added Poisson noise on channel red') + self.image_r += residual + + +################################################################################################################### + + def applyCosmetics(self): + """ + Apply cosmetic defects described in the input file. + + Warning:: This method does not work if the input file has exactly one line. + """ + cosmetics = np.loadtxt(self.information['cosmeticsfile_b']) + + x = np.round(cosmetics[:, 0]).astype(int) + y = np.round(cosmetics[:, 1]).astype(int) + value = cosmetics[:, 2] + + cosmetics_b=np.zeros((3712,6784)) + cosmetics_r=np.zeros((3712,6784)) + self.log.info('Adding cosmetic defects to blue channel:' ) + for xc, yc, val in zip(x, y, value): + if 0 <= xc <= 6784 and 0 <= yc <= 3712: + #self.image[yc, xc] = val + self.image_b[yc, xc] = val + cosmetics_b[yc,xc]=val + self.log.info('x=%i, y=%i, value=%f' % (xc, yc, val)) +###################################################################################################### + + cosmetics = np.loadtxt(self.information['cosmeticsfile_r']) + + x = np.round(cosmetics[:, 0]).astype(int) + y = np.round(cosmetics[:, 1]).astype(int) + value = cosmetics[:, 2] + + + self.log.info('Adding cosmetic defects to red channel:' ) + + for xc, yc, val in zip(x, y, value): + if 0 <= xc <= 6784 and 0 <= yc <= 3712: + #self.image[yc, xc] = val + self.image_r[yc, xc] = val + cosmetics_r[yc,xc]=val + self.log.info('x=%i, y=%i, value=%f' % (xc, yc, val)) + + +################################################################################ + + + def applyRadiationDamage(self): + """ + Applies CDM03 radiation model to the image being constructed. + + .. seealso:: Class :`CDM03` + """ + + + self.log.debug('Starting to apply radiation damage model...') + #at this point we can give fake data... + cti = CTI.CDM03bidir(self.information, [], log=self.log) + #here we need the right input data + self.image_b = cti.applyRadiationDamage(self.image_b.copy().transpose(), iquadrant=self.information['quadrant']).transpose() + self.log.info('Radiation damage added.') + + + self.log.debug('Starting to apply radiation damage model...') + #at this point we can give fake data... + cti = 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.') +################################################################################## + + + 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_b = IFSinstrumentModel.CCDnonLinearityModel(self.image_b.copy()) + + self.log.info('Non-linearity effects included.') + + + self.log.debug('Starting to apply non-linearity model...') + self.image_r = IFSinstrumentModel.CCDnonLinearityModel(self.image_r.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. + """ + + + self.log.info('readnoise added in blue channel') + ########## blue zone 1 + np.random.seed() + self.image_b[0:1856,0:3392] += np.random.normal(loc=0.0, scale=self.information['rn1_b'], size=(1856,3392)) + + ########## zone 4 ################# + np.random.seed() + self.image_b[1856:3712,0:3392] += np.random.normal(loc=0.0, scale=self.information['rn4_b'], size=(1856,3392)) + + ########## zone 2 ################### + np.random.seed() + self.image_b[0:1856,3392:6784] += np.random.normal(loc=0.0, scale=self.information['rn2_b'], size=(1856,3392)) + + ########## zone 3 + np.random.seed() + self.image_b[1856:3712,3392:6784]+= np.random.normal(loc=0.0, scale=self.information['rn3_b'], size=(1856,3392)) + + + ############################################################################ + + self.log.info('readnoise added in blue channel') + + + ########## red zone 1 + np.random.seed() + self.image_r[0:1856, 0:3392] += np.random.normal(loc=0.0, scale=self.information['rn1_r'], size=(1856,3392)) + + ########## zone 4 ################# + np.random.seed() + self.image_r[1856:3712,0:3392] += np.random.normal(loc=0.0, scale=self.information['rn4_r'], size=(1856,3392)) + + ########## zone 2 ################### + np.random.seed() + self.image_r[0:1856,3392:6784] += np.random.normal(loc=0.0, scale=self.information['rn2_r'], size=(1856,3392)) + + ########## zone 3 + np.random.seed() + self.image_r[1856:3712,3392:6784]+= np.random.normal(loc=0.0, scale=self.information['rn3_r'], size=(1856,3392)) + +########################################################################################## + + def electrons2ADU(self): + """ + Convert from electrons to ADUs using the value read from the configuration file. + """ + + + ############################################################### + + self.log.info('Converting from electrons to ADUs using a factor of gain' ) + + ########## blue zone 1 + self.image_b[0:1856,0:3392] /= self.information['gain1_b'] + + ########## zone 4 ################# + self.image_b[1856:3712,0:3392] /= self.information['gain4_b'] + + ########## zone 2 ################### + self.image_b[0:1856,3392:6784] /= self.information['gain2_b'] + + ########## zone 3 + self.image_b[1856:3712,3392:6784]/= self.information['gain3_b'] + + ############################################################################ + + ########## red zone 1 + self.image_r[0:1856, 0:3392] /= self.information['gain1_r'] + + ########## zone 4 ################# + self.image_r[1856:3712,0:3392] /= self.information['gain4_r'] + + ########## zone 2 ################### + self.image_r[0:1856,3392:6784] /= self.information['gain2_r'] + + ########## zone 3 + self.image_r[1856:3712,3392:6784] /= self.information['gain3_r'] + + ##########################################################################3 + + + #################################################################################### + + 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). + """ + + + ########## blue zone 1 + self.image_b[0:1856,0:3392] += self.information['bias1_b'] + + ########## zone 4 ################# + self.image_b[1856:3712,0:3392] += self.information['bias4_b'] + + ########## zone 2 ################### + self.image_b[0:1856,3392:6784] += self.information['bias2_b'] + + ########## zone 3 + self.image_b[1856:3712,3392:6784] += self.information['bias3_b'] + + self.log.info('Bias counts were added to the blue image' ) + + ############################################################################ + + ########## red zone 1 + self.image_r[0:1856, 0:3392] += self.information['bias1_r'] + + ########## zone 4 ################# + self.image_r[1856:3712,0:3392] += self.information['bias4_r'] + + ########## zone 2 ################### + self.image_r[0:1856,3392:6784] += self.information['bias2_r'] + + ########## zone 3 + self.image_r[1856:3712,3392:6784] += self.information['bias3_r'] + + ########################################################################## + + self.log.info('Bias counts were added to the red image' ) + +############################################################################### + +############################################################################### + + def applyBleeding_yan(self): + """ + Apply bleeding along the CCD columns if the number of electrons in a pixel exceeds the full-well capacity. + + Bleeding is modelled in the parallel direction only, because the CCD273s are assumed not to bleed in + serial direction. + + :return: None + """ + + + if self.image_b.max()>self.information['fullwellcapacity']: + + self.log.info('Applying column bleeding to blue CCD image...') + + #loop over each column, as bleeding is modelled column-wise + for i, column in enumerate(self.image_b.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 + self.image_b[j, i] -= overload + + sum += overload + + elif sum > 0.: + if -overload > sum: + overload = -sum + + self.image_b[j, i] -= overload + sum += overload + + for i, column in enumerate(self.image_b.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 + self.image_b[-j-1, i] -= overload + + sum += overload + elif sum > 0.: + if -overload > sum: + overload = -sum + #self.image[-j-1, i] -= overload + self.image_b[-j-1, i] -= overload + + + sum += overload + print('Applying column bleeding to blue image finished.......') + + + ###################################################################### + if self.image_r.max()>self.information['fullwellcapacity']: + + self.log.info('Applying column bleeding to red CCD image...') + + for i, column in enumerate(self.image_r.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 + self.image_r[j, i] -= overload + + sum += overload + elif sum > 0.: + if -overload > sum: + overload = -sum + #self.image[j, i] -= overload + self.image_r[j, i] -= overload + + + sum += overload + + for i, column in enumerate(self.image_r.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 + self.image_r[-j-1, i] -= overload + + sum += overload + elif sum > 0.: + if -overload > sum: + overload = -sum + #self.image[-j-1, i] -= overload + self.image_r[-j-1, i] -= overload + + + sum += overload + print('Applying column bleeding to red image finished.......') + + ############################################################################ + + ############################################################################ + + 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_b[self.image_b < 0.0] = 0. + #cut of the values larger than max + self.image_b[self.image_b > max] = max + + self.image_b = np.rint(self.image_b).astype(int) + self.log.info('Maximum and total values of the image are %i and %i, respectively' % (np.max(self.image_b), + np.sum(self.image_b))) + + #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))) + + ################################################################################################## + + def applyImageShift(self): + + np.random.seed() + ud= np.random.random() # Choose a random rotation + dx= 2* (ud-0.5) * self.information['shiftmax'] + + np.random.seed() + ud= np.random.random() # Choose a random rotation + dy= 2* (ud-0.5) * self.information['shiftmax'] + + self.image_b= ndimage.shift(self.image_b.copy(), [dy+self.information['shift_b_y'] , dx+self.information['shift_b_x']], order=0, mode='nearest') + self.image_r= ndimage.shift(self.image_r.copy(), [dy+self.information['shift_r_y'] , dx+self.information['shift_r_x']], order=0, mode='nearest') + + self.log.info('Applied image shifting to g r i channels.') + self.information['ra'] = dx*self.information['pixel_size'] + self.information['dec'] = dy*self.information['pixel_size'] + + +######################################################################################################33 + + def applyImageRotate(self ): + np.random.seed() + ud = np.random.random() # Choose a random rotation + angle = 2* (ud-0.5) * self.information['tel_rotmax'] + + inputimg=self.image_b.copy() + rotimg=ndimage.rotate(inputimg,angle+self.information['rotate_b'],order=1,reshape=False) # here we choose reshape=False, the rotated image will + self.image_b=rotimg + + inputimg=self.image_r.copy() + rotimg=ndimage.rotate(inputimg,angle+self.information['rotate_r'],order=1,reshape=False) # here we choose reshape=False, the rotated image will + self.image_r=rotimg + + self.information['Tel_rot']=angle + + self.log.info('Applied telescope rotation with angle (in degree)= %f.', angle) +############################################################################### + + def CCDreadout(self): + + imgb=self.image_b.copy() + temp=np.zeros((3712,6784)) + ########## zone 1 + x1=0 + x2=x1+1024 + + y1=0 + y2=y1+2048 + temp[x1:x2,y1:y2]=imgb[0:1024,0:2048] + ########## zone 4 ################# + x1=2688 + x2=x1+1024 + + y1=0 + y2=y1+2048 + temp[x1:x2,y1:y2]=imgb[1024:2048,0:2048] + ########## zone 2 ################### + x1=0 + x2=x1+1024 + + y1=6784-2048 + y2=y1+2048 + temp[x1:x2,y1:y2]=imgb[0:1024,2048:4096] + ########## zone 3 + x1=2688 + x2=x1+1024 + + y1=6784-2048 + y2=y1+2048 + temp[x1:x2,y1:y2]=imgb[1024:2048,2048:4096] + + self.image_b=temp + + ############################################################################## + imgr=self.image_r.copy() + temp=np.zeros((3712,6784)) + ########## zone 1 + x1=0 + x2=x1+1536 + + y1=0 + y2=y1+3072 + temp[x1:x2,y1:y2]=imgr[0:1536, 0:3072] + ########## zone 4 ################# + x1=2176 + x2=x1+1536 + + y1=0 + y2=y1+3072 + temp[x1:x2,y1:y2]=imgr[1536:3712,0:3072] + ########## zone 2 ################### + x1=0 + x2=x1+1536 + + y1=6784-3072 + y2=y1+3072 + temp[x1:x2,y1:y2]=imgr[0:1536,3072:6144] + ########## zone 3 + x1=2176 + x2=x1+1536 + + y1=6784-3072 + y2=y1+3072 + temp[x1:x2,y1:y2]=imgr[1536:3072,3072:6144] + + self.image_r=temp + + return +############################################################################## + + def writeOutputs(self): + """ + 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 + self.source='sci' + now=datetime.now() + data_time=now.strftime("%Y-%m-%d %H:%M:%S") + exp_endtime=now.strftime("%Y%m%d%H%M%S") + start=now-timedelta(seconds=self.information['exptime']) + exp_starttime=start.strftime("%Y%m%d%H%M%S") + + #write the actual file + obsid=300000000+1 + + data_time=self.dt.strftime("%Y-%m-%d %H:%M:%S") + + 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") + t3=self.dt+timedelta(seconds=self.information['exptime'])+timedelta(seconds=self.information['readouttime']) + + filename_b='CSST_IFS_B_'+self.source+'_'+exp_starttime+'_'+exp_endtime+'_'+str(obsid)+'_X_L0_VER_'+self.information['img_ver']+'.fits' + file_b=self.result_path+'/sky_Data/'+filename_b + + filename_r='CSST_IFS_R_'+self.source+'_'+exp_starttime+'_'+exp_endtime+'_'+str(obsid)+'_X_L0_VER_'+self.information['img_ver']+'.fits' + file_r=self.result_path+'/sky_Data/'+filename_r + + + #create a new FITS file, using HDUList instance + ofd_b = fits.PrimaryHDU() + + ofd_b.header['GROUPS']=( bool(False), 'always F') + ofd_b.header['DATE'] =( data_time, 'date this file was written' ) + + ofd_b.header['FILENAME']=(filename_b, ' file name C48 ') + ofd_b.header['OBSTYPE'] =( self.source, 'observation type raw,flt, mask, bias, dark, sci') + ofd_b.header['TELESCOP']=('CSST', 'always CSST') + ofd_b.header['INSTRUME']=( 'IFS', ' ') + ofd_b.header['RADECSYS']=('ICRS', ' always ICRS ') + ofd_b.header['EQUINOX'] =( float(2000.0), 'always 2000.0') + ofd_b.header['FITSCREA']=( '4.2.1', 'FITS create software version') + + ######### Object information ############# + + ofd_b.header['OBJECT']=( self.information['name_obj'], 'object name') + ofd_b.header['TARGET']=( (self.information['target']), 'target name, hhmmss+ddmmss') + ofd_b.header['OBJ_RA'] =(np.float64(self.information['ra_obj']) , 'RA of the object in deg') + ofd_b.header['OBJ_DEC']=(np.float64(self.information['dec_obj']) , 'DEC of the object in deg') + + ofd_b.header['RA_PNT0']=( np.float64(self.information['ra_pnt0']) , 'RA of the pointing (degrees) at EXPSTART') + ofd_b.header['DEC_PNT0']=(np.float64(self.information['dec_pnt0']) , 'DEC of the pointing (degrees) at EXPSTART') + + + + ############## + ofd_b.header['OBSID'] =(str(obsid) , 'observation ID, 3+8bit') + + ######## Telescope information ############### + # ofd_b.header['COMMENT'] ='==========================================================' + # ofd_b.header['COMMENT'] ='Telescope information' + # ofd_b.header['COMMENT'] ='==========================================================' + + ofd_b.header['REFFRAME']=('CSSTGSC-1.0' , 'guide star catalog version') + ofd_b.header['DATE-OBS']=(data_time , 'date of the observation (yyyy-mm-dd hh:mm:ss)') + + ofd_b.header['EXPSTART']=(np.float64(time2jd(self.dt)), 'exposure start time') + ofd_b.header['SUNANGL0']=(np.float32(0.0) , 'angle between sun and optical axis at EXPSTART') + ofd_b.header['MOONANG0']=(np.float32(0.0) , 'angle between moon and optical axis at EXPSTART') + ofd_b.header['POS_ANG0']=(np.float64(0.0), 'angle between optical axis and the North Pole at EXPSTART in arcsec') + ofd_b.header['TEL_ALT0']=(np.float64(0.0), 'angle between optical axis and the ground- piston at EXPSTART in deg') + + ofd_b.header['HOODSTA0']=(np.float32(0.0) , 'lens hood altitude at EXPSTART') + ofd_b.header['HOODANG0']=(np.float32(0.0), 'lens hood azimuth at EXPSTART') + ofd_b.header['POSI0_X'] =(np.float64(self.information['POSI0_X']) , 'the orbital position of CSST in X direction at EXPSTART') + ofd_b.header['POSI0_Y'] =(np.float64(self.information['POSI0_Y']) , 'the orbital position of CSST in Y direction at EXPSTART') + ofd_b.header['POSI0_Z'] =(np.float64(self.information['POSI0_Z']) , 'the orbital position of CSST in Z direction at EXPSTART') + ofd_b.header['VELO0_X']=( np.float64(self.information['VELO0_X']) , 'the orbital velocity of CSST in X direction at EXPSTART') + ofd_b.header['VELO0_Y']=( np.float64(self.information['VELO0_Y']) , 'the orbital velocity of CSST in Y direction at EXPSTART') + ofd_b.header['VELO0_Z']=( np.float64(self.information['VELO0_Z']) , 'the orbital velocity of CSST in Z direction at EXPSTART') + + ofd_b.header['Euler0_1']=( np.float64(0.0), 'Euler angle 1 at EXPSTART') + ofd_b.header['Euler0_2']=( np.float64(0.0), 'Euler angle 2 at EXPSTART') + ofd_b.header['Euler0_3']=( np.float64(0.0), 'Euler angle 3 at EXPSTART') + + + ofd_b.header['EXPEND'] =(np.float64(time2jd(t2)) , 'exposure end time') + + ofd_b.header['SUNANGL1']=(np.float32(0.0), 'angle between sun and optical axis at EXPEND') + ofd_b.header['MOONANG1']=(np.float32(0.0) , 'angle between moon and optical axis at EXPEND ') + ofd_b.header['POS_ANG1']=(np.float64(0.0) , 'angle between optical axis and the North Pole at EXPEND in arcsec') + ofd_b.header['TEL_ALT1']=(np.float64(0.0) , 'angle between optical axis and the ground- piston at EXPEND in deg ') + ofd_b.header['HOODSTA1']=(np.float32(0.0), 'lens hood altitude at EXPEND ') + ofd_b.header['HOODANG1']=(np.float32(0.0), 'lens hood azimuth at EXPEND ') + + ofd_b.header['POSI1_X'] =(np.float64(self.information['POSI1_X']) , 'the orbital position of CSST in X direction at EXPEND') + ofd_b.header['POSI1_Y'] =(np.float64(self.information['POSI1_Y']) , 'the orbital position of CSST in Y direction at EXPEND') + ofd_b.header['POSI1_Z'] =(np.float64(self.information['POSI1_Z']) , 'the orbital position of CSST in Z direction at EXPEND') + + ofd_b.header['VELO1_X']=(np.float64(self.information['VELO1_X']) , 'the orbital velocity of CSST in X direction at EXPEND') + ofd_b.header['VELO1_Y']=(np.float64(self.information['VELO1_Y']) , 'the orbital velocity of CSST in Y direction at EXPEND') + ofd_b.header['VELO1_Z']=(np.float64(self.information['VELO1_Z']) , 'the orbital velocity of CSST in Z direction at EXPEND') + + + ofd_b.header['Euler1_1']=( np.float64(0.0), 'Euler angle 1 at EXPEND') + ofd_b.header['Euler1_2']=( np.float64(0.0), 'Euler angle 2 at EXPEND') + ofd_b.header['Euler1_3']=( np.float64(0.0), 'Euler angle 3 at EXPEND') + + ofd_b.header['RA_PNT1']=(np.float64(ofd_b.header['RA_PNT0']), 'RA of the pointing (degrees) at EXPEND in deg') + ofd_b.header['DEC_PNT1']=(np.float64(ofd_b.header['DEC_PNT0']), 'DEC of the pointing (degrees) at EXPEND in deg') + + ofd_b.header['EXPTIME']=(self.information['exptime'], 'exposure duration') + ofd_b.header['EPOCH'] =(np.float32(0.0), 'coordinate epoch') + ofd_b.header['CHECKSUM']=( 0 , 'hdu-checksum') + + ########## finish header for 0 layer + #############################################################################3 + ##### header + + b1= self.image_b[1856:3712,0:3392] #b4 + b2= self.image_b[1856:3712,3392:6784] #b3 + b3= self.image_b[0:1856,0:3392] #b1 + b4= self.image_b[0:1856,3392:6784] #b2 + + ####### do Flip the b2 b2 and b4 array in the up/down or left/right direction. + b2=np.fliplr(b2) ## left to right + b3=np.flipud(b3) ## down to up + + b4=np.fliplr(b4) ## left to right and down to up + b4=np.flipud(b4) + + + bb=np.hstack((b1,b2,b3,b4)) + + + #new image HDU, blue channel, layer 1 + hdu_b =fits.ImageHDU(data=np.uint16(bb)) + + ######### instrument information ###### + ##### + hdu_b.header['PMIRRPOS']=(bool(False), 'FSM pointing,T: to MCI, F: not to MCI') + + if self.source =='sci': + hdu_b.header['CMIRRPOS']=(bool(False), 'position of calibration switch mirror,T: for calibration, F: not') + else: + hdu_b.header['CMIRRPOS']=(bool(True), 'position of calibration switch mirror,T: for calibration, F: not') + + if self.source=='flat': + hdu_b.header['FLAMP'] =(int(1), 'status of flat lamp,0: off, 1: , 2: ') + else: + hdu_b.header['FLAMP'] =(int(0), 'status of flat lamp,0: off, 1: , 2: ') + + + if self.source=='lamp': + hdu_b.header['ALAMP'] =(int(1),'status of atomic emission line lamp,0: off, 1: , 2: ') + else: + hdu_b.header['ALAMP'] =(int(0),'status of atomic emission line lamp,0: off, 1: , 2: ') + + ############# + hdu_b.header['IFSMODE'] =(int(0), 'IFS working mode') + hdu_b.header['IFSTEMP'] =(float(0.0), 'IFS components temperature in degC') + hdu_b.header['IFSSTAT'] =(int(0), 'IFS components status parameter') + ############################################################################## + ################### detector information############################# + # hdu_b.header['COMMENT'] ='==========================================================' + # hdu_b.header['COMMENT'] ='Detector information' + # hdu_b.header['COMMENT'] ='==========================================================' + + hdu_b.header['CAMERA'] =('Blue','camera of IFS') + hdu_b.header['DETNAM'] =('CCD231-c4','detector name') + hdu_b.header['DETSIZE'] =('', 'detector size') + hdu_b.header['DATASEC'] =('', 'data section') + hdu_b.header['PIXSCAL1']=(1856, 'pixel scale for axis 1') + hdu_b.header['PIXSCAL2']=(3392*4, 'pixel scale for axis 2') + hdu_b.header['PIXSIZE1']=(15, 'pixel size in um') + hdu_b.header['PIXSIZE2']=(15, 'pixel size in um') + hdu_b.header['NCHAN'] =(4, 'number of readout channels') + hdu_b.header['NCHAN1'] =(2, 'number of horizontal channels') + hdu_b.header['NCHAN2'] =(2, 'number of verticalchannels') + hdu_b.header['PSCAN1'] =(0, 'horizontal prescan width, per readout channel') + hdu_b.header['PSCAN2'] =(0, 'vertical prescan width, per readout channel') + hdu_b.header['OSCAN1'] =(0,' horizontal overscan width, per readout channel') + hdu_b.header['OSCAN2'] =(0, 'vertical overscan width, per readout channel') + + ## Readout information + # hdu_b.header['COMMENT'] ='=============================================================' + # hdu_b.header['COMMENT'] ='Readout information' + # hdu_b.header['COMMENT'] ='=============================================================' + + + hdu_b.header['READT0'] =(np.float64(time2jd(t2)),'read start time (UTC)') + hdu_b.header['READT1'] =(np.float64(time2jd(t3)), 'read end time (UTC)') + + hdu_b.header['DETTEMP0']=(np.float32(0.0), 'detector temperature at READT0') + hdu_b.header['DETTEMP1']=(np.float32(0.0), 'detector temperature at READT1') + hdu_b.header['BIN_X'] =(0, 'bin number in X (wavelength)') + hdu_b.header['BIN_Y'] =(0, 'bin number in Y (spatial)') + + hdu_b.header['GAIN1'] =(self.information['gain4_b'],'CCD gain (channel 1)') + hdu_b.header['GAIN2'] =(self.information['gain3_b'],'CCD gain (channel 2)') + hdu_b.header['GAIN3'] =(self.information['gain1_b'],'CCD gain (channel 3)') + hdu_b.header['GAIN4'] =(self.information['gain2_b'],'CCD gain (channel 4)') + + hdu_b.header['DARK1'] =(self.information['dark4_b'],'CCD dark (channel 1)') + hdu_b.header['DARK2'] =(self.information['dark3_b'],'CCD dark (channel 2)') + hdu_b.header['DARK3'] =(self.information['dark1_b'],'CCD dark (channel 3)') + hdu_b.header['DARK4'] =(self.information['dark2_b'],'CCD dark (channel 4)') + + + hdu_b.header['RDNOIS1'] =(self.information['rn4_b'],'read noise (channel 1') + hdu_b.header['RDNOIS2'] =(self.information['rn3_b'],'read noise (channel 2') + hdu_b.header['RDNOIS3'] =(self.information['rn1_b'],'read noise (channel 3') + hdu_b.header['RDNOIS4'] =(self.information['rn2_b'],'read noise (channel 4') + + hdu_b.header['DETBIA1'] =(self.information['bias4_b'],'amplifier bias voltage (channel1)') + hdu_b.header['DETBIA2'] =(self.information['bias3_b'],'amplifier bias voltage (channel2)') + hdu_b.header['DETBIA3'] =(self.information['bias1_b'],'amplifier bias voltage (channel3)') + hdu_b.header['DETBIA4'] =(self.information['bias2_b'],'amplifier bias voltage (channel4)') + + hdu_b.header['RDSPEED'] =(100,'read speed (in MHz)') + + hdu_b.header['EXPTIME'] =(self.information['exptime'],'exposure time in seconds') + + hdu_b.header['Img_Ver'] =(self.information['img_ver'], 'IFS CCD image Version') + + hdu_b.header['sky_obj'] =(self.skyfilepath,'input sky fits filepath') + + + ########################################################## + #################### red camera ###################### + + + #create a new FITS file, using HDUList instance + ofd_r = fits.PrimaryHDU() + + ofd_r.header['GROUPS']=( bool(False), 'always F') + ofd_r.header['DATE'] =( data_time, 'date this file was written' ) + + ofd_r.header['FILENAME']=(filename_r, ' file name C48 ') + + ofd_r.header['OBSTYPE'] =( self.source, 'observation type raw,flt, mask, bias, dark, sci') + + ofd_r.header['TELESCOP']=('CSST' , 'always CSST') + ofd_r.header['INSTRUME']=( 'IFS', ' ') + ofd_r.header['RADECSYS']=('ICRS', ' always ICRS ') + ofd_r.header['EQUINOX'] =( float(2000.0), 'always 2000.0') + ofd_r.header['FITSCREA']=( '4.2.1' , 'FITS create software version') + ######### Object information ############# + + ######### Object information ############# + # ofd_r.header['COMMENT']='=======================================================================' + # ofd_r.header['COMMENT']='Object information' + # ofd_r.header['COMMENT']='=======================================================================' + + ofd_r.header['OBJECT']=( self.information['name_obj'], 'object name') + ofd_r.header['TARGET']=( (self.information['target']), 'target name, hhmmss+ddmmss') + ofd_r.header['OBJ_RA'] =(np.float64(self.information['ra_obj']) , 'RA of the object in deg') + ofd_r.header['OBJ_DEC']=(np.float64(self.information['dec_obj']) , 'DEC of the object in deg') + + ofd_r.header['RA_PNT0']=( np.float64(self.information['ra_pnt0']) , 'RA of the pointing (degrees) at EXPSTART') + ofd_r.header['DEC_PNT0']=(np.float64(self.information['dec_pnt0']) , 'DEC of the pointing (degrees) at EXPSTART') + + ofd_r.header['OBSID'] =(str(obsid) , 'observation ID, 3+8bit ') + + ######## Telescope information ############### + # ofd_r.header['COMMENT']='=======================================================================' + # ofd_r.header['COMMENT']='Telescope information' + # ofd_r.header['COMMENT']='=======================================================================' + + ofd_r.header['REFFRAME']=('CSSTGSC-1.0' , 'guide star catalog version ') + ofd_r.header['DATE-OBS']=(data_time , 'date of the observation (yyyy-mm-dd hh:mm:ss)') + + ofd_r.header['EXPSTART']=(np.float64(exp_starttime), 'exposure start time ') + ofd_r.header['SUNANGL0']=(np.float32(0.0) , 'angle between sun and optical axis at EXPSTART ') + ofd_r.header['MOONANG0']=(np.float32(0.0) , 'angle between moon and optical axis at EXPSTART ') + ofd_r.header['POS_ANG0']=(np.float64(0.0), 'angle between optical axis and the North Pole at EXPSTART in arcsec ') + ofd_r.header['TEL_ALT0']=(np.float64(0.0), 'angle between optical axis and the ground- piston at EXPSTART in deg ') + ofd_r.header['HOODSTA0']=(np.float32(0.0) , 'lens hood altitude at EXPSTART ') + ofd_r.header['HOODANG0']=(np.float32(0.0), 'lens hood azimuth at EXPSTART ') + ofd_r.header['POSI0_X'] =(np.float64(self.information['POSI0_X']) , 'the orbital position of CSST in X direction at EXPSTART') + ofd_r.header['POSI0_Y'] =(np.float64(self.information['POSI0_Y']) , 'the orbital position of CSST in Y direction at EXPSTART') + ofd_r.header['POSI0_Z'] =(np.float64(self.information['POSI0_Z']) , 'the orbital position of CSST in Z direction at EXPSTART') + ofd_r.header['VELO0_X']=( np.float64(self.information['VELO0_X']) , 'the orbital velocity of CSST in X direction at EXPSTART') + ofd_r.header['VELO0_Y']=( np.float64(self.information['VELO0_Y']) , 'the orbital velocity of CSST in Y direction at EXPSTART') + ofd_r.header['VELO0_Z']=( np.float64(self.information['VELO0_Z']) , 'the orbital velocity of CSST in Z direction at EXPSTART') + + ofd_r.header['Euler0_1']=( np.float64(0.0), 'Euler angle 1 at EXPSTART') + ofd_r.header['Euler0_2']=( np.float64(0.0), 'Euler angle 2 at EXPSTART') + ofd_r.header['Euler0_3']=( np.float64(0.0), 'Euler angle 3 at EXPSTART') + + ofd_r.header['EXPEND'] =(np.float64(exp_endtime) , 'exposure end time') + + ofd_r.header['SUNANGL1']=(np.float32(0.0), 'angle between sun and optical axis at EXPEND') + ofd_r.header['MOONANG1']=(np.float32(0.0) , 'angle between moon and optical axis at EXPEND ') + ofd_r.header['POS_ANG1']=(np.float64(0.0) , 'angle between optical axis and the North Pole at EXPEND in arcsec') + ofd_r.header['TEL_ALT1']=(np.float64(0.0) , 'angle between optical axis and the ground- piston at EXPEND in deg ') + ofd_r.header['HOODSTA1']=(np.float32(0.0), 'lens hood altitude at EXPEND ') + ofd_r.header['HOODANG1']=(np.float32(0.0), 'lens hood azimuth at EXPEND ') + + ofd_r.header['POSI1_X'] =(np.float64(self.information['POSI1_X']) , 'the orbital position of CSST in X direction at EXPEND') + ofd_r.header['POSI1_Y'] =(np.float64(self.information['POSI1_Y']) , 'the orbital position of CSST in Y direction at EXPEND') + ofd_r.header['POSI1_Z'] =(np.float64(self.information['POSI1_Z']) , 'the orbital position of CSST in Z direction at EXPEND') + + ofd_r.header['VELO1_X']=(np.float64(self.information['VELO1_X']) , 'the orbital velocity of CSST in X direction at EXPEND') + ofd_r.header['VELO1_Y']=(np.float64(self.information['VELO1_Y']) , 'the orbital velocity of CSST in Y direction at EXPEND') + ofd_r.header['VELO1_Z']=(np.float64(self.information['VELO1_Z']) , 'the orbital velocity of CSST in Z direction at EXPEND') + + ofd_r.header['Euler1_1']=( np.float64(0.0), 'Euler angle 1 at EXPEND') + ofd_r.header['Euler1_2']=( np.float64(0.0), 'Euler angle 2 at EXPEND') + ofd_r.header['Euler1_3']=( np.float64(0.0), 'Euler angle 3 at EXPEND') + + ofd_r.header['RA_PNT1']=(np.float64(ofd_r.header['RA_PNT0']), 'RA of the pointing (degrees) at EXPEND in deg') + ofd_r.header['DEC_PNT1']=(np.float64(ofd_r.header['DEC_PNT0']), 'DEC of the pointing (degrees) at EXPEND in deg') + + ofd_r.header['EXPTIME']=(self.information['exptime'], 'exposure duration') + ofd_r.header['EPOCH'] =(np.float32(0.0), 'coordinate epoch') + ofd_r.header['CHECKSUM']=( 0 , 'hdu-checksum') + + ### finish 0 layer header + + ########## finish header for 0 layer + + # ########## blue zone 1--to--3 + # self.image_r[0:1856,0:3392] += self.information['bias1_r'] + + # ########## zone 4 --to---1 ################# + # self.image_r[1856:3712,0:3392] += self.information['bias4_r'] + + # ########## zone 2 ----to----4 ################### + # self.image_r[0:1856,3392:6784] += self.information['bias2_r'] + + # ########## zone 3 ---to------2 + # self.image_r[1856:3712,3392:6784] += self.information['bias3_r'] + #############################################################################3 + ##### header + + b1= self.image_r[1856:3712,0:3392] + b2= self.image_r[1856:3712,3392:6784] + b3= self.image_r[0:1856,0:3392] + b4= self.image_r[0:1856,3392:6784] + + + ####### do Flip the b2 b2 and b4 array in the up/down or left/right direction. + b2=np.fliplr(b2) ## left to right + b3=np.flipud(b3) ## down to up + + b4=np.fliplr(b4) ## left to right and down to up + b4=np.flipud(b4) + + + + rr=np.hstack((b1,b2,b3,b4)) + #new image HDU, blue channel, layer 1 + hdu_r =fits.ImageHDU(data=np.uint16(rr)) + + + ######################################### + ######### instrument information ###### + hdu_r.header['PMIRRPOS']=(bool(False), 'FSM pointing,T: to MCI, F: not to MCI') + + + hdu_r.header['CMIRRPOS']=(bool(False), 'position of calibration switch mirror,T: for calibration, F: not') + + hdu_r.header['FLAMP'] =(int(0), 'status of flat lamp,0: off, 1: , 2: ') + + hdu_r.header['ALAMP'] =(int(0),'status of atomic emission line lamp,0: off, 1: , 2: ') + + hdu_r.header['IFSMODE'] =(int(0), 'IFS working mode') + hdu_r.header['IFSTEMP'] =(float(0.0), 'IFS components temperature in degC') + hdu_r.header['IFSSTAT'] =(int(0), 'IFS components status parameter') + + ################### detector information############################# + # hdu_r.header['COMMENT']='=======================================================================' + # hdu_r.header['COMMENT']='Detector information' + # hdu_r.header['COMMENT']='=======================================================================' + + hdu_r.header['CAMERA'] =('Red','camera of IFS') + hdu_r.header['DETNAM'] =('CCD231-c4','detector name') + hdu_r.header['DETSIZE'] =('', 'detector size') + hdu_r.header['DATASEC'] =('', 'data section') + hdu_r.header['PIXSCAL1']=(1856, 'pixel scale for axis 1') + hdu_r.header['PIXSCAL2']=(3392*4, 'pixel scale for axis 2') + hdu_r.header['PIXSIZE1']=(15, 'pixel size in um') + hdu_r.header['PIXSIZE2']=(15, 'pixel size in um') + hdu_r.header['NCHAN'] =(4, 'number of readout channels') + hdu_r.header['NCHAN1'] =(2, 'number of horizontal channels') + hdu_r.header['NCHAN2'] =(2, 'number of verticalchannels') + hdu_r.header['PSCAN1'] =(0, 'horizontal prescan width, per readout channel') + hdu_r.header['PSCAN2'] =(0, 'vertical prescan width, per readout channel') + hdu_r.header['OSCAN1'] =(0,' horizontal overscan width, per readout channel') + hdu_r.header['OSCAN2'] =(0, 'vertical overscan width, per readout channel') + +##################################################################################################### + ## Readout information + # hdu_r.header['COMMENT']='=======================================================================' + # hdu_r.header['COMMENT']='Readout information' + # hdu_r.header['COMMENT']='=======================================================================' + + + + hdu_r.header['READT0'] =(np.float64(time2jd(t2)),'read start time (UTC)') + hdu_r.header['READT1'] =(np.float64(time2jd(t3)), 'read end time (UTC)') + + hdu_r.header['DETTEMP0']=(np.float32(0.0), 'detector temperature at READT0') + hdu_r.header['DETTEMP1']=(np.float32(0.0), 'detector temperature at READT1') + hdu_r.header['BIN_X'] =(0, 'bin number in X (wavelength)') + hdu_r.header['BIN_Y'] =(0, 'bin number in Y (spatial)') + + hdu_r.header['GAIN1'] =(self.information['gain4_r'],'CCD gain (channel 1)') + hdu_r.header['GAIN2'] =(self.information['gain3_r'],'CCD gain (channel 2)') + hdu_r.header['GAIN3'] =(self.information['gain1_r'],'CCD gain (channel 3)') + hdu_r.header['GAIN4'] =(self.information['gain2_r'],'CCD gain (channel 4)') + + + hdu_r.header['DARK1'] =(self.information['dark4_r'],'CCD dark (channel 1)') + hdu_r.header['DARK2'] =(self.information['dark3_r'],'CCD dark (channel 2)') + hdu_r.header['DARK3'] =(self.information['dark1_r'],'CCD dark (channel 3)') + hdu_r.header['DARK4'] =(self.information['dark2_r'],'CCD dark (channel 4)') + + + hdu_r.header['RDNOIS1'] =(self.information['rn4_r'],'read noise (channel 1') + hdu_r.header['RDNOIS2'] =(self.information['rn3_r'],'read noise (channel 2') + hdu_r.header['RDNOIS3'] =(self.information['rn1_r'],'read noise (channel 3') + hdu_r.header['RDNOIS4'] =(self.information['rn2_r'],'read noise (channel 4') + + hdu_r.header['DETBIA1'] =(self.information['bias4_r'],'amplifier bias voltage (channel1)') + hdu_r.header['DETBIA2'] =(self.information['bias3_r'],'amplifier bias voltage (channel2)') + hdu_r.header['DETBIA3'] =(self.information['bias1_r'],'amplifier bias voltage (channel3)') + hdu_r.header['DETBIA4'] =(self.information['bias2_r'],'amplifier bias voltage (channel4)') + + hdu_r.header['RDSPEED'] =(100,'read speed (in MHz)') + + hdu_r.header['EXPTIME'] =(self.information['exptime'],'exposure time in seconds') + + hdu_r.header['Img_Ver'] =(self.information['img_ver'], 'IFS CCD image Version') + + hdu_r.header['sky_obj'] =(self.skyfilepath,'input sky fits filename') + + + + hdulist_b=fits.HDUList([ofd_b, hdu_b]) + hdulist_b.writeto(file_b, overwrite=True) + #print('IFS_b.fits is created ') + + hdulist_r=fits.HDUList([ofd_r, hdu_r]) + hdulist_r.writeto(file_r, overwrite=True) + #print('IFS_r.fits is created ') + ################################################################################## + + def earthshine(self,theta): + """ + For given theta angle, return the earth-shine spectrum. + + :param theta: angle (in degree) from the target to earth limb. + :return: the scaled solar spectrum + template_wave: unit in A + template_flux: unit in erg/s/cm^2/A/arcsec^2 + + """ + + # read solar template + solar_template = pd.read_csv(self.information['indata_path']+'/refs/solar_spec.dat', sep='\s+', + header=None, comment='#') + template_wave = solar_template[0].values + template_flux = solar_template[1].values + + # read earth shine surface brightness + earthshine_curve = pd.read_csv(self.information['indata_path']+'/refs/earthshine.dat', + header=None, comment='#') + angle = earthshine_curve[0].values + surface_brightness = earthshine_curve[1].values + + # read V-band throughtput + cat_filter_V = pd.read_csv(self.information['indata_path']+'/refs/filter_Bessell_V.dat', sep='\s+', + header=None, comment='#') + filter_wave = cat_filter_V[0].values + filter_response = cat_filter_V[1].values + + # interplate to the target wavelength in V-band + ind_filter = (template_wave >= np.min(filter_wave)) & (template_wave <= np.max(filter_wave)) + filter_wave_interp = template_wave[ind_filter] + filter_response_interp = np.interp(filter_wave_interp, filter_wave, filter_response) + + filter_constant = simps(filter_response_interp * filter_wave_interp, filter_wave_interp) + template_constant = simps(filter_response_interp * template_wave[ind_filter] * template_flux[ind_filter], + template_wave[ind_filter]) + dwave = filter_wave_interp[1:] - filter_wave_interp[:-1] + wave_eff = np.nansum(dwave * filter_wave_interp[1:] * filter_response_interp[1:]) / \ + np.nansum(dwave * filter_response_interp[1:]) + + # get the normalized value at theta. + u0 = np.interp(theta, angle, surface_brightness) # mag/arcsec^2 + u0 = 10**((u0 + 48.6)/(-2.5)) # target flux in erg/s/cm^2/Hz unit + u0 = u0 * 3e18 / wave_eff**2 # erg/s/cm^2/A/arcsec^2 + + factor = u0 * filter_constant / template_constant + norm_flux = template_flux * factor # erg/s/cm^2/A/arcsec^2 + + self.earthshine_wave=template_wave # A + self.earthshine_flux=norm_flux + + return + + +######################################################################################################################################################################################################################################################## + + + 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') + 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['indata_path']+'/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^−8 W m−2 sr−1 um−1 + + # read the zodical spectrum in the ecliptic + cat_spec = pd.read_csv(self.information['indata_path']+'/refs/solar_spec.dat', sep='\s+', header=None, comment='#') + wave = cat_spec[0].values # A + spec0 = cat_spec[1].values # 10^-8 W m^−2 sr^−1 μm^−1 + zodi_norm = 252 # 10^-8 W m^−2 sr^−1 μm^−1 + + spec = spec0 * (zodi_obj / zodi_norm) * 1e-8 # W m^−2 sr^−1 μm^−1 + + # 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 + ################################################################################### + ################################################################################## + def CalskyNoise(self, lam): + + + # calculate sky noise; + planckh= 6.62620*10**-27 # erg s; + cc=2.99792458*10**17 # in nm/s + fov2=0.01 # arcsec^2 + + # lam is input wavelength in nm + + ########################################## + self.earthshine_wave # A + self.earthshine_flux # erg/s/cm^2/A/arcsec^2 + earthshine_flux=np.interp(lam*10.0, self.earthshine_wave,self.earthshine_flux) # flux from zodiacal + + self.zodiacal_wave # in A + self.zodiacal_flux # erg/s/cm^2/A/arcsec^2 + zodiacal_flux=np.interp(lam*10.0, self.zodiacal_wave,self.zodiacal_flux) # flux from zodiacal + fluxlam_sky=(earthshine_flux+zodiacal_flux)*fov2 # erg/s/cm2/A + ############### + + ephoton=planckh*cc/lam # in erg/photon, cc与lambda单位需要一致； + Ns_skynoise=fluxlam_sky/ephoton # in unit of photons/cm2/s/A + + + return Ns_skynoise + + + ################################################################################# + def sim_sky_img(self,skyfitsfilename, skyRa_shift, skyDec_shift, sky_rot, exposuretime): + + + ############################################################################ + ### load fits file + indatafile=skyfitsfilename + a=fits.open(indatafile) + ####################################### + self.information['name_obj']=a[0].header['OBJECT'] + self.information['ra_obj'] =a[0].header['RA'] ### in degree + self.information['dec_obj'] =a[0].header['DEC'] ### in degree + + disRa =(skyRa_shift)/3600.0 ##convert unit of degree to arcsec + disDec=(skyDec_shift)/3600.0 ##convert unit of degree to arcsec + + self.information['ra_pnt0'] =a[0].header['RA'] + disRa/np.cos(a[0].header['DEC']/180.0*np.pi) + self.information['dec_pnt0']=a[0].header['DEC']+ disDec + + self.earthshine(self.earthshine_theta) + + self.zodiacal(self.information['ra_obj'], self.information['dec_obj'], self.zodiacal_time) + + self.information['target']=deg2HMS(self.information['ra_obj'], self.information['dec_obj']) + + + ### main input data + SpecCube=a[1].data ## spectrum data cube; + + Wave=0.1*a[2].data # the relatived wavelength which is converted from Unit A to nm + #print('Wave data header', hdr) + + ###################################################################################### + + exptime=self.information['exptime'] #exposure time + + dis_dx=disRa # image shift Ra in arcsec + dis_dy=disDec # image shift Dec in arcsec + + sizeout= len(SpecCube[:,0,0]) + + blue_img=galsim.Image(np.zeros((2048,4096)),copy=True) + blue_img.scale=self.pixelscale + blue_img.setOrigin=(0,0) + + red_img=galsim.Image(np.zeros((3072,6144)),copy=True) + red_img.scale=self.pixelscale + red_img.setOrigin=(0,0) + + blue_sensor=galsim.Sensor() + + red_sensor=galsim.Sensor() + + deltalam=np.mean(np.diff(Wave)) + + energy=0.0 + + energy_blue=0.0 + energy_red=0.0 + + width_blue=0 + ################################ + ############## doppler effect to photons.wavelength ############# + + #self.orbit_pars + 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]) + vx_sat=float(self.orbit_pars[self.orbit_exp_num,4]) + vy_sat=float(self.orbit_pars[self.orbit_exp_num,5]) + vz_sat=float(self.orbit_pars[self.orbit_exp_num,6]) + + self.information['POSI0_X']=x_sat + self.information['POSI0_Y']=y_sat + self.information['POSI0_Z']=z_sat + self.information['VELO0_X']=vx_sat + self.information['VELO0_Y']=vy_sat + self.information['VELO0_Z']=vz_sat + + + theta1=beta_angle( x_sat,y_sat,z_sat,vx_sat,vy_sat,vz_sat, self.information['ra_obj'], self.information['dec_obj']) + + v1=np.sqrt(vx_sat**2+vy_sat**2+vz_sat**2)*np.cos(theta1/180.0*np.pi) # velocity at stat exposure time + vv1=LSR_velocity(self.information['ra_obj'],self.information['dec_obj'],v1, self.TianCe_day) + + + ################################################# + + ### exposure end time is t2 ; + t2=self.dt+timedelta(seconds=self.information['exptime']) + + t2jd=time2jd(t2) + + if self.orbit_pars[-1,0]1e-3) + energy=energy+sum(photons.flux[idx0]) ### totla energy for slice image + p_num=len(idx0[0]) + + ############################################################################################### + ############### find photons for blue channel, and make the flux multiple the optical and CCD efficiency + + np.random.seed() + wavesample=lam1+(lam2-lam1)*np.random.rand(p_num) + + + if (lam>=350.0 and lam<=650.0): + ## bulue channel + photons_blue=galsim.PhotonArray(N=p_num, x=photons.x[idx0], y=photons.y[idx0], flux = Qe_blue*photons.flux[idx0],wavelength=wavesample) + + dx_blue, dy_blue = get_dx_dy_blue(wavesample) + + photons_blue.x=photons_blue.x/self.pixelscale+dx_blue+self.slice_blue['px'][k] + + photons_blue.y=photons_blue.y/self.pixelscale+dy_blue+self.slice_blue['py'][k] + + blue_sensor.accumulate(photons_blue, blue_img) + + energy_blue=energy_blue+sum(photons_blue.flux) + + width_blue=width_blue+deltalam/32.0 + + + + if (lam>=560.0) & (lam<=1000.0 ): + ## red channel + photons_red=galsim.PhotonArray(N=p_num,x=photons.x[idx0],y=photons.y[idx0], flux = Qe_red*photons.flux[idx0],wavelength=wavesample) + + dx_red, dy_red = get_dx_dy_red(wavesample) + + photons_red.x=photons_red.x/self.pixelscale+dx_red+self.slice_red['px'][k] + + photons_red.y=photons_red.y/self.pixelscale+dy_red+self.slice_red['py'][k] + + red_sensor.accumulate(photons_red, red_img) + + energy_red=energy_red+sum(photons_red.flux) + + #################################################################################### + ## stray light will cover 2% of input total light; + blue_img.array[:,:]=blue_img.array[:,:]+0.01*energy/2048/4096 + + red_img.array[:,:] =red_img.array[:,:]+ 0.01*energy/3072/6144 + + self.image_b=blue_img.array + self.image_r=red_img.array + + return + +################################################################################################# + +################################################################################################# + + def simulate(self, skyfitsin, skyRa_shift, skyDec_shift,sky_rot, exptime): + + + """ + Create a single simulated image of a quadrant defined by the configuration file. + + + """ + + #self.configure() #print the configfile name and path; + + self.dt=datetime.now() + + self.information['exptime']=exptime + + self.skyfilepath=skyfitsin + + np.random.seed() + ud = np.random.random() # Choose a random + self.earthshine_theta=ud * 60 # in degree + + ################################################################## + #### load orbit parameters ##### + flag=0 + for k in range(1,50,1): + + fn=self.information['indata_path']+'/refs/orbit20160925/'+str(k)+'.txt'; + d=np.loadtxt(fn); + self.dt_num=int((self.information['exptime']+self.information['readouttime']+125)/120) + now_dt=datetime.utcnow() + now_jd=time2jd(now_dt) + for kk in range(len(d[:,0])): + if now_jd-d[kk,0]<=0: + flag=1 + break + if flag==1: + break + #####################end for + self.orbit_pars=d + self.orbit_file_num=k + self.orbit_exp_num =kk + + exptime_start_jd=d[kk,0] #### jd time, utc format + + self.dt=julian.from_jd(exptime_start_jd, fmt='jd') + + self.TianCe_day=self.dt.strftime("%Y-%m-%d") ###str(self.dt.year)+'-'+str(self.dt.month)+'-'+str(self.dt.day) + + self.TianCe_exp_start=dt2hmd(self.dt) + + self.zodiacal_time=self.TianCe_day + + ####################################################################### + + self.sim_sky_img(skyfitsin, skyRa_shift, skyDec_shift, sky_rot, self.information['exptime']) + + self.information['sky_rot']=sky_rot + self.information['skyRa_shift']=skyRa_shift + self.information['skyRa_shift']=skyDec_shift + +############################################################################### +############ add some effect to images ####################################### + + + if self.flatfieldM: + self.applyflatfield() + #print('Applying flatfieldM finished.......') + + + if self.darknoise: + self.applyDarkCurrent() + + + if self.cosmicRays: + self.addCosmicRays() + #print('Applying cosmicRays finished.......') + + if self.bleeding: + self.applyBleeding_yan() + #print('Applying bleeding finished.......') + + self.applyPoissonNoise() + + + if self.nonlinearity: + self.applyNonlinearity() + #print('Applying nonlinearity finished.......') + + + if self.radiationDamage: + self.applyRadiationDamage() + #print('Applying radiationDamage finished.......') + + + + ##### cut original CCD image to four parts by four read out channels and zones + self.CCDreadout() + + if self.readoutNoise: + self.applyReadoutNoise() + #print('Applying readoutNoise finished.......') + + self.electrons2ADU() + + self.applyBias() + + if self.cosmetics: + self.applyCosmetics() + #print('Applying cosmetics finished.......') + + self.discretise() + + self.writeOutputs() + + 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 simulation.') + +############################################################################################## +############################################################################################## +def runIFSsim(configfile): + + + opts, args = processArgs() + + opts.configfile=configfile + + simulate= IFSsimulator(opts) + + skyfitsin=simulate.information['skyfitsin'] + + exptime=simulate.information['exptime'] + + sky_ra_dis=simulate.information['sky_ra_dis'] + + sky_dec_dis=simulate.information['sky_dec_dis'] + + sky_angle_dis=simulate.information['sky_angle_dis'] + + + + simulate.simulate(skyfitsin, sky_ra_dis, sky_dec_dis, sky_angle_dis, exptime) + + +########################## begin main fucntion ####################################### +############################################################################################## +############################################################################################## + +if __name__ == "__main__": + + + + if len(sys.argv[:]) <2: + + configfile='./ifs_data/ifs_sim_example.config' + + ########################################################################################### + + if len(sys.argv[:]) >=2: + configfile=sys.argv[1] + if not os.path.exists(configfile): + print('The given input configfile path is wrong......') + sys.exit(1) + + ################################################ + + runIFSsim(configfile) + + print('---The CSST-IFS simulation is successful!---') + + #'/home/yan/IFS_FabuCode/InputData/IFS_sim_Fabu.config' + ############################################################ + + + diff --git a/csst_ifs_sim/ifs_so/__init__.py b/csst_ifs_sim/ifs_so/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/csst_ifs_sim/ifs_so/cdm03.cpython-37m-x86_64-linux-gnu.so b/csst_ifs_sim/ifs_so/cdm03.cpython-37m-x86_64-linux-gnu.so new file mode 100644 index 0000000000000000000000000000000000000000..85a4cfc31920218e8b634b8d53cf9db62547f63f Binary files /dev/null and b/csst_ifs_sim/ifs_so/cdm03.cpython-37m-x86_64-linux-gnu.so differ diff --git a/csst_ifs_sim/ifs_so/cdm03.cpython-38-x86_64-linux-gnu.so b/csst_ifs_sim/ifs_so/cdm03.cpython-38-x86_64-linux-gnu.so new file mode 100644 index 0000000000000000000000000000000000000000..7bfa02f3f33ab5d5a16e4d35ac622d8bf7a49ed1 Binary files /dev/null and b/csst_ifs_sim/ifs_so/cdm03.cpython-38-x86_64-linux-gnu.so differ diff --git a/csst_ifs_sim/ifs_so/cdm03bidir.cpython-37m-x86_64-linux-gnu.so b/csst_ifs_sim/ifs_so/cdm03bidir.cpython-37m-x86_64-linux-gnu.so new file mode 100644 index 0000000000000000000000000000000000000000..a63494c33d3e3f2c0dd42085e549d38415c682c6 Binary files /dev/null and b/csst_ifs_sim/ifs_so/cdm03bidir.cpython-37m-x86_64-linux-gnu.so differ diff --git a/csst_ifs_sim/ifs_so/cdm03bidir.cpython-38-x86_64-linux-gnu.so b/csst_ifs_sim/ifs_so/cdm03bidir.cpython-38-x86_64-linux-gnu.so new file mode 100644 index 0000000000000000000000000000000000000000..4ddf6dd453c91a1ed8b3448956bc127596d748f8 Binary files /dev/null and b/csst_ifs_sim/ifs_so/cdm03bidir.cpython-38-x86_64-linux-gnu.so differ diff --git a/csst_ifs_sim/support/IFSinstrumentModel.py b/csst_ifs_sim/support/IFSinstrumentModel.py new file mode 100644 index 0000000000000000000000000000000000000000..52550f9722521599ed5fa9948184efa96b98936b --- /dev/null +++ b/csst_ifs_sim/support/IFSinstrumentModel.py @@ -0,0 +1,55 @@ +""" +VIS Instrument Model +==================== + +The file provides a function that returns VIS related information such as pixel +size, dark current, gain, zeropoint, and sky background. + +:requires: NumPy +:requires: numexpr + +:version: 0.7 +""" + +# import datetime, math +# import numpy as np +# import numexpr as ne + + +def IFSinformation(): + """ + Returns a dictionary describing VIS. The following information is provided (id: value - reference):: + + + + :return: instrument model parameters + :rtype: dict + """ + + ######################################################################################################### + #out = dict(readnoise=4, pixel_size=0.1, dark=0.0008333, fullwellcapacity=90000, bluesize=4000, redsize=6000, readtime=300.) + out=dict() + out.update({'dob' : 0, 'rdose' : 8.0e9, + 'parallelTrapfile' : 'cdm_euclid_parallel.dat', 'serialTrapfile' : 'cdm_euclid_serial.dat', + 'beta_s' : 0.6, 'beta_p': 0.6, 'fwc' : 90000, 'vth' : 1.168e7, 't' : 20.48e-3, 'vg' : 6.e-11, + 'st' : 5.0e-6, 'sfwc' : 730000., 'svg' : 1.0e-10}) + return out + + +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 + + +if __name__ == '__main__': + print('IFSinstrument') diff --git a/csst_ifs_sim/support/__init__.py b/csst_ifs_sim/support/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/csst_ifs_sim/support/cosmicrays.py b/csst_ifs_sim/support/cosmicrays.py new file mode 100644 index 0000000000000000000000000000000000000000..e6c7464b7529deab9442071625632a9a7a7547fa --- /dev/null +++ b/csst_ifs_sim/support/cosmicrays.py @@ -0,0 +1,424 @@ +""" +Cosmic Rays +=========== + +This simple class can be used to include cosmic ray events to an image. +By default the cosmic ray events are drawn from distributions describing +the length and energy of the events. Such distributions can be generated +for example using Stardust code (http://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=04636917). +The energy of the cosmic ray events can also be set to constant for +testing purposes. The class can be used to draw a single cosmic ray +event or up to a covering fraction. + +:requires: NumPy +:requires: SciPy + +:version: 0.2 + +""" +import numpy as np +from scipy.interpolate import InterpolatedUnivariateSpline + + +class cosmicrays(): + """ + Cosmic ray generation class. Can either draw events from distributions or + set the energy of the events to a constant. + + :param log: logger instance + :param image: image to which cosmic rays are added to (a copy is made not to change the original numpy array) + :param crInfo: column information (cosmic ray file) + :param information: cosmic ray track information (file containing track length and energy information) and + exposure time. + """ + def __init__(self, log, image, crInfo=None, information=None): + """ + Cosmic ray generation class. Can either draw events from distributions or + set the energy of the events to a constant. + + :param log: logger instance + :param image: image to which cosmic rays are added to (a copy is made not to change the original numpy array) + :param crInfo: column information (cosmic ray file) + :param information: cosmic ray track information (file containing track length and energy information) and + exposure time. + """ + #setup logger + self.log = log + + #image and size + self.image = image.copy() + self.ysize, self.xsize = self.image.shape + + #set up the information dictionary, first with defaults and then overwrite with inputs if given + self.information = (dict(cosmicraylengths='/home/yan/csst-master/data/cdf_cr_length.dat', + cosmicraydistance='/home/yan/csst-master/data/cdf_cr_total.dat', + exptime=565)) + if information is not None: + self.information.update(information) + + if crInfo is not None: + self.cr = crInfo + else: + self._readCosmicrayInformation() + + + def _readCosmicrayInformation(self): + self.log.info('Reading in cosmic ray information from %s and %s' % (self.information['cosmicraylengths'], + self.information['cosmicraydistance'])) + #read in the information from the files + crLengths = np.loadtxt(self.information['cosmicraylengths']) + crDists = np.loadtxt(self.information['cosmicraydistance']) + + #set up the cosmic ray information dictionary + 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 self.cr + + + def _cosmicRayIntercepts(self, lum, x0, y0, l, phi): + """ + Derive cosmic ray streak intercept points. + + :param lum: luminosities of the cosmic ray tracks + :param x0: central positions of the cosmic ray tracks in x-direction + :param y0: central positions of the cosmic ray tracks in y-direction + :param l: lengths of the cosmic ray tracks + :param phi: orientation angles of the cosmic ray tracks + + :return: cosmic ray map (image) + :rtype: nd-array + """ + #create empty array + crImage = np.zeros((self.ysize, self.xsize), dtype=np.float64) + + #x and y shifts + dx = l * np.cos(phi) / 2. + dy = l * np.sin(phi) / 2. + mskdx = np.abs(dx) < 1e-8 + mskdy = np.abs(dy) < 1e-8 + dx[mskdx] = 0. + dy[mskdy] = 0. + + #pixels in x-direction + ilo = np.round(x0.copy() - dx) + msk = ilo < 0. + ilo[msk] = 0 + ilo = ilo.astype(np.int) + + ihi = 1 + np.round(x0.copy() + dx) + msk = ihi > self.xsize + ihi[msk] = self.xsize + ihi = ihi.astype(np.int) + + #pixels in y-directions + jlo = np.round(y0.copy() - dy) + msk = jlo < 0. + jlo[msk] = 0 + jlo = jlo.astype(np.int) + + jhi = 1 + np.round(y0.copy() + dy) + msk = jhi > self.ysize + jhi[msk] = self.ysize + jhi = jhi.astype(np.int) + + #loop over the individual events + for i, luminosity in enumerate(lum): + n = 0 # count the intercepts + + u = [] + x = [] + y = [] + + #Compute X intercepts on the pixel grid + if ilo[i] < ihi[i]: + for xcoord in range(ilo[i], ihi[i]): + ok = (xcoord - x0[i]) / dx[i] + if np.abs(ok) <= 0.5: + n += 1 + u.append(ok) + x.append(xcoord) + y.append(y0[i] + ok * dy[i]) + else: + for xcoord in range(ihi[i], ilo[i]): + ok = (xcoord - x0[i]) / dx[i] + if np.abs(ok) <= 0.5: + n += 1 + u.append(ok) + x.append(xcoord) + y.append(y0[i] + ok * dy[i]) + + #Compute Y intercepts on the pixel grid + if jlo[i] < jhi[i]: + for ycoord in range(jlo[i], jhi[i]): + ok = (ycoord - y0[i]) / dy[i] + if np.abs(ok) <= 0.5: + n += 1 + u.append(ok) + x.append(x0[i] + ok * dx[i]) + y.append(ycoord) + else: + for ycoord in range(jhi[i], jlo[i]): + ok = (ycoord - y0[i]) / dy[i] + if np.abs(ok) <= 0.5: + n += 1 + u.append(ok) + x.append(x0[i] + ok * dx[i]) + y.append(ycoord) + + #check if no intercepts were found + if n < 1: + xc = int(np.floor(x0[i])) + yc = int(np.floor(y0[i])) + crImage[yc, xc] += luminosity + + #Find the arguments that sort the intersections along the track + u = np.asarray(u) + x = np.asarray(x) + y = np.asarray(y) + + args = np.argsort(u) + + u = u[args] + x = x[args] + y = y[args] + + #Decide which cell each interval traverses, and the path length + for i in range(1, n - 1): + w = u[i + 1] - u[i] + cx = int(1 + np.floor((x[i + 1] + x[i]) / 2.)) + cy = int(1 + np.floor((y[i + 1] + y[i]) / 2.)) + + if 0 <= cx < self.xsize and 0 <= cy < self.ysize: + crImage[cy, cx] += (w * luminosity) + + return crImage + + + def _drawCosmicRays(self, limit=None): + """ + 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. + """ + #estimate the number of cosmics + cr_n = self.xsize * self.ysize * 0.014 / 43.263316 * 2. + #scale with exposure time, the above numbers are for the nominal 565s exposure + cr_n *= (self.information['exptime'] / 565.0) + + #assume a power-law intensity distribution for tracks + fit = dict(cr_lo=1.0e3, cr_hi=1.0e5, cr_q=2.0e0) + fit['q1'] = 1.0e0 - fit['cr_q'] + fit['en1'] = fit['cr_lo'] ** fit['q1'] + fit['en2'] = fit['cr_hi'] ** fit['q1'] + + #pseudo-random numbers taken from a uniform distribution between 0 and 1 + np.random.seed() + luck = np.random.rand(int(np.floor(cr_n))) + + #draw the length of the tracks + if self.cr['cr_cdfn'] > 1: + ius = InterpolatedUnivariateSpline(self.cr['cr_cdf'], self.cr['cr_u']) + self.cr['cr_l'] = ius(luck) + else: + self.cr['cr_l'] = np.sqrt(1.0 - luck ** 2) / luck + + #draw the energy of the tracks + if self.cr['cr_cden'] > 1: + ius = InterpolatedUnivariateSpline(self.cr['cr_cde'], self.cr['cr_v']) + self.cr['cr_e'] = ius(luck) + else: + np.random.seed() + self.cr['cr_e'] = (fit['en1'] + (fit['en2'] - fit['en1']) * + np.random.rand(int(np.floor(cr_n)))) ** (1.0 / fit['q1']) + + #Choose the properties such as positions and an angle from a random Uniform dist + np.random.seed() + cr_x = self.xsize * np.random.rand(int(np.floor(cr_n))) + + np.random.seed() + cr_y = self.ysize * np.random.rand(int(np.floor(cr_n))) + + np.random.seed() + cr_phi = np.pi * np.random.rand(int(np.floor(cr_n))) + + #find the intercepts + if limit is None: + self.cosmicrayMap = self._cosmicRayIntercepts(self.cr['cr_e'], cr_x, cr_y, self.cr['cr_l'], cr_phi) + print ('Number of cosmic ray events:', len(self.cr['cr_e'])) + else: + #limit to electron levels < limit + msk = self.cr['cr_e'] < limit + print ('Number of cosmic ray events: %i / %i' % (len(self.cr['cr_e'][msk]), int(np.floor(cr_n)))) + self.cosmicrayMap = self._cosmicRayIntercepts(self.cr['cr_e'][msk], cr_x[msk], cr_y[msk], + self.cr['cr_l'][msk], cr_phi[msk]) + + #count the covering factor + area_cr = np.count_nonzero(self.cosmicrayMap) + text = 'The cosmic ray covering factor is %i pixels i.e. %.3f per cent' \ + % (area_cr, 100.*area_cr / (self.xsize*self.ysize)) + self.log.info(text) + print (text) + + + def _drawSingleEvent(self, limit=1000, cr_n=1): + """ + Generate a single cosmic ray event and include it to a cosmic ray map (self.cosmicrayMap). + + :param limit: limiting energy for the cosmic ray event + :type limit: float + :param cr_n: number of cosmic ray events to include + :type cr_n: int + + :return: None + """ + #pseudo-random numbers taken from a uniform distribution between 0 and 1 + np.random.seed() + luck = np.random.rand(cr_n) + + #draw the length of the tracks + ius = InterpolatedUnivariateSpline(self.cr['cr_cdf'], self.cr['cr_u']) + self.cr['cr_l'] = ius(luck) + + #set the energy directly to the limit + self.cr['cr_e'] = np.asarray([limit, ]*cr_n) + + #Choose the properties such as positions and an angle from a random Uniform dist + np.random.seed() + cr_x = self.xsize * np.random.rand(int(np.floor(cr_n))) + + np.random.seed() + cr_y = self.ysize * np.random.rand(int(np.floor(cr_n))) + + np.random.seed() + cr_phi = np.pi * np.random.rand(int(np.floor(cr_n))) + + #find the intercepts + self.cosmicrayMap = self._cosmicRayIntercepts(self.cr['cr_e'], cr_x, cr_y, self.cr['cr_l'], cr_phi) + + #count the covering factor + area_cr = np.count_nonzero(self.cosmicrayMap) + text = 'The cosmic ray covering factor is %i pixels i.e. %.3f per cent' \ + % (area_cr, 100.*area_cr / (self.xsize*self.ysize)) + self.log.info(text) + print( text) + + + def _drawEventsToCoveringFactor(self, coveringFraction=3.0, limit=1000, verbose=False): + """ + Generate cosmic ray events up to a covering fraction and include it to a cosmic ray map (self.cosmicrayMap). + + :param coveringFraction: covering fraction of cosmic rya events in per cent of total number of pixels + :type coveringFraction: float + :param limit: limiting energy for the cosmic ray event [None = draw from distribution] + :type limit: None or float + :param verbose: print out information to stdout + :type verbose: bool + + + :return: None + """ + self.cosmicrayMap = np.zeros((self.ysize, self.xsize)) + + #how many events to draw at once, too large number leads to exceeding the covering fraction + cr_n = int(295 * self.information['exptime'] / 565. * coveringFraction / 1.4) + + covering = 0.0 + + while covering < coveringFraction: + #pseudo-random numbers taken from a uniform distribution between 0 and 1 + np.random.seed() + luck = np.random.rand(cr_n) + + #draw the length of the tracks + ius = InterpolatedUnivariateSpline(self.cr['cr_cdf'], self.cr['cr_u']) + self.cr['cr_l'] = ius(luck) + + if limit is None: + ius = InterpolatedUnivariateSpline(self.cr['cr_cde'], self.cr['cr_v']) + self.cr['cr_e'] = ius(luck) + else: + #set the energy directly to the limit + self.cr['cr_e'] = np.asarray([limit,]) + + #Choose the properties such as positions and an angle from a random Uniform dist + np.random.seed() + cr_x = self.xsize * np.random.rand(int(np.floor(cr_n))) + + np.random.seed() + cr_y = self.ysize * np.random.rand(int(np.floor(cr_n))) + + np.random.seed() + cr_phi = np.pi * np.random.rand(int(np.floor(cr_n))) + + #find the intercepts + self.cosmicrayMap += self._cosmicRayIntercepts(self.cr['cr_e'], cr_x, cr_y, self.cr['cr_l'], cr_phi) + + #count the covering factor + area_cr = np.count_nonzero(self.cosmicrayMap) + covering = 100.*area_cr / (self.xsize*self.ysize) + + text = 'The cosmic ray covering factor is %i pixels i.e. %.3f per cent' % (area_cr, covering) + self.log.info(text) + + if verbose: + print( text) + + + def addCosmicRays(self, limit=None): + """ + Include cosmic rays to the image given. + + :return: image with cosmic rays + :rtype: ndarray + """ + self._drawCosmicRays(limit=limit) + + #paste cosmic rays + self.image += self.cosmicrayMap + + return self.image + + + def addSingleEvent(self, limit=None): + """ + Include a single cosmic ray event to the image given. + + :return: image with cosmic rays + :rtype: ndarray + """ + self._drawSingleEvent(limit=limit) + + #paste cosmic rays + self.image += self.cosmicrayMap + + return self.image + + + def addUpToFraction(self, coveringFraction, limit=None, verbose=False): + """ + Add cosmic ray events up to the covering Fraction. + + :param coveringFraction: covering fraction of cosmic rya events in per cent of total number of pixels + :type coveringFraction: float + :param limit: limiting energy for the cosmic ray event [None = draw from distribution] + :type limit: None or float + :param verbose: print out information to stdout + :type verbose: bool + + :return: image with cosmic rays + :rtype: ndarray + """ + self._drawEventsToCoveringFactor(coveringFraction, limit=limit, verbose=verbose) + + #paste cosmic rays + self.image += self.cosmicrayMap + + return self.image + + +if __name__ == "__main__": + + + print() + diff --git a/csst_ifs_sim/support/logger.py b/csst_ifs_sim/support/logger.py new file mode 100644 index 0000000000000000000000000000000000000000..b09ef822fe697c173e3d8d4b1970fe6af1e6f39a --- /dev/null +++ b/csst_ifs_sim/support/logger.py @@ -0,0 +1,40 @@ +""" +These functions can be used for logging information. + +.. Warning:: logger is not multiprocessing safe. + +:version: 0.3 +""" +import logging +import logging.handlers + + +def setUpLogger(log_filename, loggername='logger'): + """ + Sets up a logger. + + :param: log_filename: name of the file to save the log. + :param: loggername: name of the logger + + :return: logger instance + """ + # create logger + logger = logging.getLogger(loggername) + logger.setLevel(logging.DEBUG) + # Add the log message handler to the logger + handler = logging.handlers.RotatingFileHandler(log_filename) + #maxBytes=20, backupCount=5) + # create formatter + formatter = logging.Formatter('%(asctime)s - %(module)s - %(funcName)s - %(levelname)s - %(message)s') + # add formatter to ch + handler.setFormatter(formatter) + # add handler to logger + + if (logger.hasHandlers()): + logger.handlers.clear() + + logger.addHandler(handler) + + return logger + + diff --git a/ifs_data/__init__.py b/ifs_data/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/setup.py b/setup.py new file mode 100644 index 0000000000000000000000000000000000000000..f65e2b384bdc95954408133449d4fdc1562b5262 --- /dev/null +++ b/setup.py @@ -0,0 +1,36 @@ +import setuptools + +with open("README.md", "r") as fh: + long_description = fh.read() + +setuptools.setup( + name='csst_ifs_sim', + version='2.0.0-IFS1.0.0', + author='CSST Team', + author_email='zhaojunyan@shao.ac.cn', + description='The CSST - ifs - sim', # short description + long_description=long_description, + long_description_content_type="text/markdown", + url='https://csst-tb.bao.ac.cn/', + # project_urls={ + # 'Source': 'https://csst-tb.bao.ac.cn/code/csst-l1/ifs/csst_ifs_sim', + # }, + packages=setuptools.find_packages(), + license='MIT', + classifiers=["Development Status :: 5 - Production/Stable", + "Intended Audience :: Science/Research", + "License :: OSI Approved :: MIT License", + "Operating System :: OS Independent", + "Programming Language :: Python :: 3.8", + "Topic :: Scientific/Engineering :: Physics", + "Topic :: Scientific/Engineering :: Astronomy"], + package_dir={'csst_ifs_sim': 'csst_ifs_sim'}, + # include_package_data=True, + package_data={"": ["LICENSE", "README.md"], + "csst_ifs_sim": ["ifs_so/*", "ifs_data/*", "ifs_data/refs/*", "ifs_data/refs/orbit20160925/*"]}, + # install_requires=['sphinx', + # 'numpy', + # 'scipy', 'matplotlib', + # 'astropy', 'healpy', 'ccdproc', 'deepCR', 'photutils'], + python_requires='>=3.8', +)