debug

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'],

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
+                             
+                            sum += overload
+                        elif sum > 0.:
+                            if -overload > sum:
+                                overload = -sum
+                            #self.image[-j-1, i] -= overload
+                            data[-j-1, i] -= overload                
+                            sum += 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
-                        
-            ######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['result_path']    = result_path
+    
     
     sim[iLoop].simulate(iLoop)   
        

csst_mci_sim/mci_data/mci_all_9K.config

 [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

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
 

tests/test_mci_sim.py

@@ -49,7 +49,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,
@@ -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,