test

csst_mci_sim/CTI/CTI.py

@@ -159,6 +159,8 @@ class 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

@@ -71,7 +71,7 @@ from scipy import ndimage
 
 
 #sys.path.append('./csst_mci_sim')
-from .CTI import CTI
+#from .CTI import CTI
 from .support import logger as lg
 from .support import cosmicrays
 from .support import shao
@@ -80,7 +80,6 @@ from .support import MCIinstrumentModel
 from .mci_so import cdm03bidir
 
 
-
 from joblib import Parallel, delayed
 from astropy.coordinates import SkyCoord
 from scipy import interpolate
@@ -91,6 +90,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)
@@ -2599,7 +2772,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.')
@@ -2611,7 +2784,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.')
@@ -2622,7 +2795,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.')