diff --git a/csst_mci_sim/CTI/CTI.py b/csst_mci_sim/CTI/CTI.py index 13e69410474ec6ce2a5817cac97b79824dd4dd09..3abb3d88cbbfaab33216deb322b330808211dcd5 100644 --- a/csst_mci_sim/CTI/CTI.py +++ b/csst_mci_sim/CTI/CTI.py @@ -156,9 +156,11 @@ class CDM03bidir(): ################################################################################# ###modify #sys.path.append('../so') - from mci_so import cdm03bidir + # from ifs_so.cdm03.cpython-38-x86_64-linux-gnu import cdm03bidir # import cdm03bidir + from .mci_so import cdm03bidir + CTIed = cdm03bidir.cdm03(np.asfortranarray(data), jflip, iflip, self.values['dob'], self.values['rdose'], diff --git a/csst_mci_sim/csst_mci_sim.py b/csst_mci_sim/csst_mci_sim.py index 9221c10267dcc56ca564ae201545b21c9a2f5a18..2afefc3436ae4b8b51248ab1d444d7e046a45659 100644 --- a/csst_mci_sim/csst_mci_sim.py +++ b/csst_mci_sim/csst_mci_sim.py @@ -66,15 +66,17 @@ from astropy.io import fits from astropy import units as u import os, sys, math import configparser as ConfigParser -from matplotlib import pyplot as plt +#from matplotlib import pyplot as plt + from scipy import ndimage -sys.path.append('./csst_mci_sim') -from CTI import CTI from support import logger as lg from support import cosmicrays from support import shao from support import sed from support import MCIinstrumentModel +from mci_so import cdm03bidir + + from joblib import Parallel, delayed from astropy.coordinates import SkyCoord from scipy import interpolate @@ -85,6 +87,180 @@ import astropy.coordinates as coord from scipy.interpolate import interp1d ########################### functions ######################### + +""" +Charge Transfer Inefficiency +============================ + +This file contains a simple class to run a CDM03 CTI model developed by Alex Short (ESA). + +This now contains both the official CDM03 and a new version that allows different trap +parameters in parallel and serial direction. + +:requires: NumPy +:requires: CDM03 (FORTRAN code, f2py -c -m cdm03bidir cdm03bidir.f90) + +:version: 0.35 +""" +import numpy as np + + +#CDM03bidir +class CDM03bidir(): + """ + Class to run CDM03 CTI model, class Fortran routine to perform the actual CDM03 calculations. + + :param settings: input parameters + :type settings: dict + :param data: input data to be radiated + :type data: ndarray + :param log: instance to Python logging + :type log: logging instance + """ + def __init__(self, settings, data, log=None): + """ + Class constructor. + + :param settings: input parameters + :type settings: dict + :param data: input data to be radiated + :type data: ndarray + :param log: instance to Python logging + :type log: logging instance + """ + self.data = data + self.values = dict(quads=(0,1,2,3), xsize=2048, ysize=2066, dob=0.0, rdose=8.0e9) + self.values.update(settings) + self.log = log + self._setupLogger() + + #default CDM03 settings + self.params = dict(beta_p=0.6, beta_s=0.6, fwc=200000., vth=1.168e7, vg=6.e-11, t=20.48e-3, + sfwc=730000., svg=1.0e-10, st=5.0e-6, parallel=1., serial=0.0) + #update with inputs + self.params.update(self.values) + + #read in trap information + trapdata = np.loadtxt(self.values['dir_path']+self.values['paralleltrapfile']) + if trapdata.ndim > 1: + self.nt_p = trapdata[:, 0] + self.sigma_p = trapdata[:, 1] + self.taur_p = trapdata[:, 2] + else: + #only one trap species + self.nt_p = [trapdata[0],] + self.sigma_p = [trapdata[1],] + self.taur_p = [trapdata[2],] + + trapdata = np.loadtxt(self.values['dir_path']+self.values['serialtrapfile']) + if trapdata.ndim > 1: + self.nt_s = trapdata[:, 0] + self.sigma_s = trapdata[:, 1] + self.taur_s = trapdata[:, 2] + else: + #only one trap species + self.nt_s = [trapdata[0],] + self.sigma_s = [trapdata[1],] + self.taur_s = [trapdata[2],] + + #scale thibaut's values + if 'thibaut' in self.values['parallelTrapfile']: + self.nt_p /= 0.576 #thibaut's values traps / pixel + self.sigma_p *= 1.e4 #thibaut's values in m**2 + if 'thibaut' in self.values['serialTrapfile']: + self.nt_s *= 0.576 #thibaut's values traps / pixel #should be division? + self.sigma_s *= 1.e4 #thibaut's values in m**2 + + + def _setupLogger(self): + """ + Set up the logger. + """ + self.logger = True + # if self.log is None: + # self.logger = False + + + def applyRadiationDamage(self, data, iquadrant=0): + """ + Apply radian damage based on FORTRAN CDM03 model. The method assumes that + input data covers only a single quadrant defined by the iquadrant integer. + + :param data: imaging data to which the CDM03 model will be applied to. + :type data: ndarray + + :param iquandrant: number of the quadrant to process + :type iquandrant: int + + cdm03 - Function signature:: + + sout = cdm03(sinp,iflip,jflip,dob,rdose,in_nt,in_sigma,in_tr,[xdim,ydim,zdim]) + Required arguments: + sinp : input rank-2 array('d') with bounds (xdim,ydim) + iflip : input int + jflip : input int + dob : input float + rdose : input float + in_nt : input rank-1 array('d') with bounds (zdim) + in_sigma : input rank-1 array('d') with bounds (zdim) + in_tr : input rank-1 array('d') with bounds (zdim) + Optional arguments: + xdim := shape(sinp,0) input int + ydim := shape(sinp,1) input int + zdim := len(in_nt) input int + Return objects: + sout : rank-2 array('d') with bounds (xdim,ydim) + + .. Note:: Because Python/NumPy arrays are different row/column based, one needs + to be extra careful here. NumPy.asfortranarray will be called to get + an array laid out in Fortran order in memory. Before returning the + array will be laid out in memory in C-style (row-major order). + + :return: image that has been run through the CDM03 model + :rtype: ndarray + """"" + #return data + + iflip = iquadrant / 2 + jflip = iquadrant % 2 + + params = [self.params['beta_p'], self.params['beta_s'], self.params['fwc'], self.params['vth'], + self.params['vg'], self.params['t'], self.params['sfwc'], self.params['svg'], self.params['st'], + self.params['parallel'], self.params['serial']] + + if self.logger: + self.log.info('nt_p=' + str(self.nt_p)) + self.log.info('nt_s=' + str(self.nt_s)) + self.log.info('sigma_p= ' + str(self.sigma_p)) + self.log.info('sigma_s= ' + str(self.sigma_s)) + self.log.info('taur_p= ' + str(self.taur_p)) + self.log.info('taur_s= ' + str(self.taur_s)) + self.log.info('dob=%f' % self.values['dob']) + self.log.info('rdose=%e' % self.values['rdose']) + self.log.info('xsize=%i' % data.shape[1]) + self.log.info('ysize=%i' % data.shape[0]) + self.log.info('quadrant=%i' % iquadrant) + self.log.info('iflip=%i' % iflip) + self.log.info('jflip=%i' % jflip) + +################################################################################# + + + CTIed = cdm03bidir.cdm03(np.asfortranarray(data), + jflip, iflip, + self.values['dob'], self.values['rdose'], + self.nt_p, self.sigma_p, self.taur_p, + self.nt_s, self.sigma_s, self.taur_s, + params, + [data.shape[0], data.shape[1], len(self.nt_p), len(self.nt_s), len(self.params)]) + + + return np.asanyarray(CTIed) + +################################################################################################################# +################################################################################################################# + + def transRaDec2D(ra, dec): # radec转为竞天程序里的ob, 赤道坐标系下的笛卡尔三维坐标xyz. x1 = np.cos(dec / 57.2957795) * np.cos(ra / 57.2957795) @@ -835,7 +1011,7 @@ class MCIsimulator(): ############################################################################### ############################################################################### - def configure(self,simnumber,sourcein,dir_path): + def configure(self,simnumber,sourcein,dir_path,result_path): """ Configures the simulator with input information and creates and empty array to which the final image will be build on. @@ -848,6 +1024,8 @@ class MCIsimulator(): self.information['dir_path']=dir_path + self.information['result_path']=result_path + self.source=sourcein ##print('print information:', self.information) @@ -858,18 +1036,21 @@ class MCIsimulator(): #data_time=now.strftime("%Y-%m-%d-%H-%M-%S") result_day=now.strftime("%Y-%m-%d") - if self.information['dir_path'] =='/nfsdata/share/simulation-unittest/mci_sim/': - self.result_path = self.information['dir_path'] +'mci_sim_result/'+result_day - else: + # if self.information['dir_path'] =='/nfsdata/share/simulation-unittest/mci_sim/': + # self.result_path = self.information['dir_path'] +'mci_sim_result/'+result_day + # else: - home_path = os.environ['HOME'] + # home_path = os.environ['HOME'] - if home_path == '/home/yan': + # if home_path == '/home/yan': - self.result_path = '../MCI_simResult/'+self.source+"_"+result_day - else: - self.result_path = '/data/mcisimdata/'+result_day - + # self.result_path = '../MCI_simResult/'+self.source+"_"+result_day + # else: + # self.result_path = '/data/mcisimdata/'+result_day + + print(self.information['result_path']) + + self.result_path = self.information['result_path']+self.source+"_"+result_day if os.path.isdir(self.result_path)==False: os.mkdir(self.result_path) @@ -2593,7 +2774,7 @@ class MCIsimulator(): 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) + cti = CDM03bidir(self.information, [], log=self.log) #here we need the right input data self.image_g = cti.applyRadiationDamage(self.image_g.copy().transpose(), iquadrant=self.information['quadrant']).transpose() self.log.info('Radiation damage added.') @@ -2605,7 +2786,7 @@ class MCIsimulator(): 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) + 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.') @@ -2616,7 +2797,7 @@ class MCIsimulator(): 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) + cti = CDM03bidir(self.information, [], log=self.log) #here we need the right input data self.image_i = cti.applyRadiationDamage(self.image_i.copy().transpose(), iquadrant=self.information['quadrant']).transpose() self.log.info('Radiation damage added.') @@ -2809,7 +2990,7 @@ class MCIsimulator(): #################################################################################### - def applyBleeding(self, img, direction='horizon'): + def applyBleeding(self, img, direction='not_horizon'): """ Apply bleeding along the CCD readout direction if the number of electrons in a pixel exceeds the full-well capacity. @@ -2867,45 +3048,45 @@ class MCIsimulator(): data[i,-j-1,] -= overload sum += overload - # else: + else: - # #loop over each column, as bleeding is modelled column-wise - # for i, column in enumerate(data.T): - # sum = 0. - # for j, value in enumerate(column): - # #first round - from bottom to top (need to half the bleeding) - # overload = value - self.information['fullwellcapacity'] - # if overload > 0.: - # overload /= 2. - # #self.image[j, i] -= overload - # data[j, i] -= overload - - # sum += overload - # elif sum > 0.: - # if -overload > sum: - # overload = -sum - # #self.image[j, i] -= overload - # data[j, i] -= overload - # sum += overload - # ################################ - # for i, column in enumerate(data.T): - # sum = 0. - # for j, value in enumerate(column[::-1]): - # #second round - from top to bottom (bleeding was half'd already, so now full) - # overload = value - self.information['fullwellcapacity'] - # if overload > 0.: - # #self.image[-j-1, i] -= overload - # data[-j-1, i] -= overload + #loop over each column, as bleeding is modelled column-wise + for i, column in enumerate(data.T): + sum = 0. + for j, value in enumerate(column): + #first round - from bottom to top (need to half the bleeding) + overload = value - self.information['fullwellcapacity'] + if overload > 0.: + overload /= 2. + #self.image[j, i] -= overload + data[j, i] -= overload + + sum += overload + elif sum > 0.: + if -overload > sum: + overload = -sum + #self.image[j, i] -= overload + data[j, i] -= overload + sum += overload + ################################ + for i, column in enumerate(data.T): + sum = 0. + for j, value in enumerate(column[::-1]): + #second round - from top to bottom (bleeding was half'd already, so now full) + overload = value - self.information['fullwellcapacity'] + if overload > 0.: + #self.image[-j-1, i] -= overload + data[-j-1, i] -= overload - # sum += overload - # elif sum > 0.: - # if -overload > sum: - # overload = -sum - # #self.image[-j-1, i] -= overload - # data[-j-1, i] -= overload - # sum += overload + sum += overload + elif sum > 0.: + if -overload > sum: + overload = -sum + #self.image[-j-1, i] -= overload + data[-j-1, i] -= overload + sum += overload - ######print('Applying column bleeding finished.......') + #####print('Applying column bleeding finished.......') return data ############################################################################ @@ -4732,7 +4913,7 @@ class MCIsimulator(): ################################################################################################ -def runMCIsim(sourcein,configfile,dir_path, debug, iLoop): +def runMCIsim(sourcein,configfile,dir_path, result_path, debug, iLoop): print('Path Test：dir_path', dir_path) @@ -4740,10 +4921,12 @@ def runMCIsim(sourcein,configfile,dir_path, debug, iLoop): sim= dict() sim[iLoop] = MCIsimulator(configfile) - sim[iLoop].configure(iLoop,sourcein,dir_path) # load the configfile; + sim[iLoop].configure(iLoop,sourcein,dir_path,result_path) # load the configfile; sim[iLoop].information['sourcein'] =sourcein - sim[iLoop].information['debug'] =debug + sim[iLoop].information['debug'] =debug + sim[iLoop].information['result_path'] = result_path + sim[iLoop].simulate(iLoop) diff --git a/csst_mci_sim/help/CSST-MCI-Cycle9.pdf b/csst_mci_sim/help/CSST-MCI-Cycle9.pdf new file mode 100644 index 0000000000000000000000000000000000000000..780fad013c3cf6d36e398a3c22f47b20c735b29a Binary files /dev/null and b/csst_mci_sim/help/CSST-MCI-Cycle9.pdf differ diff --git a/csst_mci_sim/mci_data/mci_all_9K.config b/csst_mci_sim/mci_data/mci_all_9K.config index 06806bdc6b2576246ca257b012eb75d161e7728f..789cafc19f6b656f2334988e74c8d9e3b49a29e5 100755 --- a/csst_mci_sim/mci_data/mci_all_9K.config +++ b/csst_mci_sim/mci_data/mci_all_9K.config @@ -1,9 +1,5 @@ [TEST] -dir_path=mci_sim/MCI_inputData/ - -result_path=mci_sim/mci_sim_result/ - #size of the output image array, xsize is column, ysize is row, xsize = 9216,ysize = 9232 xsize =9216 @@ -79,9 +75,9 @@ sim_star = yes sim_galaxy = yes -save_starpsf = yes +save_starpsf = no -save_cosmicrays = yes +save_cosmicrays = no ############################################## ############################################## @@ -92,7 +88,7 @@ fullwellcapacity = 90000 dark = 0.001 #exposure to simulate, exposure time -exptime = 100.0 +exptime = 300.0 ###PNRU matrix sigma flatsigma=0.001 diff --git a/csst_mci_sim/mci_so/__pycache__/__init__.cpython-311.pyc b/csst_mci_sim/mci_so/__pycache__/__init__.cpython-311.pyc new file mode 100644 index 0000000000000000000000000000000000000000..ee21e1d4e6a56d65b280e1a4797d2e030b258fd5 Binary files /dev/null and b/csst_mci_sim/mci_so/__pycache__/__init__.cpython-311.pyc differ diff --git a/csst_mci_sim/support/__pycache__/MCIinstrumentModel.cpython-311.pyc b/csst_mci_sim/support/__pycache__/MCIinstrumentModel.cpython-311.pyc new file mode 100644 index 0000000000000000000000000000000000000000..d34909d1f594354d1de22bf15b997700ec0cfdc0 Binary files /dev/null and b/csst_mci_sim/support/__pycache__/MCIinstrumentModel.cpython-311.pyc differ diff --git a/csst_mci_sim/support/__pycache__/__init__.cpython-311.pyc b/csst_mci_sim/support/__pycache__/__init__.cpython-311.pyc new file mode 100644 index 0000000000000000000000000000000000000000..71b21f61bdb1a1a9db5f769739a1c4615ec21e6b Binary files /dev/null and b/csst_mci_sim/support/__pycache__/__init__.cpython-311.pyc differ diff --git a/csst_mci_sim/support/__pycache__/cosmicrays.cpython-311.pyc b/csst_mci_sim/support/__pycache__/cosmicrays.cpython-311.pyc new file mode 100644 index 0000000000000000000000000000000000000000..f7de180a574374496317402a0311d36d2d4c1681 Binary files /dev/null and b/csst_mci_sim/support/__pycache__/cosmicrays.cpython-311.pyc differ diff --git a/csst_mci_sim/support/__pycache__/logger.cpython-311.pyc b/csst_mci_sim/support/__pycache__/logger.cpython-311.pyc new file mode 100644 index 0000000000000000000000000000000000000000..36380654cf82d1ab0b8ad24e9810738b8d53be53 Binary files /dev/null and b/csst_mci_sim/support/__pycache__/logger.cpython-311.pyc differ diff --git a/csst_mci_sim/support/__pycache__/sed.cpython-311.pyc b/csst_mci_sim/support/__pycache__/sed.cpython-311.pyc new file mode 100644 index 0000000000000000000000000000000000000000..817b124eb28ff644f99b716f7c2244c0a40f2aa9 Binary files /dev/null and b/csst_mci_sim/support/__pycache__/sed.cpython-311.pyc differ diff --git a/csst_mci_sim/support/__pycache__/shao.cpython-311.pyc b/csst_mci_sim/support/__pycache__/shao.cpython-311.pyc new file mode 100644 index 0000000000000000000000000000000000000000..40873b02fb010fdbb9de69c408974963f8663252 Binary files /dev/null and b/csst_mci_sim/support/__pycache__/shao.cpython-311.pyc differ diff --git a/csst_mci_sim/support/cosmicrays.py b/csst_mci_sim/support/cosmicrays.py index c1386fca517a4e7437d8f245e706607a374913eb..161be49d2c7df66e8668c786143deb2af4cf41d2 100644 --- a/csst_mci_sim/support/cosmicrays.py +++ b/csst_mci_sim/support/cosmicrays.py @@ -199,7 +199,9 @@ class cosmicrays(): 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.exptime / 565. * coveringFraction / 1.4) + ####cr_n = int(295 * self.exptime / 565. * coveringFraction / 1.4) + + cr_n = int(5000 * self.exptime / 565. * coveringFraction) covering = 0.0 @@ -236,8 +238,8 @@ class cosmicrays(): 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) + text = 'The cosmic ray covering factor is %i pixels i.e. %.3f per cent' % (area_cr, covering) + self.log.info(text) ###################################################33 diff --git a/tests/test_mci_sim.py b/tests/test_mci_sim.py index 95a644feaea5086f12c26a90ff08f8a6ed072585..b840ba10d6ad2f052cf5be9f67b1dc598f46461a 100644 --- a/tests/test_mci_sim.py +++ b/tests/test_mci_sim.py @@ -48,8 +48,10 @@ class TestDemoFunction(unittest.TestCase): print(configfile) debug=True + + result_path = dir_path +'mci_sim_result/' - csst_mci_sim.runMCIsim(sourcein, configfile, dir_path, debug, 1) + csst_mci_sim.runMCIsim(sourcein, configfile, dir_path, result_path, debug, 1) self.assertEqual( 1 , 1, @@ -89,7 +91,9 @@ class TestDemoFunction(unittest.TestCase): debug=True - csst_mci_sim.runMCIsim(sourcein, configfile, dir_path, debug, 1) + result_path = dir_path +'mci_sim_result/' + + csst_mci_sim.runMCIsim(sourcein, configfile, dir_path, result_path, debug, 1) self.assertEqual( 1 , 1, @@ -130,7 +134,9 @@ class TestDemoFunction(unittest.TestCase): debug=True - csst_mci_sim.runMCIsim(sourcein, configfile, dir_path, debug, 1) + result_path = dir_path +'mci_sim_result/' + + csst_mci_sim.runMCIsim(sourcein, configfile, dir_path, result_path, debug, 1) self.assertEqual( 1 , 1, @@ -170,7 +176,9 @@ class TestDemoFunction(unittest.TestCase): debug=True - csst_mci_sim.runMCIsim(sourcein, configfile, dir_path, debug, 1) + result_path = dir_path +'mci_sim_result/' + + csst_mci_sim.runMCIsim(sourcein, configfile, dir_path, result_path, debug, 1) self.assertEqual( 1 , 1, @@ -210,7 +218,9 @@ class TestDemoFunction(unittest.TestCase): debug=True - csst_mci_sim.runMCIsim(sourcein, configfile, dir_path, debug, 1) + result_path = dir_path +'mci_sim_result/' + + csst_mci_sim.runMCIsim(sourcein, configfile, dir_path, result_path, debug, 1) self.assertEqual( 1 , 1,