build

LICENSE

+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.

README.md

-# mci
+# 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/mci.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/mci/-/settings/integrations)
+- [ ] [Set up project integrations](http://10.3.10.28/code/csst-l1/ifs/csst_ifs_common/-/settings/integrations)
 
 ## Collaborate with your team
 

csst_mci_sim/CTI/CTI.py

+"""
+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)
+       
+
+#################################################################################################################
+      

csst_mci_sim/csst_mci_sim.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+
+@author: yan, zhaojun"""
+
+
+"""
+The CSST MCI Image Simulator
+=============================================
+
+This file contains an image simulator for the CSST MCI.
+
+The approximate sequence of events in the simulator is as follows:
+
+      #. Read in a configuration file, which defines for example,
+         detector characteristics (bias, dark and readout noise, gain,
+         plate scale and pixel scale, exposure time etc.).
+      #. Read in another file containing charge trap definitions (for CTI modelling).
+      #. Read in a file defining the cosmic rays (trail lengths and cumulative distributions).
+     
+      #. Loop over the number of exposures to co-add and for each object in the galaxy fits file:
+
+      #. Apply a multiplicative flat-field map to emulate pixel-to-pixel non-uniformity [optional].
+
+      #. Add cosmic ray tracks onto the CCD with random positions but known distribution [optional].
+      #. Apply detector charge bleeding in column direction [optional].
+      #. Add constant dark current  [optional].
+
+      #. Add photon (Poisson) noise [optional]
+      #. Add cosmetic defects from an input file [optional].
+      #. Add prescan and overscan regions in the serial direction [optional].
+      #. Apply the CDM03 radiation damage model [optional].
+      #. Apply non-linearity model to the pixel data [optional].
+      #. Add readout noise selected [optional].
+      #. Convert from electrons to ADUs using a given gain factor.
+      #. Add a given bias level and discretise the counts (the output is going to be in 16bit unsigned integers).
+      #. Finally the simulated image is converted to a FITS file, a WCS is assigned
+         and the output is saved to the current working directory.
+
+.. Warning:: The code is still work in progress and new features are being added.
+             The code has been tested, but nevertheless bugs may be lurking in corners, so
+             please report any weird or inconsistent simulations to the author.
+
+
+Contact Information: zhaojunyan@shao.ac.cn
+    
+-------
+
+"""
+import galsim
+import pandas as pd
+from scipy.integrate import simps
+import julian
+from datetime import datetime, timedelta
+from astropy.time import Time
+from astropy.coordinates import get_sun
+
+import numpy as np
+
+from tqdm import tqdm
+
+from astropy.wcs import WCS as WCS
+
+from  astropy.io import fits
+
+from astropy import units as u
+
+import os, sys, math 
+
+import configparser as ConfigParser
+
+from optparse import OptionParser
+
+from scipy import ndimage
+
+
+from astropy.coordinates import SkyCoord
+
+from scipy import interpolate
+
+from scipy.signal import fftconvolve
+
+
+
+sys.path.append('..')
+from CTI import CTI
+from support import logger as lg
+from support import cosmicrays
+from support import MCIinstrumentModel
+from support import shao 
+from support import sed 
+
+# set the folder
+FOLDER ='../'
+
+
+#####################################################################################
+
+
+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()
+
+
+###############################################################################
+
+def make_c_coor(fov, step):
+    """
+    Draw the mesh grids for a fov * fov box with given step .
+    """
+    
+    nc=int(fov/step)
+    
+    bs=fov
+    ds = bs / nc
+    xx01 = np.linspace(-bs / 2.0, bs / 2.0 - ds, nc) + 0.5 * ds
+    xx02 = np.linspace(-bs / 2.0, bs / 2.0 - ds, nc) + 0.5 * ds
+    xg2, xg1 = np.meshgrid(xx01, xx02)
+    return xg1, xg2
+
+##############################################################################    
+
+###############################################################################
+
+def str2time(strTime):
+    if len(strTime)>20:#
+        msec=int(float('0.'+strTime[20:])*1000000) # microsecond
+    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 cut_radius(rcutstar, mstar, mag):
+    return rcutstar * 10**(0.4*2*(mstar-mag)/3.125)
+
+def core_radius(rcorestar, mstar, mag):
+    return rcorestar * 10**(0.4*(mstar-mag)/2)
+
+def v_disp(sigstar, mstar, mag):
+    return sigstar * 10**(0.4*(mstar-mag)/3.57)
+
+###############################################################################
+#####################################################################################
+#### distortion function 
+
+def distortField(ra, dec, ch):
+  ###% ra ,dec are the idea position in arcsec , and the center position is ra=0, dec=0
+
+    '''MCI geometry distortion effect in channel G  R and I'''
+    distortA=dict()
+    
+    distortB=dict()      
+    
+    distortA['g']=np.array([0.000586224851462036,0.999778507355332,-0.000377391397200834,-4.40403573019393e-06,0.000147444528630966,-1.93281825050268e-06,5.73520238714447e-07,-1.05194725549731e-08,7.55124671251098e-07,-3.85623216175251e-09,-1.36254798567168e-11,1.36983654024573e-09,-4.53104347730230e-11,1.50641840169747e-09,-1.05226890727366e-11,1.06471556810228e-12,-1.42299485385165e-16,5.90008722878028e-12,2.36749977009538e-17,4.11061448730916e-12,3.39287277469805e-17])
+    distortB['g']=np.array([1.00113878750680,-0.000495107770623867,0.999162436209357,9.67138672907941e-05,-3.88808001621361e-06,0.000267940195644359,-3.03504629416955e-09,7.59508802931395e-07,-1.13952684052500e-08,9.58672893217047e-07,2.96397490595616e-10,-3.01506961685930e-11,2.22570315950441e-09,-4.26077028191289e-11,2.07336026666512e-09,1.71450079597168e-16,2.94455108664049e-12,8.18246797239659e-17,8.32235684990184e-12,1.05203957047210e-17,5.58541337682320e-12])
+    
+    distortA['r']=np.array([0.000530712822898294,0.999814397089242,-0.000366783811357609,-1.08650906916023e-05,0.000146500108063480,-4.48994291741769e-06,5.66992328511032e-07,-2.10826363791224e-08,7.51050898843171e-07,-7.74459106063614e-09,-5.78712436084704e-11,1.36389366137393e-09,-9.50057131576629e-11,1.49666815592076e-09,-2.16074939825761e-11,2.07124539578673e-12,-1.32787329448528e-13,5.78542850850562e-12,-5.82328830750271e-14,4.04881815088026e-12,2.46095156355282e-15])
+    distortB['r']=np.array([1.00110972320988,-0.000483253233767592,0.999181377618578,9.60316928622784e-05,-9.06574382424422e-06,0.000266147076486786,-6.64513480754565e-09,7.53794491639207e-07,-2.25340150955077e-08,9.53032549996219e-07,2.93844965469408e-10,-6.80521155195455e-11,2.21272371816576e-09,-9.13097728847483e-11,2.06225316601308e-09,9.41096324034730e-14,3.00253083795091e-12,-1.40935245750903e-13,8.27122551593984e-12,-2.41500808188601e-13,5.52074803508046e-12])
+    
+    distortA['i']=np.array([0.000532645905256854,0.999813271238129,-0.000366606839338804,-1.06782246147986e-05,0.000146529811926289,-4.54488847969004e-06,5.68028930846942e-07,-2.11055358580630e-08,7.51187628288423e-07,-7.93885390157909e-09,-5.63104333653064e-11,1.36559362569725e-09,-9.64353839395137e-11,1.50231201562648e-09,-2.04159956020294e-11,1.91408535007684e-12,-7.46369323634455e-14,5.86071470639700e-12,-1.65328163262914e-13,4.07648238374705e-12,3.40880674871652e-14])
+    distortB['i']=np.array([1.00111276400037,-0.000484310769918440,0.999180286636823,9.61240720951134e-05,-8.76526511019577e-06,0.000266192433565878,-6.59462433108999e-09,7.54391611102465e-07,-2.30957980169170e-08,9.53612393886641e-07,2.95558631849358e-10,-7.08307092409029e-11,2.21952307845185e-09,-9.03816603147003e-11,2.06450141393973e-09,-3.77836686311826e-14,3.13358046393416e-12,-5.20417845240954e-14,8.24487152802918e-12,-1.17846469609206e-13,5.35830020655191e-12])
+    
+    x=ra/0.05*10/1000  # convert ra in arcsec to mm
+    
+    y=dec/0.05*10/1000 # convert ra in arcsec to mm
+    
+    dd=np.array([1, x, y, x*x, x*y, y*y, x**3, x**2*y, x*y**2, y**3, x**4, x**3*y, x**2*y**2, x*y**3, y**4,  x**5, x**4*y, x**3*y**2, x**2*y**3 , x*y**4, y**5])
+    
+    xc= np.dot(distortA[ch],dd)  
+    yc= np.dot(distortB[ch],dd)  
+    
+    distortRa =xc*1000/10*0.05-0.00265
+    
+    distortDec=yc*1000/10*0.05-5.0055
+    
+    return distortRa, distortDec
+
+#####################################################################################
+def world_to_pixel(sra,
+                   sdec,
+                   rotsky,
+                   tra,
+                   tdec,
+                   x_center_pixel,
+                   y_center_pixel,
+                   pixelsize=0.05):
+    """_summary_
+
+    Parameters
+    ----------
+    sra : np.array
+         star ra such as np.array([22,23,24]),unit:deg
+    sdec : np.array
+         stat dec such as np.array([10,20,30]),unit:deg
+    rotsky : float
+        rotation angel of the telescope,unit:deg
+    tra : float
+        telescope pointing ra,such as 222.2 deg
+    rdec : float
+        telescope pointing dec,such as 33.3 deg
+    x_center_pixel : float
+        x refer point of MCI,usually the center of the CCD,which start from (0,0)
+    y_center_pixel : float
+        y refer point of MCI,usually the center of the CCD,which start from(0,0)
+    pixelsize : float
+        pixelsize for MCI ccd, default :0.05 arcsec / pixel
+
+    Returns
+    -------
+    pixel_position:list with two np.array
+        such as [array([124.16937605,  99.        ,  99.        ]),
+                array([ 97.52378661,  99.        , 100.50014483])]  
+    """
+    theta_r = rotsky / 180 * np.pi
+    w = WCS(naxis=2)
+    w.wcs.crpix = [x_center_pixel + 1, y_center_pixel + 1]  # pixel from (1, 1)
+    w.wcs.cd = np.array([[
+        -np.cos(-theta_r) * pixelsize / 3600,
+        np.sin(-theta_r) * pixelsize / 3600
+    ],
+                         [
+                             np.sin(-theta_r) * pixelsize / 3600,
+                             np.cos(-theta_r) * pixelsize / 3600
+                         ]])
+    w.wcs.crval = [tra, tdec]
+    w.wcs.ctype = ["RA---TAN", "DEC--TAN"]
+    pixel_position = w.all_world2pix(sra, sdec, 0)
+    return pixel_position
+
+###############################################################################
+
+def cal_pos(center_ra, center_dec, rotTel, rotSky, sRa, sDec):
+# ''' center_ra: telescope pointing in ra direction, unit: degree
+#     center_dec: telescope pointing in dec direction, unit: degree    
+#     rotTel:  telescope totation in degree    
+#     rotSky:  telescope rotation in degree    
+#     sRa:  source ra in arcsec
+#     sDec: source dec in arcsec
+
+    boresight = galsim.CelestialCoord(ra=-center_ra*galsim.degrees, dec=center_dec*galsim.degrees)
+
+    q = rotTel - rotSky    
+    cq, sq = np.cos(q), np.sin(q)
+    jac    = [cq, -sq, sq, cq]
+    affine = galsim.AffineTransform(*jac) 
+    wcs    = galsim.TanWCS(affine, boresight, units=galsim.arcsec)    
+    objCoord  = galsim.CelestialCoord(ra =-sRa*galsim.arcsec, dec=sDec*galsim.arcsec)
+    field_pos = wcs.toImage(objCoord)    
+    return -field_pos.x, field_pos.y  ## return the field position
+
+#####################################################################################
+
+def krebin(a, sample):
+    """ Fast Rebinning with flux conservation
+    New shape must be an integer divisor of the current shape.
+    This algorithm is much faster than rebin_array
+    Parameters
+    ----------
+    a : array_like
+        input array
+    sample :   rebinned sample 2 or 3 or 4 or other int value
+    """
+    # Klaus P's fastrebin from web
+    size=a.shape
+    shape=np.zeros(2,dtype=int)
+    shape[0]=int(size[0]/sample)
+    shape[1]=int(size[1]/sample)
+    
+    sh = shape[0], a.shape[0] // shape[0], shape[1], a.shape[1] // shape[1]
+    
+    return a.reshape(sh).sum(-1).sum(1)
+
+#####################################################################
+#########################################################
+
+def centroid(data):
+    '''
+    calculate the centroid of the input two-dimentional image data
+
+    Parameters
+    ----------
+    data : input image.
+
+    Returns
+    -------
+    cx: the centroid column number,  in horizontal direction definet in python image show 
+    cy: the centroid row number  ,   in vertical direction 
+
+    '''
+    ###
+    from scipy import ndimage
+    cy,cx=ndimage.center_of_mass(data)
+    return cx, cy
+
+############################################################################### 
+
+
+
+###############################################################################
+
+
+#########################################################################
+
+'''
+PSF interpolation for MCI_sim
+'''
+
+import scipy.spatial as spatial
+
+###find neighbors-KDtree  ###
+def findNeighbors(tx, ty, px, py, dn=5):
+    """
+    find nearest neighbors by 2D-KDTree
+
+    Parameters:
+        tx, ty (float, float): a given position
+        px, py (numpy.array, numpy.array): position data for tree
+        dr (float-optional): distance in arcsec
+        dn (int-optional): nearest-N
+        OnlyDistance (bool-optional): only use distance to find neighbors. Default: True
+
+    Returns:
+        indexq (numpy.array): index
+    """
+    datax = px
+    datay = py
+    tree = spatial.KDTree(list(zip(datax.ravel(), datay.ravel())))
+    indexq=[]   
+    dist, indexq = tree.query([tx, ty], dn)
+    return indexq
+
+###############################################################################
+###PSF-IDW###
+def psfMaker_IDW(px, py, PSFMat, cen_col, cen_row, IDWindex=5, OnlyNeighbors=True):
+    """
+    psf interpolation by IDW
+
+    Parameters:
+        px, py (float, float): position of the target
+        PSFMat (numpy.array): PSF matrix array at given position and wavelength
+        cen_col, cen_row (numpy.array, numpy.array): potions of the psf centers in arcsec
+        IDWindex (int-optional): the power index of IDW
+        OnlyNeighbors (bool-optional): only neighbors are used for psf interpolation
+
+    Returns:
+        psfMaker (numpy.array)
+    """
+
+    minimum_psf_weight = 1e-8
+    ref_col = px  ### target position in x 
+    ref_row = py  ### target position in y
+
+    ngy, ngx = PSFMat[:, :, 0, 0].shape  ### PSF data size
+    npsf = PSFMat[:, :, :,:].shape[3]  ### total psf number in one wavelength
+    psfWeight = np.zeros([npsf])
+
+    if OnlyNeighbors == True:
+        
+        neigh = findNeighbors(px, py, cen_col, cen_row, dn=5)
+        # if hoc is not None:
+        #     neigh = findNeighbors_hoclist(cen_col, cen_row, tx=px,ty=py, dn=4, hoc=hoc, hoclist=hoclist)
+
+        neighFlag = np.zeros(npsf)
+        neighFlag[neigh] = 1
+
+    for ipsf in range(npsf):
+        if OnlyNeighbors == True:
+            if neighFlag[ipsf] != 1:
+                continue
+
+        dist = np.sqrt((ref_col - cen_col[ipsf])**2 + (ref_row - cen_row[ipsf])**2)
+        if IDWindex == 1:
+            psfWeight[ipsf] = dist
+        if IDWindex == 2:
+            psfWeight[ipsf] = dist**2
+        if IDWindex == 3:
+            psfWeight[ipsf] = dist**3
+        if IDWindex == 4:
+            psfWeight[ipsf] = dist**4
+            
+        psfWeight[ipsf] = max(psfWeight[ipsf], minimum_psf_weight)
+        psfWeight[ipsf] = 1./psfWeight[ipsf]
+        
+    psfWeight /= np.sum(psfWeight)
+
+    psfMaker  = np.zeros([ngy, ngx, 7 ], dtype=np.float32)
+    
+    
+    for waven in range(7):
+    
+        for ipsf in range(npsf):
+            if OnlyNeighbors == True:
+                if neighFlag[ipsf] != 1:
+                    continue
+    
+            iPSFMat = PSFMat[:, :, waven,ipsf].copy()
+            ipsfWeight = psfWeight[ipsf]
+    
+            psfMaker[:, :, waven] += iPSFMat * ipsfWeight
+            
+        psfMaker[:, :, waven] /= np.nansum(psfMaker[:, :, waven])
+
+    return psfMaker
+
+###############################################################################
+###############################################################################    
+###############################################################################
+
+class MCIsimulator():
+    """
+    CSST MCI Simulator
+
+    The image that is being build is in::
+
+        self.image
+        self.image_g    g channel
+        self.image_r    r channel
+        self.image_i    i channel
+        
+
+    :param opts: OptionParser instance
+    :type opts: OptionParser instance
+    """
+
+    def __init__(self, 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
+
+
+    ###########################################################################
+    
+
+        self.information = MCIinstrumentModel.MCIinformation()
+        #update settings with defaults
+        self.information.update(dict(quadrant=int(0),
+                                     ccdx=int(0),
+                                     ccdy=int(0),                                
+                                     ccdxgap=1.643,
+                                     ccdygap=8.116,
+                                     xsize=2000,
+                                     ysize=2000,
+                                     fullwellcapacity=90000,
+                                     dark=0.001,
+                                     exptime=300.0,
+                                     readouttime=4.,
+                                     rdose=8.0e9,                                  
+                                     ghostCutoff=22.0,
+                                     ghostRatio=5.e-5,
+                                     coveringfraction=1.0,
+                                     pixel_size=0.05))
+        
+        
+        self.configure(1)    #print the configfile name and path; 
+        
+        self.information.update(dict(
+                                                                     
+        cosmicraylengths =self.information['indata_path']+'/data/cdf_cr_length.dat',
+        cosmicraydistance=self.information['indata_path']+'/data/cdf_cr_total.dat',                                
+        parallelTrapfile =self.information['indata_path']+'/data/cdm_euclid_parallel.dat',
+        serialTrapfile   =self.information['indata_path']+'/data/cdm_euclid_serial.dat',
+        cosmeticsfile_g  =self.information['indata_path']+'/data/Cosmetics_g.txt',
+        cosmeticsfile_r  =self.information['indata_path']+'/data/Cosmetics_r.txt' ,
+        cosmeticsfile_i  =self.information['indata_path']+'/data/Cosmetics_i.txt'     ))
+        
+###############################################################################
+    def readConfigs(self,simnumber):
+        """
+        Reads the config file information using configParser and sets up a logger.
+        """
+        self.config = ConfigParser.RawConfigParser()
+        
+      
+        
+        #self.config.readfp(open(self.configfile))
+        self.config.read_file(open(self.configfile))
+        
+
+       
+################################################################################################3
+
+    def processConfigs(self):
+        """
+        Processes configuration information and save the information to a dictionary self.information.
+        For explanation of each field, see /data/test.config. Note that if an input field does not exist,
+        then the values are taken from the default instrument model as described in
+        support.MCIinstrumentModel.VISinformation(). Any of the defaults can be overwritten by providing
+        a config file with a correct field name.
+        """
+        #parse options and update the information dictionary
+        options = self.config.options(self.section)
+        settings = {}
+        for option in options:
+            try:
+                settings[option] = self.config.getint(self.section, option)
+            except ValueError:
+                try:
+                    settings[option] = self.config.getfloat(self.section, option)
+                except ValueError:
+                    settings[option] = self.config.get(self.section, option)
+
+        self.information.update(settings)
+
+
+
+        #ghost ratio can be in engineering format, so getfloat does not capture it...
+        try:
+            self.information['ghostRatio'] = float(self.config.get(self.section, 'ghostRatio'))
+        except:
+            pass
+        #########
+
+        #booleans to control the flow
+        
+        self.cosmicRays      = self.config.getboolean(self.section, 'cosmicRays')
+        self.darknoise       = self.config.getboolean(self.section, 'darknoise')
+        self.cosmetics       = self.config.getboolean(self.section, 'cosmetics')
+        self.radiationDamage = self.config.getboolean(self.section, 'radiationDamage')      
+        self.bleeding        = self.config.getboolean(self.section, 'bleeding')
+        self.overscans       = self.config.getboolean(self.section, 'overscans')
+        self.nonlinearity    = self.config.getboolean(self.section, 'nonlinearity')      
+        self.readoutNoise    = self.config.getboolean(self.section, 'readoutNoise')            
+        self.skyback      = self.config.getboolean(self.section, 'skyback')
+        self.TianceEffect = self.config.getboolean(self.section, 'TianceEffect')
+        self.intscale     = self.config.getboolean(self.section, 'intscale')          
+        self.ghosts       = self.config.getboolean(self.section, 'ghosts')            
+        self.shutterEffect        =self.config.getboolean(self.section, 'shutterEffect') 
+        self.flatfieldM           =self.config.getboolean(self.section, 'flatfieldm') 
+        self.PRNUeffect           =self.config.getboolean(self.section, 'PRNUeffect')        
+        self.appFatt              =self.config.getboolean(self.section, 'appFatt')        
+        self.sky_shift_rot        =self.config.getboolean(self.section, 'sky_shift_rot')        
+        self.distortion           =self.config.getboolean(self.section, 'distortion')
+        self.sim_star             =self.config.getboolean(self.section, 'sim_star')
+        self.sim_galaxy           =self.config.getboolean(self.section, 'sim_galaxy')
+        
+        
+        ###############################################3################################# 
+        self.booleans = dict(cosmicRays     =self.cosmicRays,
+                             darknoise      =self.darknoise,
+                             cosmetics      =self.cosmetics,
+                             radiationDamage=self.radiationDamage,
+                             bleeding       =self.bleeding,
+                             overscans      =self.overscans,
+                             nonlinearity   =self.nonlinearity,                             
+                             readoutNoise   =self.readoutNoise,                             
+                             skyback        =self.skyback,
+                             TianceEffect   =self.TianceEffect,
+                             intscale       =self.intscale,
+                             ghosts         =self.ghosts ,
+                             shutterEffect      =self.shutterEffect,                             
+                             flatfieldM         =self.flatfieldM,
+                             PRNUeffect         =self.PRNUeffect,
+                             appFatt            =self.appFatt ,
+                             sky_shift_rot      =self.sky_shift_rot,
+                             distortion         =self.distortion,
+                             sim_star           =self.sim_star,
+                             sim_galaxy         =self.sim_galaxy)        
+        ##################################################################### 
+        
+        
+        now=datetime.now()
+        #data_time=now.strftime("%Y-%m-%d-%H-%M-%S")
+        result_day=now.strftime("%Y-%m-%d")
+             
+        self.result_path=self.information['result_path']+'/'+result_day ### CSST1      
+
+        
+        if os.path.isdir(self.result_path)==False:
+            os.mkdir(self.result_path)
+            os.mkdir(self.result_path+'/cali_Data')
+            os.mkdir(self.result_path+'/log_Data')           
+            os.mkdir(self.result_path+'/sky_Data')
+        #############################################################           
+            
+               
+        data_time=now.strftime("%Y-%m-%d-%H-%M-%S")
+        
+
+        self.log = lg.setUpLogger(self.result_path+'/log_Data/MCIsim_'+'_'+data_time+'.log')
+        
+        
+        self.log.info('-------STARTING A NEW SIMULATION------------')
+        
+        
+        self.log.info(self.information)
+        
+        ##############################################
+                #load instrument model, these values are also stored in the FITS header
+        
+        self.information['G_filters']=["F275W",  "F280N","NUV",   "WU",   "CBU",   "F343N",  "u",      "F373N", "F395N",  "F336W"]
+        self.information['R_filters']=["F487N", "F502N", "CBV",   "r",    "F656N", "F658N",  "F467M",  "F555W", "F606W",  "F673N"]        
+        self.information['I_filters']=["z",      "y",    "F815N", "CBI",  "F925N", "F960M",  "F968N",  "F845M" ,"F850LP" ,"F814W"]
+        
+        #### load telescope efficiency data
+        self.tel_eff=np.load(self.information['indata_path']+'/tel_eff/tel_eff.npy',allow_pickle=True).item()
+        
+        #### load MCI filter data
+        self.filterP=np.load(self.information['indata_path']+'/MCI_filters/mci_filterPWTC.npy',allow_pickle=True).item()    
+        
+                
+        #####################################################################
+
+        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))
+            
+        return
+            
+ #########################################################################################           
+    def make_c_coor(self, bs, nc):
+        '''
+            Draw the mesh grids for a bs*bs box with nc*nc pixels
+        '''
+
+        ds=bs/nc
+        xx01 = np.linspace(-bs/2.0,bs/2.0-ds,nc)+0.5*ds
+        xx02 = np.linspace(-bs/2.0,bs/2.0-ds,nc)+0.5*ds
+        xg2,xg1 = np.meshgrid(xx01,xx02)
+        return xg1,xg2
+##########################################################################################
+
+    def _createEmpty(self):
+        """
+        Creates and empty array of a given x and y size full of zeros.
+        add g r i channel images;
+        
+        Creates lensing parameters;  
+    
+        """        
+        self.image_g=np.zeros((self.information['ysize'], self.information['xsize']), dtype=np.float64)
+        self.image_r=np.zeros((self.information['ysize'], self.information['xsize']), dtype=np.float64)
+        self.image_i=np.zeros((self.information['ysize'], self.information['xsize']), dtype=np.float64)
+        return
+        
+#########################################################################################################################
+
+#########################################################################################################
+
+
+    def _loadGhostModel(self):
+        """
+        Reads in a ghost model from a FITS file and stores the data to self.ghostModel.
+
+        Currently assumes that the ghost model has already been properly scaled and that the pixel
+        scale of the input data corresponds to the nominal VIS pixel scale. Futhermore, assumes that the
+        distance to the ghost from y=0 is appropriate (given current knowledge, about 750 VIS pixels).
+        """
+
+        self.ghostOffsetX=self.information['ghostoffsetx']
+        self.ghostOffsetY=self.information['ghostoffsety']
+        
+        self.ghostMax = self.information['ghostratio']
+        self.log.info('Maximum in the ghost model %e' % self.ghostMax)
+        return
+   ###############################################
+
+    def readCosmicRayInformation(self):
+        """
+        Reads in the cosmic ray track information from two input files.
+
+        Stores the information to a dictionary called cr.
+        """
+        self.log.info('Reading in cosmic ray information from %s and %s' % (self.information['cosmicraylengths'],
+                                                                            self.information['cosmicraydistance']))
+
+        crLengths = np.loadtxt(self.information['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])
+        
+        return
+
+###############################################################################
+
+###############################################################################
+    def configure(self,simnumber):
+        """
+        Configures the simulator with input information and creates and empty array to which the final image will
+        be build on.
+        """
+        self.readConfigs(simnumber)
+        
+        self.processConfigs()
+        
+        self._createEmpty()
+        
+        self.log.info('Read in the configuration files and created an empty array')
+        
+        return
+ 
+################################################################################
+    def load_filter_PSF(self):
+        ##### load filter PSF in three channels ,G , R and I       
+        ###
+        # load filter name in three channels from information
+        self.filter_g=self.information['filter_g']
+        self.filter_r=self.information['filter_r']
+        self.filter_i=self.information['filter_i']
+        
+        self.filter_psf=dict()
+        
+        filtername=self.filter_g
+        self.filter_psf['g']=np.load(self.information['indata_path']+'/PSF/'+filtername+'_PSF.npy', allow_pickle=True).item()
+        
+        filtername=self.filter_r
+        self.filter_psf['r']=np.load(self.information['indata_path']+'/PSF/'+filtername+'_PSF.npy', allow_pickle=True).item()
+        
+        filtername=self.filter_i
+        self.filter_psf['i']=np.load(self.information['indata_path']+'/PSF/'+filtername+'_PSF.npy', allow_pickle=True).item()
+       
+        return  
+    
+############################################################################################### 
+
+    def get_PSF(self, px, py, ch):
+        """
+        Get the PSF at a given image position
+
+        Parameters:
+        px,py        :  target position
+        filternmae   :  the selected filter name
+        ch           :  the MCI channle , g or r or i 
+        Returns:
+            PSF:     PSF array at given 7 wavelength and at target position
+        """
+        PSFMat = self.filter_psf[ch]['psf_mat'] 
+        cen_col= self.filter_psf[ch]['psf_field_X'] 
+        cen_row= self.filter_psf[ch]['psf_field_Y'] 
+        imPSF = psfMaker_IDW(px, py, PSFMat, cen_col, cen_row, IDWindex=2)
+                
+        return imPSF
+
+###############################################################################
+
+############################################################################### 
+   
+    def cal_SED_photons(self, flux):        
+        ## flux_arr,                input flux array with unit of  1e-17 erg/s/A/cm^2  
+
+        wave = np.linspace(2500,10000,8501)     # set the wavelenth for MCI from 250nm to 1100nm
+        wave=wave/10 
+            
+        # calcutle photons from original spec cube data,
+        # data: input data, two-dimentional, unit :  1e-17 erg/s/A/cm^2
+      
+        planckh= 6.62620*10**-27                #  % erg s;   
+        cc=2.99792458*10**17                    #  in nm/s       
+        telarea=3.1415926*100*100               # in cm^2，望远镜聚光面积�? 
+        fluxlam=1e-17*(flux)                    # convert original unit to unit of erg/s/A/cm^2
+                                                # wave in nm   ;;      
+        ephoton=planckh*cc/wave                 # single photon energy in erg/photon;  cc与lambda单位需要一致；
+        
+        Nphoton =fluxlam/ephoton*telarea*self.information['exptime']        # in unit of  photons       
+                    
+        #### load telescope efficiency data
+        teleff=dict()
+        teleff['g']=self.tel_eff['G'] ### fisrt colunmm  is wavelength in nm
+        teleff['r']=self.tel_eff['R'] ### second colunmm is  efficiency
+        teleff['i']=self.tel_eff['I']
+        #####   
+        
+        ft=dict()    
+        ft['g']=self.filter_g
+        ft['r']=self.filter_r
+        ft['i']=self.filter_i
+        
+        
+        # load filter name in three channels
+        ft_wave=dict()
+        ft_wave['g']=self.filterP[ft['g']]['wave_nm']
+        ft_wave['r']=self.filterP[ft['r']]['wave_nm']
+        ft_wave['i']=self.filterP[ft['i']]['wave_nm']
+        
+        ft_eff=dict()
+        ft_eff['g']=self.filterP[ft['g']]['throughout']
+        ft_eff['r']=self.filterP[ft['r']]['throughout']
+        ft_eff['i']=self.filterP[ft['i']]['throughout']
+        
+        chlist=['g','r','i']
+        
+        Sumphoton=dict()
+        PSF_eff_weight=dict()
+        
+        for k in range(3):
+            ch=chlist[k]
+            tel_wavearr=teleff[ch][:,0]
+            tel_effarr =teleff[ch][:,1]        
+            tel_effnew=np.interp(wave, tel_wavearr, tel_effarr)
+            
+            ft_wavearr=ft_wave[ch]
+            ft_effarr =ft_eff[ch]
+            ft_effnew =np.interp(wave, ft_wavearr, ft_effarr)           
+                       
+            
+            #### SED flux*filter_flux*tel_eff
+            eff_arr=Nphoton*ft_effnew*tel_effnew            
+            
+            Sumphoton[ch]=np.sum(eff_arr)  
+            
+            iwave=self.filter_psf[ch]['psf_iwave'][:]  ### seven wavelength for the seleced filter
+            
+            psf_coeff =np.interp(iwave, wave, eff_arr)
+            
+            if psf_coeff.sum()==0:                
+                psf_coeff=1.0/7*np.ones(7)                
+            else:
+                psf_coeff=psf_coeff/psf_coeff.sum()              
+              
+            PSF_eff_weight[ch]=psf_coeff          
+            
+        return Sumphoton, PSF_eff_weight      
+
+    ##################################################
+############################################################################### 
+   
+    def cal_sky_noise(self):
+        #######       add earthshine       #######
+        #######       add earthshine       #######
+        #self.earthshine_wave        # A
+        #self.earthshine_flux        # erg/s/cm^2/A/arcsec^2
+        
+        wave = np.linspace(2500,11000, 8501)     # set the wavelenth for MCI from 250nm to 1100nm
+        
+        wavearr =self.earthshine_wave    # A
+        fluxarr =self.earthshine_flux    # erg/s/cm^2/A/arcsec^2  
+        earthshinefluxarr=np.interp(wave, wavearr, fluxarr)  # erg/s/cm^2/A/arcsec^2  
+        
+        wavearr =self.zodiacal_wave     # A
+        fluxarr =self.zodiacal_flux   # erg/s/cm^2/A/arcsec^2  
+        zodiacalfluxarr=np.interp(wave, wavearr, fluxarr)  # erg/s/cm^2/A/arcsec^2         
+
+        skynoise_flux=earthshinefluxarr+zodiacalfluxarr   # erg/s/cm^2/A/arcsec^2 
+                                 
+        wave = np.linspace(2500,10000,8501)     # set the wavelenth for MCI from 250nm to 1100nm
+        wave=wave/10 
+            
+        # data: input data, two-dimentional, unit :  1e-17 erg/s/A/cm^2      
+        planckh= 6.62620*10**-27                #  % erg s;   
+        cc=2.99792458*10**17                    #  in nm/s       
+        telarea=3.1415926*100*100               # in cm^2，   
+        ephoton=planckh*cc/wave                 # in erg/photon,  
+        Nphoton_flux =skynoise_flux/ephoton*telarea*0.05*0.05*self.information['exptime']    # each pixel will have the photons in unit of  photons/s/A        
+          
+        #### load telescope efficiency data
+        teleff=dict()
+        teleff['g']=self.tel_eff['G'] ### fisrt colunmm  is wavelength in nm
+        teleff['r']=self.tel_eff['R'] ### second colunmm is  efficiency
+        teleff['i']=self.tel_eff['I']
+        #####   
+        
+        ft=dict()    
+        ft['g']=self.filter_g
+        ft['r']=self.filter_r
+        ft['i']=self.filter_i
+                
+        # load filter name in three channels
+        ft_wave=dict()
+        ft_wave['g']=self.filterP[ft['g']]['wave_nm']
+        ft_wave['r']=self.filterP[ft['r']]['wave_nm']
+        ft_wave['i']=self.filterP[ft['i']]['wave_nm']
+        
+        ft_eff=dict()
+        ft_eff['g']=self.filterP[ft['g']]['throughout']
+        ft_eff['r']=self.filterP[ft['r']]['throughout']
+        ft_eff['i']=self.filterP[ft['i']]['throughout']
+        
+        chlist=['g','r','i']
+        
+        self.Sky_Noise=dict()
+
+        
+        for k in range(3):
+            ch=chlist[k]
+            tel_wavearr=teleff[ch][:,0]
+            tel_effarr =teleff[ch][:,1]        
+            tel_effnew=np.interp(wave, tel_wavearr, tel_effarr)
+            
+            ft_wavearr=ft_wave[ch]
+            ft_effarr =ft_eff[ch]
+            ft_effnew =np.interp(wave, ft_wavearr, ft_effarr)            
+                        
+            self.Sky_Noise[ch]=np.sum(Nphoton_flux*ft_effnew*tel_effnew)  
+            
+        return self.Sky_Noise  
+
+    ##################################################
+    ##############################################################################
+
+    def cal_StarSED_img(self): 
+        #####################
+        
+        pixelscale=self.information['pixel_size']    ## arcsec, pixel scale size
+        ghostOffsetX =self.information['ghostoffsetx']
+        ghostOffsetY =self.information['ghostoffsety']
+        ghostMx      =self.information['ghostratio'] 
+
+        rotTelPos=self.information['rotTelPos']
+        rotSkyPos=self.information['rotSkyPos']
+        
+        theta = rotTelPos - rotSkyPos          
+        
+        ############ load star data catlog #####################    
+     
+        df0=pd.read_csv(self.information['indata_path']+'/star_input/GaiaSource_675688-675713.csv')   
+        
+        ####
+        df=df0.fillna(1000)
+        df2=df[ (df['phot_rp_mean_mag']>16)&(df['phot_bp_mean_mag']>16) & (df['phot_bp_mean_mag']<31) & (df['pmra']<100) & (df['phot_rp_mean_mag']<31)]
+        
+        df2.index = range(len(df2))      
+        ## limit the filed of view to 600*600 arcsec^2 square zone;
+        df3=df2[ (abs(df2['ra']-df2['ra'].mean())<300/3600.0) & (abs(df2['dec']-df2['dec'].mean())<300/3600.0) ]
+        
+        df3.index = range(len(df3))
+        
+        self.star=df3
+        del df0,df,df2,df3
+        
+        self.information['ra_obj']  = self.star['ra'].mean()
+        self.information['dec_obj'] = self.star['dec'].mean()
+        
+        center_ra =self.information['T_disRa']  +self.information['ra_obj']
+        center_dec=self.information['T_disDec'] +self.information['dec_obj']      
+        
+        self.information['ra_pnt0'] =center_ra
+        self.information['dec_pnt0']=center_dec  
+        
+        channel=['g','r','i']  
+        #######################################################################
+     
+        nsrcs=len(self.star) 
+        
+        ################################################################## 
+        obj=dict()   
+        obj['g']=np.zeros(3)            
+        obj['r']=np.zeros(3)
+        obj['i']=np.zeros(3)
+        
+        fullimg=dict()   
+        star_input=dict()
+        star_input['g']=np.zeros((nsrcs,3))
+        star_input['r']=np.zeros((nsrcs,3))
+        star_input['i']=np.zeros((nsrcs,3))
+        
+        star_output=dict()
+        star_output['g']=np.zeros((nsrcs,3))
+        star_output['r']=np.zeros((nsrcs,3))
+        star_output['i']=np.zeros((nsrcs,3))
+        
+
+        final_image=dict()      
+        
+        final_image['g'] = galsim.ImageF(int(self.information['xsize']), int(self.information['ysize']))
+        final_image['r'] = galsim.ImageF(int(self.information['xsize']), int(self.information['ysize']))
+        final_image['i'] = galsim.ImageF(int(self.information['xsize']), int(self.information['ysize']))
+        
+        final_image['g'].setOrigin(0,0)
+        final_image['r'].setOrigin(0,0)
+        final_image['i'].setOrigin(0,0)   
+        #######################################################################
+        nlayccd = 0    
+          
+        #### calculate sky noise #####
+        self.earthshine(self.earthshine_theta)
+        
+        self.zodiacal(self.information['ra_obj'], self.information['dec_obj'], self.dt.strftime("%Y-%m-%d"))
+        
+        self.cal_sky_noise()  
+
+        #################################   
+
+        self.information['target']  =deg2HMS(self.information['ra_obj'], self.information['dec_obj'])  
+        
+        self.information['CRVAL1']=self.information['ra_pnt0']
+        self.information['CRVAL2']=self.information['dec_pnt0']
+        
+        self.information['CRPIX1']=self.information['xsize']/2-0.5 #### 
+        self.information['CRPIX2']=self.information['ysize']/2-0.5 #### 
+        
+        self.information['CD1_1']=-np.cos(-theta)*self.information['pixel_size']/3600.0 #### 
+        self.information['CD1_2']= np.sin(-theta)*self.information['pixel_size']/3600.0 #### 
+        
+        self.information['CD2_1']= np.sin(-theta)*self.information['pixel_size']/3600.0 #### 
+        self.information['CD2_2']= np.cos(-theta)*self.information['pixel_size']/3600.0 ####  
+        
+        ####################################################################### 
+        
+        if self.TianceEffect:          
+            ra_list = self.star['ra'].tolist()
+            dec_list = self.star['dec'].tolist()
+            pmra_list = self.star['pmra'].tolist()
+            pmdec_list = self.star['pmdec'].tolist()
+            parallax_list = self.star['parallax'].tolist()
+            rv_list = [0.0 for i in range(len(ra_list))]            
+            
+            ################################################  
+            
+            newRa, newDec = shao.onOrbitObsPosition(ra_list, dec_list, pmra_list, \
+                                              pmdec_list, rv_list, parallax_list, len(ra_list), \
+                                              self.information['pos_x'], self.information['pos_y'], self.information['pos_z'], self.information['velocity_x'],self.information['velocity_y'],self.information['velocity_z'], "J2000", self.TianCe_day, self.TianCe_exp_start)
+        else:
+            newRa  =self.star['ra']
+            newDec =self.star['dec']
+        ######################################################
+
+        ####################   generate star image  ##########
+        for j in tqdm(range(nsrcs)): ###   nsrcs length       
+            
+            
+            starRa  = 3600.0*( newRa[j]   )  # ra of star,  arcsecond
+            starDec = 3600.0*( newDec[j]  )  # dec of star, arcsecond
+              
+            ###################################################################
+            fsx,fsy=cal_pos(center_ra,center_dec, rotTelPos, rotSkyPos, starRa, starDec)   
+            row=  fsy/pixelscale+self.information['ysize']/2-0.5 ### row number on CCD image
+            col=  fsx/pixelscale+self.information['xsize']/2-0.5 ### col number on CCD image
+            
+            
+            if row>=self.information['ysize']+100 or row<=-100 or col>=self.information['xsize']+100 or col <=-100:
+                continue
+                 
+            nlayccd=nlayccd+1
+            
+            ################################################################
+            for i in range(3):
+                ###                        
+                ch=channel[i]  
+                star_input[ch][int(nlayccd)-1, 0]   =fsx       #ra
+                star_input[ch][int(nlayccd)-1, 1]   =fsy       #dec                 
+ 
+            #############  do field distottion  ##########
+         
+            if self.distortion:                   
+                
+                for i in range(3):
+                    
+                    ch=channel[i]                     
+                    fpx,fpy=distortField(fsx, fsy, ch)                        
+                    obj[ch][0]=fpx
+                    obj[ch][1]=fpy                                      
+                    
+                    star_output[ch][int(nlayccd)-1, 0] =fpx       #ra
+                    star_output[ch][int(nlayccd)-1, 1] =fpy       #dec                            
+             
+            else:
+                
+                fpx=fsx
+                fpy=fsy                   
+                
+                for i in range(3):  
+                    
+                    ch=channel[i] 
+                   
+                    obj[ch][0]=fpx
+                    obj[ch][1]=fpy                    
+                    
+                    star_output[ch][int(nlayccd)-1, 0]   =fpx       #ra
+                    star_output[ch][int(nlayccd)-1, 1]   =fpy       #dec                    
+                
+            #######################################################################
+         
+            ####### use  SED_code to generate star SED   #######
+            # SED of j-th star        
+
+            if j<nsrcs-8:
+                bp  =  self.star['phot_bp_mean_mag'][j]
+                rp  =  self.star['phot_rp_mean_mag'][j] 
+            else:
+
+                bp  =  binary_star[j-nsrcs, 2]
+                rp  =  bp
+                
+            wave = np.linspace(2500, 11000, 8501)
+            
+            # Loading stellar SED template
+            t=sed.Star_Temp(self.information['indata_path'])
+            # Calculating the magnitude (Bp, Rp) of each template
+            t.toMag()
+            # Calculating flux
+            star_flux = sed.Model_Stellar_SED(wave, bp, rp, t, self.information['indata_path'])  # unit of 10-17 erg/s/A/cm2             
+            #########################################
+            
+            ##cal_SED_photons(self, flux_arr):
+            intscales,PSF_eff_weight=self.cal_SED_photons(star_flux) 
+            
+            if j<nsrcs-8:
+                self.log.info('star number=%i, original Ra(in deg)=%f, original Dec(in deg)=%f, bp=%f,rp=%f, TotPhotons_G=%f, TotPhotons_R=%f,TotPhotons_I=%f' % (nlayccd, newRa[j], newDec[j],self.star['phot_bp_mean_mag'][j], self.star['phot_rp_mean_mag'][j],intscales['g'],intscales['r'],intscales['i']))
+   
+            
+            gx=dict()
+            gy=dict()
+            for i in range(3):                
+                ch=channel[i]                   
+                gx[ch]= obj[ch][1]/pixelscale+self.information['ysize']/2-0.5 ### row number on CCD image
+                gy[ch]= obj[ch][0]/pixelscale+self.information['xsize']/2-0.5 ### col number on CCD image
+                
+                self.log.info('Channel in =%s, PosX(in pixel)=%f, PosY(in pixel)=%f' % (channel[i], gx[ch], gy[ch]))
+                
+            ######################################################################
+            psf=dict() 
+            
+            for i in range(3):
+
+                ch=channel[i]  
+                #######
+                psfmat=self.get_PSF(fsx, fsy, ch)                    
+                                
+                for iwave in range(7):
+                    if iwave==0:
+                        temp=PSF_eff_weight[ch][iwave]*psfmat[:,:,iwave]
+                    else:
+                        temp=temp+PSF_eff_weight[ch][iwave]*psfmat[:,:,iwave]
+                
+                temp=temp/temp.sum()                        
+                             
+                ####rotate the PSF data
+                if abs(theta.deg>0):
+                    psf[ch]=ndimage.rotate(temp, theta.deg, order=1, reshape=True) # here we choose reshape=False, the rotated image will
+                else:
+                    psf[ch]=temp                    
+
+                conv = psf[ch]                         
+                conv=conv/conv.sum()
+                
+                conv=conv*intscales[ch]
+                #################################################                                                
+                                                        
+                stamp_img = galsim.Image(conv, copy=True)                        
+                stamp_img.setOrigin(0,0)                            
+                stamp_img.scale=0.025                             
+            
+                #################################################   
+                    
+                photons=galsim.PhotonArray.makeFromImage(stamp_img, max_flux=max(1.0, stamp_img.array.max()/10000.0))
+                
+                cx0,cy0=centroid(stamp_img.array)
+                
+                if self.appFatt : 
+                    
+                ### apply treering and bright fatter and diffusion;                        
+                    SimpleTreeRing = galsim.SiliconSensor().simple_treerings(amplitude=self.information['treering'])                  
+                    
+                    cx, cy = centroid(conv)
+                        
+                    pos = galsim.PositionD(x=int(cx), y=int(cy)) ### set subimge center as the pos 
+       
+                    siliconsensor = galsim.SiliconSensor(strength=self.information['fatter'], diffusion_factor=self.information['difusion'], treering_func=SimpleTreeRing, treering_center=pos)
+                                           
+                    image = galsim.Image(conv.shape[0], conv.shape[1])                            
+                    image.scale = 0.025 
+                    image.setOrigin(0,0)
+                    
+                    photons.x=photons.x/image.scale
+                    photons.y=photons.y/image.scale
+                    
+                    siliconsensor.accumulate(photons, image)   
+                    
+                    photons=galsim.PhotonArray.makeFromImage(image, max_flux=1.0)
+                    
+                    cx0,cy0=centroid(image.array)
+                        
+                ############################################
+                ############################################
+                photons.x=photons.x-cx0*0.025
+                photons.y=photons.y-cy0*0.025
+                
+                ##########################################################
+                #### add ghost image ####
+                if self.ghosts:
+                    ghostphotons=galsim.PhotonArray(len(photons.x),x=photons.x, y=photons.y, flux=photons.flux)
+                    ghostphotons.flux=ghostMx*ghostphotons.flux
+                    ghostphotons.x=photons.x/0.05+ghostOffsetY
+                    ghostphotons.y=photons.y/0.05+ghostOffsetX
+                    
+                    if ghostphotons.flux.max()>0.1/1500:
+                        ghostphotons.addTo(final_image[ch]) 
+                                                                    
+                #########################################################                
+                photons.x=photons.x/0.05+gy[ch]
+                photons.y=photons.y/0.05+gx[ch]
+                
+                photons.addTo(final_image[ch]) 
+            
+
+                ###############################################################
+                ###############################################################
+        #print('total star nummber on CCD is: ',nlayccd )
+    
+        #################################################################  
+
+        ####################################################################### 
+                
+        fullimg['g']=final_image['g'].array
+        fullimg['r']=final_image['r'].array
+        fullimg['i']=final_image['i'].array
+    
+    
+        return fullimg
+    
+################################################################################
+################################################################################
+#################################################################################
+    ########################################################################
+    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)
+        #xx, yy = np.meshgrid(beta_angle, lamda_angle,indexing='ij', sparse=True)
+        
+        f = interpolate.interp2d(xx, yy, zodi, kind='linear')
+        #f = interpolate.RegularGridInterpolator((xx, yy), zodi, method='linear')
+        
+        zodi_obj = f(beta, lamda)       # 
+    
+        # 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      # 
+        zodi_norm = 252                 # 
+    
+        spec = spec0 * (zodi_obj / zodi_norm) * 1e-8  # 
+    
+        # 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 cal_Lensing_galaxy_img(self): 
+        #####################
+         
+        
+        self.dsx_arc = self.information['pixel_size'] # arcsec, pixel scale size;
+        self.bsz_arc = max(self.information['ysize'] ,self.information['xsize'])*self.dsx_arc  # arcsec, field size;
+        self.bsz_deg = self.bsz_arc/3600.
+        
+        #############################################################
+        ### galaxy catlog information
+        ra =   74.0563449685417
+        dec = -34.155241212932374
+        # halo_id = 229600100382
+        
+        self.information['star_ra']  = 249.70093286177536
+        self.information['star_dec'] =-44.60485406119162
+        
+        self.information['gal_ra'] = ra
+        self.information['gal_dec']= dec
+        
+        self.information['ra_obj']  =self.information['star_ra']
+        self.information['dec_obj'] =self.information['star_dec']
+        
+        #####################################################################
+        self.earthshine(self.earthshine_theta)
+        
+        self.zodiacal(self.information['star_ra'], self.information['star_dec'], self.dt.strftime("%Y-%m-%d"))   
+       
+        self.cal_sky_noise()  
+               
+        ############ load galaxy data with SED  ############################        
+    
+        losscale =self.information['pixel_size']
+        rotTelPos=self.information['rotTelPos']
+        rotSkyPos=self.information['rotSkyPos']
+        
+        theta = rotTelPos - rotSkyPos 
+        
+        center_ra  =self.information['T_disRa']  +self.information['star_ra']
+        center_dec =self.information['T_disDec'] +self.information['star_dec']     
+        
+        self.information['ra_pnt0'] =center_ra
+        self.information['dec_pnt0']=center_dec
+        
+        self.information['target']  =deg2HMS(self.information['star_ra'], self.information['star_dec']) 
+        
+        self.information['CRVAL1']=self.information['ra_pnt0']
+        self.information['CRVAL2']=self.information['dec_pnt0']
+        
+        self.information['CRPIX1']=self.information['xsize']/2-0.5 #### 
+        self.information['CRPIX2']=self.information['ysize']/2-0.5 #### 
+        
+        self.information['CD1_1'] =-np.cos(-theta)*self.information['pixel_size']/3600.0 #### 
+        self.information['CD1_2'] = np.sin(-theta)*self.information['pixel_size']/3600.0 #### 
+        
+        self.information['CD2_1'] = np.sin(-theta)*self.information['pixel_size']/3600.0 #### 
+        self.information['CD2_2'] = np.cos(-theta)*self.information['pixel_size']/3600.0 ####  
+        
+        obj=dict()   
+        obj['g']=np.zeros(3)            
+        obj['r']=np.zeros(3)
+        obj['i']=np.zeros(3)
+        
+        fullimg=dict()   
+        galaxy_input=dict()
+        
+        nsrcs = 72000
+        
+        galaxy_input['g']=np.zeros((nsrcs,3))
+        galaxy_input['r']=np.zeros((nsrcs,3))
+        galaxy_input['i']=np.zeros((nsrcs,3))
+        
+        galaxy_output=dict()
+        galaxy_output['g']=np.zeros((nsrcs,3))
+        galaxy_output['r']=np.zeros((nsrcs,3))
+        galaxy_output['i']=np.zeros((nsrcs,3))
+        
+        channel=['g','r','i']  
+        final_image=dict()      
+        
+        final_image['g'] = galsim.ImageF(int(self.information['xsize']), int(self.information['ysize']))
+        final_image['r'] = galsim.ImageF(int(self.information['xsize']), int(self.information['ysize']))
+        final_image['i'] = galsim.ImageF(int(self.information['xsize']), int(self.information['ysize']))
+        
+        final_image['g'].setOrigin(0,0)
+        final_image['r'].setOrigin(0,0)
+        final_image['i'].setOrigin(0,0)   
+        
+
+        ####################################################################        
+
+        nlayccd = 0 
+        
+        ####################generate loc image  
+        for k2 in (range(13)): ###   
+            print('k2=',k2)
+            
+            #  
+            filename=self.information['indata_path']+'/galaxy_Input/Lens_SED_IMG_V3_0.025/Lens_img_cut_IMG_'+str(k2+1)+'.fits'  
+           
+            if not os.path.exists(filename):
+                #print('finish load all the input galaxy image fits files already')
+                break
+                
+            srcs_cat=fits.open(filename)
+                      
+            #### load galaxy SED fitsfile  ###
+            filename=self.information['indata_path']+'/galaxy_Input/Lens_SED_IMG_V3_0.025/Lens_img_cut_SED_'+str(k2+1)+'.fits'  
+            srcs_sed=fits.open(filename)  
+            
+            ###################################################################
+            for k1 in tqdm(range(1,len(srcs_cat))):
+                
+                
+                galRa  = 3600*(srcs_cat[k1].header['new_ra'] -self.information['gal_ra'] +self.information['star_ra'])   # ra of galaxies, arcsecond
+                galDec = 3600*(srcs_cat[k1].header['new_dec']-self.information['gal_dec']+self.information['star_dec'])  # dec of galaxies, arcsecond
+
+                #################################################################
+                fsx,fsy=cal_pos(center_ra,center_dec, rotTelPos, rotSkyPos, galRa, galDec)
+                row=  fsy/losscale+self.information['ysize']/2-0.5 ### row number on CCD image, dec direction ,y direction
+                col=  fsx/losscale+self.information['xsize']/2-0.5 ### col number on CCD image, ra  direction, x direction
+                
+                if row>=self.information['ysize']+100 or row<=-100 or col>=self.information['xsize']+100 or col <=-100:
+                    continue
+                
+
+                # # SED of j-th galaxy  ,# unit of 10-17 erg/s/A/cm2         
+                gal_flux=srcs_sed[k1].data
+                
+                ################################
+                ### rotate the lensed_images_g  ###
+                if abs(theta.deg)>0:
+                    lensed_images_g=ndimage.rotate(srcs_cat[k1].data, -theta.deg, order=1, reshape=True) # here we choose reshape=False, the rotated image will
+                else:
+                    lensed_images_g=srcs_cat[k1].data
+                
+                if lensed_images_g.sum()<=0:
+                    continue
+
+                nlayccd=nlayccd+1
+
+                ################################################################
+                for i in range(3):
+                    ###                        
+                    ch=channel[i]  
+                    galaxy_input[ch][int(nlayccd)-1, 0]   =fsx       #ra
+                    galaxy_input[ch][int(nlayccd)-1, 1]   =fsy       #dec
+                    #  
+    
+                #############  do field distottion  ##########
+            
+                if self.distortion:                     
+                    
+                    for i in range(3):
+                        
+                        ch=channel[i]                    
+                        
+                        fpx,fpy=distortField(fsx, fsy, ch)                        
+                        obj[ch][0]=fpx
+                        obj[ch][1]=fpy
+                                          
+                        galaxy_output[ch][int(nlayccd)-1, 0]   =fpx       #ra
+                        galaxy_output[ch][int(nlayccd)-1, 1]   =fpy       #dec
+                                              
+                else:
+                    
+                    fpx=fsx
+                    fpy=fsy                   
+                    
+                    for i in range(3):  
+                        
+                        ch=channel[i] 
+                    
+                        obj[ch][0]=fpx
+                        obj[ch][1]=fpy                        
+                    
+                        galaxy_output[ch][int(nlayccd)-1, 0]   =fpx       #ra
+                        galaxy_output[ch][int(nlayccd)-1, 1]   =fpy       #dec
+                        
+                    
+                ####################################################################### 
+            
+                self.log.info('Galaxy number=%i, Ra(in arcsec)=%f, Dec(in arcsec)=%f' % (nlayccd, fpx, fpy))
+             
+                ##cal_SED_photons(self, flux_arr):
+                intscales,PSF_eff_weight=self.cal_SED_photons(gal_flux) 
+                            
+                img=dict()
+                
+                lensed_images_g=lensed_images_g/lensed_images_g.sum()
+                img['g']=lensed_images_g*intscales['g']
+                img['r']=lensed_images_g*intscales['r']
+                img['i']=lensed_images_g*intscales['i']
+                
+                ################################################
+                
+                gx=dict()
+                gy=dict()
+                
+                for i in range(3):
+                    
+                    ch=channel[i]                    
+                                        
+                    gy[ch]=   obj[ch][1]/losscale+self.information['ysize']/2-0.5 ### row number on CCD image
+                    gx[ch]=   obj[ch][0]/losscale+self.information['xsize']/2-0.5 ### col number on CCD image
+                        
+                ######################################################################                             
+                
+                psf=dict()     
+                
+                
+                for i in range(3):
+
+                    ch=channel[i]  
+                                     
+                    psfmat=self.get_PSF(fsx, fsy, ch)
+                    
+                    temp=0
+                    
+                    for iwave in range(7):
+                        temp=temp+PSF_eff_weight[ch][iwave]*psfmat[:,:,iwave]
+                    
+                    temp=temp/temp.sum()                        
+                                
+                    ##### rotate the PSF data
+                    if abs(theta.deg)>0:
+                        psf[ch]=ndimage.rotate(temp, theta.deg, order=1, reshape=False) # here we choose reshape=False, the rotated image will
+                    else:
+                        psf[ch]=temp
+                        
+                    conv = fftconvolve(img[ch], psf[ch], mode='full')    
+                    #suppress negative numbers
+                    conv[conv < 0.0] = 0.0                          
+                    
+                    #################################################                                                    
+                    stamp_img = galsim.Image(conv, copy=True)                        
+                    stamp_img.setOrigin(0,0)                            
+                    stamp_img.scale=0.025                             
+                                     
+                    photons=galsim.PhotonArray.makeFromImage(stamp_img, max_flux=1.0)
+                    
+                    cx0,cy0=centroid(stamp_img.array)
+                    ########### apply fat effect #####    
+                    if self.appFatt =='yes': 
+                        
+                    ### apply treering and bright fatter and diffusion;   
+                        
+                        SimpleTreeRing = galsim.SiliconSensor().simple_treerings(amplitude=self.information['treering'])                  
+                        
+                        cx, cy = centroid(conv)
+                        
+                        pos = galsim.PositionD(x=int(cx), y=int(cy)) ### set subimge center as the pos 
+                        
+                        siliconsensor = galsim.SiliconSensor(strength=self.information['fatter'], diffusion_factor=self.information['difusion'], treering_func=SimpleTreeRing, treering_center=pos)
+                                            
+                        image = galsim.Image(conv.shape[1],conv.shape[0])                            
+                        image.scale = 0.025                            
+                        image.setOrigin(0,0)
+                        
+                        photons.x=photons.x/image.scale
+                        photons.y=photons.y/image.scale
+                        
+                        siliconsensor.accumulate(photons, image)   
+                        
+                        photons=galsim.PhotonArray.makeFromImage(image, max_flux=1.0)
+                        
+                        cx0,cy0=centroid(image.array)
+                        ##############################################################
+                        ##############################################################
+                    ##################################################################    
+                    ##################################################################   
+
+                    photons.x=photons.x-cx0*0.025
+                    photons.y=photons.y-cy0*0.025
+                    
+                    photons.x=photons.x/0.05+gx[ch]
+                    photons.y=photons.y/0.05+gy[ch]
+                    
+                    photons.addTo(final_image[ch]) 
+                    
+                    ####################################################
+                ########################################################################
+
+        ##print('total los nummber on CCD is: ',nlayccd )
+    
+        #################################################################  
+        #################################################################
+                
+        fullimg['g']=final_image['g'].array
+        fullimg['r']=final_image['r'].array
+        fullimg['i']=final_image['i'].array
+    
+        return fullimg
+
+
+###############################################################################   
+    
+###############################################################################      
+
+    def generatePRNU(self, ave=1.0, sigma=0.01):
+        """
+        Creates a PRNU field image with given properties.
+
+        :return: PRNU field image
+        :rtype: ndarray
+        """
+        self.log.info('Generating a flat field...')
+        self.log.info('The flat field has mean value of 1 and a given fluctuations, usually either 1 or 2 percent defined by sigma= %d...' % sigma)
+                             
+        self.PRNU=dict()
+        
+        np.random.seed(self.filterP[self.filter_g]['seed'])
+        self.PRNU['g'] = np.random.normal(loc=ave, scale=sigma, size=(self.information['ysize'], self.information['xsize']))
+
+        np.random.seed(self.filterP[self.filter_r]['seed'])
+        self.PRNU['r'] = np.random.normal(loc=ave, scale=sigma, size=(self.information['ysize'], self.information['xsize']))
+
+        np.random.seed(self.filterP[self.filter_i]['seed'])
+        self.PRNU['i'] = np.random.normal(loc=ave, scale=sigma, size=(self.information['ysize'], self.information['xsize']))
+
+        
+        return self.PRNU
+###############################################################################
+###############################################################################
+    def MakeFlatMatrix(self, img, seed):
+        ####
+        ysize, xsize=img.shape       
+        np.random.seed(seed)
+        r1,r2,r3,r4 = np.random.random(4)
+        a1 = -0.5 + 0.2*r1
+        a2 = -0.5 + 0.2*r2
+        a3 = r3+5
+        a4 = r4+5
+        xmin,xmax,ymin,ymax = 0, xsize,0, ysize
+        Flty, Fltx = np.mgrid[ymin:ymax, xmin:xmax]        
+        np.random.seed(seed)
+        p1,p2,bg=np.random.poisson(1000, 3)
+        Fltz = 1e-6*(a1 * (Fltx-p1) ** 2 + a2 * (Flty-p2) ** 2 - a3*Fltx - a4*Flty) + bg*20
+        FlatMat = Fltz/np.mean(Fltz)
+        
+        return FlatMat
+####################################################################################################    
+    
+    def applyFlat(self):
+        
+        #### apply Flat field to image
+        self.FlatMat=dict()
+        ###
+        self.FlatMat['g']=self.MakeFlatMatrix(self.image_g, self.filterP[self.filter_g]['seed'])
+        self.image_g*=self.FlatMat['g']
+        ###
+        self.FlatMat['r']=self.MakeFlatMatrix(self.image_r, self.filterP[self.filter_r]['seed'])
+        self.image_r*=self.FlatMat['r']
+        ###
+        self.FlatMat['i']=self.MakeFlatMatrix(self.image_i, self.filterP[self.filter_i]['seed'])
+        self.image_g*=self.FlatMat['i']
+        return
+    
+######################################################################################################
+    def applyPRNUeffect(self):
+        """
+        Applies multiplicative flat field to emulate pixel-to-pixel non-uniformity.
+
+        Because the pixel-to-pixel non-uniformity effect (i.e. multiplicative) flat fielding takes place
+        before CTI and other effects, the flat field file must be the same size as the pixels that see
+        the sky. 
+        """
+        ### generate flatfiledfile, with 
+        prnu_sigma=self.information['prnu_sigma']      
+        
+        self.generatePRNU(ave=1.0, sigma=prnu_sigma) 
+        self.image_g *= self.PRNU['g']
+        self.image_r *= self.PRNU['r']
+        self.image_i *= self.PRNU['i']
+        self.log.info('Applied flatfield to images.')
+        
+##################################################################################################
+
+##################################################################################################
+
+    def addCosmicRays(self):
+        """
+        Add cosmic rays to the arrays based on a power-law intensity distribution for tracks.
+        Cosmic ray properties (such as location and angle) are chosen from random Uniform distribution.
+        For details, see the documentation for the cosmicrays class in the support package.
+        """
+        self.readCosmicRayInformation()
+        self.cr['exptime'] = self.information['exptime']  #to scale the number of cosmics with exposure time
+
+        #cosmic ray image
+        crImage = np.zeros((self.information['ysize'], self.information['xsize']), dtype=np.float64)
+
+        #cosmic ray instance
+        cosmics_g = cosmicrays.cosmicrays(self.log, crImage, crInfo=self.cr)
+        cosmics_r = cosmicrays.cosmicrays(self.log, crImage, crInfo=self.cr)
+        cosmics_i = cosmicrays.cosmicrays(self.log, crImage, crInfo=self.cr)
+
+        #add cosmic rays up to the covering fraction
+   
+        CCD_cr_g = cosmics_g.addUpToFraction(self.information['coveringfraction'], limit=None)
+        CCD_cr_r = cosmics_r.addUpToFraction(self.information['coveringfraction'], limit=None)
+        CCD_cr_i = cosmics_i.addUpToFraction(self.information['coveringfraction'], limit=None)
+
+        #paste the information
+      
+        self.image_g += CCD_cr_g
+        self.image_r += CCD_cr_r
+        self.image_i += CCD_cr_i
+
+        #save cosmic ray image map
+ 
+        self.cosmicMap_g = CCD_cr_g
+        self.cosmicMap_r = CCD_cr_r
+        self.cosmicMap_i = CCD_cr_i
+
+        #count the covering factor
+    
+        area_cr_g = np.count_nonzero(self.cosmicMap_g)
+        area_cr_r = np.count_nonzero(self.cosmicMap_r)
+        area_cr_i = np.count_nonzero(self.cosmicMap_i)
+        
+        #self.log.info('The cosmic ray covering factor is %i pixels ' % area_cr)
+        self.log.info('The cosmic ray in G channel covering factor is %i pixels ' % area_cr_g)
+        self.log.info('The cosmic ray in R channel covering factor is %i pixels ' % area_cr_r)
+        self.log.info('The cosmic ray in I channel covering factor is %i pixels ' % area_cr_i)
+
+
+########################################################################################
+
+    def ShutterEffectMat(self,  img, t_exp, t_shutter=1.3, dist_bearing=4000, dt=0.001):
+        # Generate Shutter-Effect normalized image
+        # t_shutter: time of shutter movement
+        # dist_bearing: distance between two bearings of shutter leaves
+        # dt: delta_t of sampling   
+        
+        SampleNumb = int(t_shutter/dt+1)
+        
+        DistHalf = dist_bearing/2
+        
+        t = np.arange(SampleNumb)*dt
+        a_arr = 5.84*np.sin(2*math.pi/t_shutter*t)
+        v = np.zeros(SampleNumb)
+        theta = np.zeros(SampleNumb)
+        x = np.arange(SampleNumb)/(SampleNumb-1)*dist_bearing
+        s = np.zeros(SampleNumb)
+        s1 = np.zeros(SampleNumb)
+        s2 = np.zeros(SampleNumb)
+        brt  = np.zeros(SampleNumb)
+        idx = np.arange(SampleNumb)
+        sidx = np.zeros(SampleNumb)
+        s1idx = np.zeros(SampleNumb)
+        s2idx = np.zeros(SampleNumb)
+        
+        v[0] = 0
+        theta[0] = 0
+        
+        for i in range(SampleNumb-1):
+            v[i+1] = v[i]+a_arr[i]*dt
+            theta[i+1] = theta[i]+v[i]*dt
+            s1[i] = DistHalf*np.cos(theta[i])
+            s2[i] = dist_bearing-DistHalf*np.cos(theta[i])
+            s1idx[i] = int(s1[i]/dist_bearing*(SampleNumb))
+            s2idx[i] = int(s2[i]/dist_bearing*(SampleNumb))
+            brt[(idx>s1idx[i]) & (idx<s2idx[i])] += dt
+        
+        if t_exp>t_shutter*2:
+            brt = brt*2+(t_exp-t_shutter*2)
+        else:
+            brt = brt*2
+        
+        x = (x-dist_bearing/2)*100
+        
+        intp = interpolate.splrep(x, brt, s=0)
+        
+        xmin = 0    # 
+        xmax = img.shape[1]  # row
+        ymin = 0    # 
+        ymax = img.shape[0]  # column
+        if xmin<np.min(x) or xmax>np.max(x):
+            raise LookupError("Out of focal-plane bounds in X-direction.")
+        if ymin<0 or ymax>25331:
+            raise LookupError("Out of focal-plane bounds in Y-direction.")
+            
+        sizex = xmax-xmin
+        sizey = ymax-ymin
+        xnewgrid = np.mgrid[xmin:(xmin+sizex)]
+        expeffect = interpolate.splev(xnewgrid, intp, der=0)
+        expeffect /= np.max(expeffect)
+        exparrnormal = np.tile(expeffect, (sizey,1))
+        
+        # Image *= exparrnormal
+        
+        return exparrnormal
+    
+########################################################################################    
+    def addshuttereffect(self):
+        ''' apply shutter effect to image'''
+        #          
+        self.shutterMat=self.ShutterEffectMat(self.image_g, self.information['exptime'])
+        self.image_g *= self.shutterMat
+        self.image_r *= self.shutterMat
+        self.image_i *= self.shutterMat
+        return
+        
+
+###############################################################################
+###############################################################################       
+
+
+    def applyDarkCurrent(self):
+        """
+        Apply dark current. Scales the dark with the exposure time.
+
+        Additionally saves the image without noise to a FITS file.
+        """
+
+        #add dark
+        dark = self.information['exptime'] * self.information['dark']
+
+        self.image_g += dark
+        self.image_r += dark
+        self.image_i += dark
+        
+        self.log.info('Added dark current = %f' % dark)
+        return
+
+    ###########################################################################
+
+
+    ###########################################################################
+    def applyPoissonSkyNoise(self):
+        """
+        Add Poisson sky noise to the image.
+        """
+        
+        skynoise_g=self.Sky_Noise['g']*np.ones_like(self.image_g)
+        np.random.seed()
+        self.image_g =self.image_g+ np.random.poisson(lam=skynoise_g)
+        self.log.info('Added Poisson sky noise on channel g')
+
+        skynoise_r=self.Sky_Noise['r']*np.ones_like(self.image_r)
+        np.random.seed()
+        self.image_r =self.image_r+ np.random.poisson(lam=skynoise_r)
+        self.log.info('Added Poisson sky noise on channel r')
+        
+        skynoise_i=self.Sky_Noise['i']*np.ones_like(self.image_i)
+        np.random.seed()
+        self.image_i =self.image_i+ np.random.poisson(lam=skynoise_i)
+        self.log.info('Added Poisson sky noise on channel i')
+        
+        return
+
+
+###############################################################################
+
+###############################################################################
+    def applyCosmetics(self):
+        """
+        Apply cosmetic defects described in the input file.
+       
+        #Number of hot and dead pixels from MSSL/Euclid/TR/12003 Issue 2 Draft b
+
+        .. Warning:: This method does not work if the input file has exactly one line.
+        """
+        #######################################################################
+        cosmetics = np.loadtxt(self.information['indata_path']+'/data/Cosmetics_g.txt')
+
+        x = np.round(cosmetics[:, 0]).astype(int) ## row number
+        y = np.round(cosmetics[:, 1]).astype(int) ## col number
+        value = np.round(cosmetics[:, 1]).astype(int)
+
+        #cosmetics_g=np.zeros((4636,235526)) 
+        
+        self.log.info('Adding cosmetic defects to G channel:' )
+        for xc, yc, val in zip(x, y, value):
+            if 0 < xc < 4616  and  27 < (yc % 1499) <1179:
+                
+                self.image_g[xc, yc] = val  
+                #cosmetics_g[yc,xc]=val
+                self.log.info('x=%i, y=%i, value=%f' % (xc, yc, val))
+        ######################################################################
+        #######################################################################
+        cosmetics = np.loadtxt(self.information['indata_path']+'/data/Cosmetics_r.txt')
+
+        x = np.round(cosmetics[:, 0]).astype(int) ## row number
+        y = np.round(cosmetics[:, 1]).astype(int) ## col number
+        value = np.round(cosmetics[:, 1]).astype(int)
+
+        #cosmetics_r=np.zeros((4636,235526)) 
+        
+        self.log.info('Adding cosmetic defects to R channel:' )
+        for xc, yc, val in zip(x, y, value):
+            if 0 < xc < 4616  and  27 < (yc % 1499) <1179:
+                
+                self.image_r[xc, yc] = val  
+                #cosmetics_r[yc,xc]=val
+                self.log.info('x=%i, y=%i, value=%f' % (xc, yc, val))
+##############################################################################
+        #######################################################################
+        cosmetics = np.loadtxt(self.information['indata_path']+'/data/Cosmetics_i.txt')
+
+        x = np.round(cosmetics[:, 0]).astype(int) ## row number
+        y = np.round(cosmetics[:, 1]).astype(int) ## col number
+        value = np.round(cosmetics[:, 1]).astype(int)
+
+        #cosmetics_i=np.zeros((4636,235526)) 
+        
+        self.log.info('Adding cosmetic defects to I channel:' )
+        for xc, yc, val in zip(x, y, value):
+            if 0 < xc < 4616  and  27 < (yc % 1499) <1179:
+                
+                self.image_i[xc, yc] = val  
+                #cosmetics_i[yc,xc]=val
+                self.log.info('x=%i, y=%i, value=%f' % (xc, yc, val))
+        return
+     
+#############################################################################################################################        
+    def applyRadiationDamage(self):
+        """
+        Applies CDM03 radiation model to the image being constructed.
+
+        .. seealso:: Class :`CDM03`
+        """
+
+            
+        #save image without CTI
+        self.noCTI_g = self.image_g.copy()
+       
+        self.log.info('Save NoCti Nonoise image fits file to  %s ' % ('noctinonoise' + 'mci_g.fits'))
+        
+        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_g = cti.applyRadiationDamage(self.image_g.copy().transpose(), iquadrant=self.information['quadrant']).transpose()
+        self.log.info('Radiation damage added.')
+
+        #################################################
+            
+        self.noCTI_r = self.image_r.copy()
+       
+        self.log.info('Save NoCti Nonoise image fits file to  %s ' % ('noctinonoise' + 'mci_r.fits'))
+        
+        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.')
+
+        ##################################################
+            
+        self.noCTI_i = self.image_i.copy()
+        
+        self.log.info('Save NoCti Nonoise image fits file to  %s ' % ('noctinonoise' + 'mci_i.fits'))
+        
+        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_i = cti.applyRadiationDamage(self.image_i.copy().transpose(), iquadrant=self.information['quadrant']).transpose()
+        self.log.info('Radiation damage added.')
+        return
+
+
+    ###########################################################################
+
+    def applyNonlinearity(self):
+        """
+        Applies a CCD273 non-linearity model to the image being constructed.
+        
+        """
+
+        self.log.debug('Starting to apply non-linearity model...')
+        self.image_g = MCIinstrumentModel.CCDnonLinearityModel(self.image_g.copy())
+        self.log.info('Non-linearity effects included.')
+
+        ########################################################################
+        self.log.debug('Starting to apply non-linearity model...')
+        self.image_r = MCIinstrumentModel.CCDnonLinearityModel(self.image_r.copy())
+        self.log.info('Non-linearity effects included.')
+
+        ########################################################################        
+        self.log.debug('Starting to apply non-linearity model...')
+        self.image_i = MCIinstrumentModel.CCDnonLinearityModel(self.image_i.copy())
+
+        self.log.info('Non-linearity effects included.')
+        return
+
+###############################################################################
+
+    def applyReadoutNoise(self):
+        
+        """
+        Applies readout noise to the image being constructed.
+
+        The noise is drawn from a Normal (Gaussian) distribution with average=0.0 and std=readout noise.
+        """       
+        if self.information['xsize']==9216 and self.information['ysize']==9232 :
+            xmin=dict()
+            xmax=dict()
+            ymin=dict()
+            ymax=dict()
+          
+            ############## calculate the non-over subimg of the final matrix of the image 
+            for k in range(1,17):
+                xmin[k]=0
+                xmax[k]=4616+self.overscan
+                 
+                ymin[k]=(1152+self.prescan+self.overscan)*(k-1)
+                ymax[k]=(1152+self.prescan+self.overscan)*(k)
+                ###############   
+                np.random.seed()               
+                noise_g = np.random.normal(loc=0.0, scale=self.information['g_rdnois'+str(k)], size=(4616+self.overscan,1152+self.overscan+self.prescan))
+                self.image_g[xmin[k]:xmax[k], ymin[k]:ymax[k]]+=noise_g
+                self.log.info('Bias of %i counts were added to the g band image' % self.information['g_rdnois'+str(k)])
+                 
+                np.random.seed()               
+                noise_r = np.random.normal(loc=0.0, scale=self.information['r_rdnois'+str(k)], size=(4616+self.overscan,1152+self.overscan+self.prescan))
+                self.image_r[xmin[k]:xmax[k], ymin[k]:ymax[k]]+=noise_r
+                self.log.info('Bias of %i counts were added to the r band image' % self.information['r_rdnois'+str(k)])
+                 
+                np.random.seed()               
+                noise_i = np.random.normal(loc=0.0, scale=self.information['i_rdnois'+str(k)], size=(4616+self.overscan,1152+self.overscan+self.prescan))
+                self.image_i[xmin[k]:xmax[k], ymin[k]:ymax[k]]+=noise_i
+                self.log.info('Bias of %i counts were added to the i band image' % self.information['i_rdnois'+str(k)])
+                
+                ### end for     
+ 
+        else:
+            
+            np.random.seed()
+            noise = np.random.normal(loc=0.0, scale=4, size=self.image_g.shape)
+            self.log.info('Sum of readnoise in g channel = %f' % np.sum(noise))    
+            #add to the image
+            self.image_g += noise   
+             
+            #######################################################    
+            np.random.seed()            
+            noise = np.random.normal(loc=0.0, scale=4, size=self.image_r.shape)
+            self.log.info('Sum of readnoise in r channel= %f' % np.sum(noise))    
+             #add to the image
+            self.image_r += noise 
+                      
+            ######################################## 
+            np.random.seed()               
+            noise = np.random.normal(loc=0.0, scale=4, size=self.image_i.shape)
+            self.log.info('Sum of readnoise in i channel= %f' % np.sum(noise))    
+            #add to the image
+            self.image_i += noise
+            
+        return
+        
+##########################################################################################################        
+            
+    def electrons2ADU(self):
+        """
+        Convert from electrons to ADUs using the value read from the configuration file.
+        """
+        #####
+        if self.overscans:
+            ######################         
+            xmin=dict()
+            xmax=dict()
+            ymin=dict()
+            ymax=dict()
+          
+           ############## calculate the non-over subimg of the final matrix of the image 
+            for k in range(1,17):
+                xmin[k]=0
+                xmax[k]=4616+self.overscan
+                
+                ymin[k]=(1152+self.overscan+self.prescan)*(k-1)
+                ymax[k]=(1152+self.overscan+self.prescan)*(k)
+                ###############
+                               
+                self.image_g[xmin[k]:xmax[k], ymin[k]:ymax[k]]/=self.information['g_gain'+str(k)]
+                                
+                self.image_r[xmin[k]:xmax[k], ymin[k]:ymax[k]]/=self.information['r_gain'+str(k)]
+                
+                self.image_i[xmin[k]:xmax[k], ymin[k]:ymax[k]]/=self.information['i_gain'+str(k)]
+                
+            ### end for    
+            return
+                  
+ ####################################################################################     
+            
+    def applyBias(self):
+        """
+        Adds a bias level to the image being constructed.
+
+        The value of bias is read from the configure file and stored
+        in the information dictionary (key bias).
+        
+        """
+       
+        xmin=dict()
+        xmax=dict()
+        ymin=dict()
+        ymax=dict()
+    
+        ############## calculate the non-over subimg of the final matrix of the image 
+        for k in range(1,17):
+            xmin[k]=0
+            xmax[k]=4616+self.overscan
+             
+            ymin[k]=(1152+self.prescan+self.overscan)*(k-1)
+            ymax[k]=(1152+self.prescan+self.overscan)*(k)
+             ###############                               
+            self.image_g[xmin[k]:xmax[k], ymin[k]:ymax[k]]+=self.information['g_detbia'+str(k)]
+            self.log.info('Bias of %i counts were added to the g band image' % self.information['g_detbia'+str(k)])
+
+            self.image_r[xmin[k]:xmax[k], ymin[k]:ymax[k]]+=self.information['r_detbia'+str(k)]
+            self.log.info('Bias of %i counts were added to the g band image' % self.information['r_detbia'+str(k)])
+    
+            self.image_i[xmin[k]:xmax[k], ymin[k]:ymax[k]]+=self.information['i_detbia'+str(k)]
+            self.log.info('Bias of %i counts were added to the g band image' % self.information['i_detbia'+str(k)])
+      
+            ### end for     
+        return
+                      
+##############################################################################      
+####################################################################################
+
+    def applyBleeding(self, img, direction='horizon'):
+        """
+        Apply bleeding along the CCD readout direction if the number of electrons in a pixel exceeds the full-well capacity.
+    
+       direction: this is the CCD readout direction, and default value is horizon direction
+    
+        :return: the bleeding image
+        """
+     
+        data=img.copy()
+        
+        if data.max()<self.information['fullwellcapacity']:
+            return img
+        
+        else:
+            
+            if direction=='horizon':
+                
+                #loop over each column, as bleeding is modelled column-wise
+                for i, column in enumerate(data):  ### select one solumnn
+                    if column.max() <=self.information['fullwellcapacity']:
+                        continue
+                    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[i, j] -= overload    
+                            sum += overload
+                            
+                        elif sum > 0.:
+                            if -overload > sum:
+                                overload = -sum
+                            #self.image[j, i] -= overload
+                            data[i,j] -= overload                
+                            sum += overload
+                ######################################################
+                for i, column in enumerate(data):
+                    if column.max() <=self.information['fullwellcapacity']:
+                        continue
+                    sum = 0.
+                    for j, value in enumerate(column[::-1]):
+                        #second round - from top to bottom (bleeding was half'd already, so now full)
+                        overload = value - self.information['fullwellcapacity']
+                        if overload > 0.:
+                            #self.image[-j-1, i] -= overload
+                            data[ i,-j-1] -= overload                     
+                            sum += overload
+                            
+                        elif sum > 0.:
+                            if -overload > sum:
+                                overload = -sum
+                            #self.image[-j-1, i] -= overload
+                            data[i,-j-1,] -= overload                
+                            sum += overload
+                        
+            else:        
+           
+                #loop over each column, as bleeding is modelled column-wise
+                for i, column in enumerate(data.T):
+                    sum = 0.
+                    for j, value in enumerate(column):
+                        #first round - from bottom to top (need to half the bleeding)
+                        overload = value - self.information['fullwellcapacity']
+                        if overload > 0.:
+                            overload /= 2.
+                            #self.image[j, i] -= overload
+                            data[j, i] -= overload
+            
+                            sum += overload
+                        elif sum > 0.:
+                            if -overload > sum:
+                                overload = -sum
+                            #self.image[j, i] -= overload
+                            data[j, i] -= overload                
+                            sum += overload
+            ################################
+                for i, column in enumerate(data.T):
+                    sum = 0.
+                    for j, value in enumerate(column[::-1]):
+                        #second round - from top to bottom (bleeding was half'd already, so now full)
+                        overload = value - self.information['fullwellcapacity']
+                        if overload > 0.:
+                            #self.image[-j-1, i] -= overload
+                            data[-j-1, i] -= overload
+                             
+                            sum += overload
+                        elif sum > 0.:
+                            if -overload > sum:
+                                overload = -sum
+                            #self.image[-j-1, i] -= overload
+                            data[-j-1, i] -= overload                
+                            sum += overload
+                        
+            ######print('Applying column bleeding finished.......')
+            return data
+        
+    ############################################################################
+
+    def discretise(self, max=2**16-1):
+        """
+        Converts a floating point image array (self.image) to an integer array with max values
+        defined by the argument max.
+
+        :param max: maximum value the the integer array may contain [default 65k]
+        :type max: float
+
+        :return: None
+        """
+
+        #avoid negative numbers in case bias level was not added
+        self.image_g[self.image_g < 0.0] = 0.
+        #cut of the values larger than max
+        self.image_g[self.image_g > max] = max
+
+        self.image_g = np.rint(self.image_g).astype(int)
+        self.log.info('Maximum and total values of the image are %i and %i, respectively' % (np.max(self.image_g),
+                                                                                             np.sum(self.image_g)))
+
+        #avoid negative numbers in case bias level was not added
+        self.image_r[self.image_r < 0.0] = 0.
+        #cut of the values larger than max
+        self.image_r[self.image_r > max] = max
+
+        self.image_r = np.rint(self.image_r).astype(int)
+        self.log.info('Maximum and total values of the image are %i and %i, respectively' % (np.max(self.image_r),
+                                                                                             np.sum(self.image_r)))
+  
+        ###############################################################################
+
+        #avoid negative numbers in case bias level was not added
+        self.image_i[self.image_i < 0.0] = 0.
+        #cut of the values larger than max
+        self.image_i[self.image_i > max] = max
+
+        self.image_i = np.rint(self.image_i).astype(int)
+        self.log.info('Maximum and total values of the image are %i and %i, respectively' % (np.max(self.image_i),
+                                                                                             np.sum(self.image_i)))
+        return
+    ##################################################################################################   
+    def addScanEffect(self, img):
+        #### read CCD image to matrix as defined data format
+        '''
+        function: add prescan and overscan empty zone to the CCD physical image
+        
+        return:  img
+        
+        '''
+        
+        ### xsize is column, ysize is row
+        
+        self.prescan =self.information['prescan'] ## horizon direction , columnn
+        self.overscan=self.information['overscan']    ## vertical direction, row
+        
+        row = self.information['ysize']
+        col = self.information['xsize']
+        if row !=9232 or  col !=9216:   #####
+            print('image size is not correct.......')
+            sys.exit(1)
+        ####
+        ###########################################
+        xmin=dict()
+        xmax=dict()
+        ymin=dict()
+        ymax=dict()
+        #### physical CCD zone 1-16
+        for n in range(1,17):
+            
+            #### physical zone index
+            if n<9:
+                xmin[n]=0           # row number
+                xmax[n]=4616
+                
+                ymin[n]=(n-1)*1152     # col number
+                ymax[n]=ymin[n]+1152
+            else:
+                xmin[n]=4616
+                xmax[n]=4616*2
+                
+                ymin[n]=(16-n)*1152
+                ymax[n]=ymin[n]+1152
+        ###########################################
+
+        temp=np.zeros((int(4616+self.overscan), int((1152+self.overscan+self.prescan)*16)))
+        
+        matxmin=dict()
+        matxmax=dict()
+        matymin=dict()
+        matymax=dict()
+        
+############## get the non-over subimg of the final matrix of the image 
+        for k in range(1,17):
+            
+            ### big new image include prescan and overscan;
+            matxmin[k]=0       ### row number
+            matxmax[k]=4616    ### row number
+            
+            matymin[k]=(1152+self.overscan+self.prescan)*(k-1)+self.prescan  ## col number start
+            matymax[k]= matymin[k]+1152   ## col number end
+            
+        ################################
+        
+        for k in range(1,17):
+            if k<=4:  ### zos1 zos2 zos3 zos4
+                subimg=img[xmin[k]:xmax[k],ymin[k]:ymax[k]] # get k-th subimg from input img
+                subimg=np.flipud(subimg)  ## down to up, up to down  
+                
+              
+                
+            if k>=5 and k<=8:  ### zos5 zos6 zos7 zos8
+                subimg=img[xmin[k]:xmax[k],ymin[k]:ymax[k]] # get k-th subimg from input img
+                subimg=np.flipud(subimg)  ## down to up, up to down 
+                subimg=np.fliplr(subimg)  ## left to right  ,right to left
+                
+              
+                
+            if k>=9 and k<=12:  ### zos5 zos6 zos7 zos8
+                subimg=img[xmin[k]:xmax[k],ymin[k]:ymax[k]] # get k-th subimg from input img
+                subimg=np.fliplr(subimg)  ## left to right  ,right to left
+                
+             
+                
+            if k>=13 and k<=16:  ### zos5 zos6 zos7 zos8
+                subimg=img[xmin[k]:xmax[k],ymin[k]:ymax[k]] # get k-th subimg from input img
+                
+             
+            ######################
+            
+            temp[matxmin[k]:matxmax[k],matymin[k]:matymax[k]]=subimg[:,:]
+        #### end for
+        return temp
+    ###########################################################################
+    
+    
+    def applyPoissonNoise(self):
+        """
+        Add Poisson noise to the image.
+        """
+
+        
+        rounded = np.rint(self.image_g)   ### round to 
+        residual = self.image_g.copy() - rounded #ugly workaround for multiple rounding operations...
+        rounded[rounded < 0.0] = 0.0
+        self.image_g = np.random.poisson(rounded).astype(np.float64)
+        self.log.info('Added Poisson noise on channel g')
+        self.image_g += 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
+        self.image_r = np.random.poisson(rounded).astype(np.float64)
+        self.log.info('Added Poisson noise on channel r')
+        self.image_r += residual
+        
+        rounded = np.rint(self.image_i)   ### round to 
+        residual = self.image_i.copy() - rounded #ugly workaround for multiple rounding operations...
+        rounded[rounded < 0.0] = 0.0
+        self.image_i = np.random.poisson(rounded).astype(np.float64)
+        self.log.info('Added Poisson noise on channel i')
+        self.image_i += residual
+        
+        
+    
+    #########################################################################################    
+
+
+    def writeOutputs(self, obnum=1):
+        """
+        Writes out a FITS file using PyFITS and converts the image array to 16bit unsigned integer as
+        appropriate for VIS.
+
+        Updates header with the input values and flags used during simulation.
+        """
+
+                ## Readout information
+
+        ##########################################################################
+        obsid=200000000+obnum
+        #create a new FITS file, using HDUList instance
+        ofd_g = fits.PrimaryHDU()
+        #### add primary header
+
+        exp_starttime=self.dt.strftime("%Y%m%d%H%M%S")
+        
+        ### exposure end time is t2 ;
+        t2=self.dt+timedelta(seconds=self.information['exptime'])
+        
+        exp_endtime=t2.strftime("%Y%m%d%H%M%S")
+        
+        
+        ### data read time is the exposure end time ###
+        ### data read end time is t3
+        t3=self.dt+timedelta(seconds=self.information['exptime'])+timedelta(seconds=self.information['readouttime'])
+        
+        filename_g='CSST_MCI_C1_'+self.source+'_'+exp_starttime+'_'+exp_endtime+'_'+str(obsid)+'_'+ \
+                   str(self.filterP[self.filter_g]['number'])+'_L0_'+self.information['ver']+'.fits'
+        
+         
+        ofd_g.header['GROUPS']=( bool(False),   'always F')
+        ofd_g.header['DATE']  =( t3.strftime("%Y-%m-%d %H:%M:%S"),   'date this file was written' )   
+        ofd_g.header['FILENAME']=(filename_g,   ' file name C48  ')        
+        ofd_g.header['OBSTYPE'] =( self.source,   'observation type raw,flt, mask, bias, dark, sci')        
+        ofd_g.header['TELESCOP']=('CSST',   'always CSST')
+        ofd_g.header['INSTRUME']=( 'MCI',   '   ')
+        ofd_g.header['RADECSYS']=('ICRS',   ' always ICRS ')
+        ofd_g.header['EQUINOX'] =( float(2000.0),   'always 2000.0')
+        ofd_g.header['FITSCREA']=( '4.2.1',   'FITS create software version')
+        
+        #########   Object information  #############             
+        
+        ofd_g.header['OBJECT']=( self.information['name_obj'],   'object name')
+        ofd_g.header['TARGET']=( self.information['target'],     'target name, hhmmss+ddmmss')
+        ofd_g.header['OBJ_RA'] =(np.float64(self.information['ra_obj']) ,   'RA of the object in deg')
+        ofd_g.header['OBJ_DEC']=(np.float64(self.information['dec_obj']) ,   'DEC of the object in deg')
+        
+        ofd_g.header['RA_PNT0']=( np.float64(self.information['ra_pnt0']) ,     'RA of the pointing (degrees) at EXPSTART')
+        ofd_g.header['DEC_PNT0']=(np.float64(self.information['dec_pnt0']) ,    'DEC of the pointing (degrees) at EXPSTART')
+    
+
+        
+        ##############
+        ofd_g.header['OBSID'] =(str(obsid) ,   'observation ID, 3+8bit')
+        
+        ########  Telescope information ###############
+        # ofd_g.header['COMMENT'] ='=========================================================='
+        # ofd_g.header['COMMENT'] ='Telescope information'
+        # ofd_g.header['COMMENT'] ='=========================================================='       
+        
+        ofd_g.header['REFFRAME']=('CSSTGSC-1.0' ,   'guide star catalog version')
+        ofd_g.header['DATE-OBS']=(self.dt.strftime("%Y-%m-%d %H:%M:%S") ,   'date of the observation (yyyy-mm-dd hh:mm:ss)')
+        
+                
+        ofd_g.header['EXPSTART']=(time2jd(self.dt),   'exposure start time, UTC')
+        ofd_g.header['SUNANGL0']=(np.float32(0.0) ,   'angle between sun and optical axis at EXPSTART')
+        ofd_g.header['MOONANG0']=(np.float32(0.0) ,   'angle between moon and optical axis at EXPSTART')
+        ofd_g.header['POS_ANG0']=(np.float64(0.0),   'angle between optical axis and the North Pole at EXPSTART in arcsec')
+        ofd_g.header['TEL_ALT0']=(np.float64(0.0),   'angle between optical axis and the ground- piston at EXPSTART in deg')
+        
+        ofd_g.header['HOODSTA0']=(np.float32(0.0) ,   'lens hood altitude at EXPSTART')
+        ofd_g.header['HOODANG0']=(np.float32(0.0),    'lens hood azimuth at EXPSTART')
+        ofd_g.header['POSI0_X'] =(np.float64(self.information['pos_x']) ,   'the orbital position of CSST in X direction at EXPSTART')
+        ofd_g.header['POSI0_Y'] =(np.float64(self.information['pos_y']) ,   'the orbital position of CSST in Y direction at EXPSTART')
+        ofd_g.header['POSI0_Z'] =(np.float64(self.information['pos_z']) ,   'the orbital position of CSST in Z direction at EXPSTART')
+        ofd_g.header['VELO0_X']=( np.float64(self.information['velocity_x']) ,   'the orbital velocity of CSST in X direction at EXPSTART')
+        ofd_g.header['VELO0_Y']=( np.float64(self.information['velocity_y']) ,   'the orbital velocity of CSST in Y direction at EXPSTART')
+        ofd_g.header['VELO0_Z']=( np.float64(self.information['velocity_z']) ,   'the orbital velocity of CSST in Z direction at EXPSTART')
+        
+        ofd_g.header['Euler0_1']=( np.float64(0.0),   'Euler angle 1 at EXPSTART')
+        ofd_g.header['Euler0_2']=( np.float64(0.0),   'Euler angle 2 at EXPSTART')
+        ofd_g.header['Euler0_3']=( np.float64(0.0),   'Euler angle 3 at EXPSTART')
+        
+   
+        ofd_g.header['EXPEND']  =(time2jd(t2) ,       'exposure end time,UTC')
+        
+        ofd_g.header['SUNANGL1']=(np.float32(0.0),    'angle between sun and optical axis at EXPEND')
+        ofd_g.header['MOONANG1']=(np.float32(0.0) ,   'angle between moon and optical axis at EXPEND   ')
+        ofd_g.header['POS_ANG1']=(np.float64(0.0) ,   'angle between optical axis and the North Pole at EXPEND in arcsec')
+        ofd_g.header['TEL_ALT1']=(np.float64(0.0) ,   'angle between optical axis and the ground- piston at EXPEND in deg  ')
+        ofd_g.header['HOODSTA1']=(np.float32(0.0),    'lens hood altitude at EXPEND   ')
+        ofd_g.header['HOODANG1']=(np.float32(0.0),    'lens hood azimuth at EXPEND   ')
+        
+        
+        t2jd=time2jd(t2)
+        
+        if self.orbit_pars[-1,0]<t2jd: ## orbit parameters are not in currenct txt file
+            self.orbit_file_num=self.orbit_file_num+1
+            fn=self.information['indata_path']+'/TianCe/orbit20160925/'+str(self.orbit_file_num)+'.txt';
+            self.orbit_pars=np.loadtxt(fn);
+            self.orbit_exp_num=0
+            
+            
+            
+        for k in range(self.orbit_exp_num, len(self.orbit_pars),1):
+                        
+            if t2jd-self.orbit_pars[k,0]<0:
+                break
+        if k==0:
+            deltaT=jd2time(self.orbit_pars[k,0])-t2
+            p1x=self.orbit_pars[k,1]-(self.orbit_pars[k+1,1]-self.orbit_pars[k,1])*deltaT.seconds/120
+            p1y=self.orbit_pars[k,2]-(self.orbit_pars[k+1,2]-self.orbit_pars[k,2])*deltaT.seconds/120
+            p1z=self.orbit_pars[k,3]-(self.orbit_pars[k+1,3]-self.orbit_pars[k,3])*deltaT.seconds/120
+            
+            p1vx=self.orbit_pars[k,4]-(self.orbit_pars[k+1,4]-self.orbit_pars[k,4])*deltaT.seconds/120
+            p1vx=self.orbit_pars[k,5]-(self.orbit_pars[k+1,5]-self.orbit_pars[k,5])*deltaT.seconds/120
+            p1vx=self.orbit_pars[k,6]-(self.orbit_pars[k+1,6]-self.orbit_pars[k,6])*deltaT.seconds/120
+            
+        else:
+            deltaT=jd2time(self.orbit_pars[k,0])-t2
+            p1x=self.orbit_pars[k-1,1]+(self.orbit_pars[k,1]-self.orbit_pars[k-1,1])*deltaT.seconds/120
+            p1y=self.orbit_pars[k-1,2]+(self.orbit_pars[k,2]-self.orbit_pars[k-1,2])*deltaT.seconds/120
+            p1z=self.orbit_pars[k-1,3]+(self.orbit_pars[k,3]-self.orbit_pars[k-1,3])*deltaT.seconds/120
+            
+            p1vx=self.orbit_pars[k-1,4]+(self.orbit_pars[k,4]-self.orbit_pars[k-1,4])*deltaT.seconds/120
+            p1vy=self.orbit_pars[k-1,5]+(self.orbit_pars[k,5]-self.orbit_pars[k-1,5])*deltaT.seconds/120
+            p1vz=self.orbit_pars[k-1,6]+(self.orbit_pars[k,6]-self.orbit_pars[k-1,6])*deltaT.seconds/120
+            
+        
+        ofd_g.header['POSI1_X'] =(np.float64(p1x) ,   'the orbital position of CSST in X direction at EXPEND')
+        ofd_g.header['POSI1_Y'] =(np.float64(p1y) ,   'the orbital position of CSST in Y direction at EXPEND')
+        ofd_g.header['POSI1_Z'] =(np.float64(p1z) ,   'the orbital position of CSST in Z direction at EXPEND')
+        
+        ofd_g.header['VELO1_X']=(np.float64(p1vx) ,   'the orbital velocity of CSST in X direction at EXPEND')
+        ofd_g.header['VELO1_Y']=(np.float64(p1vy) ,   'the orbital velocity of CSST in Y direction at EXPEND')
+        ofd_g.header['VELO1_Z']=(np.float64(p1vz) ,   'the orbital velocity of CSST in Z direction at EXPEND')
+        
+        ofd_g.header['Euler1_1']=( np.float64(0.0),   'Euler angle 1 at EXPEND')
+        ofd_g.header['Euler1_2']=( np.float64(0.0),   'Euler angle 2 at EXPEND')
+        ofd_g.header['Euler1_3']=( np.float64(0.0),   'Euler angle 3 at EXPEND')
+        
+        ofd_g.header['RA_PNT1']=(np.float64(0.0),    'RA of the pointing (degrees) at EXPEND in deg')
+        ofd_g.header['DEC_PNT1']=(np.float64(0.0),   'DEC of the pointing (degrees) at EXPEND in deg')
+        
+        ofd_g.header['EXPTIME']=(self.information['exptime'],   'exposure duration')
+        ofd_g.header['EPOCH']  =(np.float32(0.0),   'coordinate epoch')
+        ofd_g.header['CHECKSUM']=(   0           ,   'hdu-checksum')  
+    
+        ########## finish header for 0 layer ####################        
+    
+        #new image HDU
+        hdu_g = fits.ImageHDU(data=np.uint16(self.image_g))       
+        
+        
+        ###  MCI  header for immage layer ###
+        #### Instrument information (at EXPSTART time)
+        hdu_g.header['PMIRRPOS']=( bool(True) ,  'FSM pointing,T: point to MCI, F: not point to MCI   ')
+        hdu_g.header['DIFFUSER']=( bool(True) ,  'insert diffuser for flat calibration,T: insert diffuser, F: not insert diffuser   ')
+        hdu_g.header['FLAMP']   =( np.int16(0)  ,   'status of flat lamp,0: off, 1: 1, 2: 2, 3: 3')
+        hdu_g.header['MCITEMP'] =( np.float32(0.0), 'MCI components temperature')
+        hdu_g.header['MCISTAT'] =( np.int16(0)  ,   'MCI components status parameters')
+        
+        #####  Filter information  #####
+        hdu_g.header['FILTER'] =(self.filter_g , 'filter band name')
+        
+        #####  Detector information ####
+        hdu_g.header['DETNAM'] =('C1'     , 'detector name for channle 1')
+        hdu_g.header['DETTEMP'] =(np.float32(0)     , 'detector temperature')
+        hdu_g.header['DETSIZE'] =(str(self.image_g.shape[0])+'*'+str(self.image_g.shape[1]), 'detector size')
+        hdu_g.header['DATASEC'] =(str(self.information['ysize'])+'*'+str( self.information['xsize']), 'data section')
+        hdu_g.header['PIXSCAL1']=(np.float32(0.05)     , 'pixel scale for axis 1')
+        hdu_g.header['PIXSCAL2']=(np.float32(0.05)     , 'pixel scale for axis 2')
+        hdu_g.header['PIXSIZE1']=(np.float32(10.0)     , 'pixel size for axis 1')
+        hdu_g.header['PIXSIZE2']=(np.float32(10.0)     , 'pixel size for axis 2')
+        hdu_g.header['NCHAN']   =(np.int16(16)     , 'number of readout channels')
+        hdu_g.header['NCHAN1']  =(np.int16(8)       , 'number of horizontal channels')
+        hdu_g.header['NCHAN2']  =(np.int16(2)       , 'number of vertical channels')
+        hdu_g.header['PSCAN1']  =(np.int32(27)      , 'horizontal prescan width, per readout channel')
+        hdu_g.header['PSCAN2']  =(np.int32(8)       , 'vertical prescan height, per readout channel')
+        hdu_g.header['OSCAN1']  =(np.int32(16)      , 'horizontal overscan width, per readout channel')
+        hdu_g.header['OSCAN2']  =(np.int32(8)       , 'vertical overscan width, per readout channel')
+        hdu_g.header['BIN_X']   =(np.int16(0)       , 'bin number in X (wavelength)')
+        hdu_g.header['BIN_Y']   =(np.int16(0)       , 'bin number in Y (spatial)')
+        
+        ####   World coordinate system and related parameters  #####
+        
+        hdu_g.header['WCSAXES']=(np.int16(2)        ,    'number of World Coordinate System axes')
+        hdu_g.header['CRPIX1'] =( np.float64(self.information['CRPIX1'])       , 'x-coordinate of reference pixel')
+        hdu_g.header['CRPIX2'] =( np.float64(self.information['CRPIX2'])       , 'y-coordinate of reference pixel')
+        hdu_g.header['CRVAL1'] =( np.float64(self.information['CRVAL1'])       , 'first axis value at reference pixel')
+        hdu_g.header['CRVAL2'] =( np.float64(self.information['CRVAL2'])       , 'second axis value at reference pixel')
+        hdu_g.header['CTYPE1'] =( 'RA---TAN'                               , 'the coordinate type for the first axis')
+        hdu_g.header['CTYPE2'] =( 'DEC---TAN'                              , 'the coordinate type for the second axis')
+        hdu_g.header['CD1_1']  =( np.float64(self.information['CD1_1'])        , 'partial of first axis coordinate w.r.t. x')
+        hdu_g.header['CD1_2']  =( np.float64(self.information['CD1_2'])        , 'partial of first axis coordinate w.r.t. y')
+        hdu_g.header['CD2_1']  =( np.float64(self.information['CD2_1'])        , 'partial of second axis coordinate w.r.t. x')
+        hdu_g.header['CD2_2']  =( np.float64(self.information['CD2_2'])             , 'partial of second axis coordinate w.r.t. y')
+        hdu_g.header['others'] =( 0.0       ,         'other specific parameters')
+        
+        ######  Readout information  #######
+        hdu_g.header['GAIN1'] =( np.float32(self.information['g_gain1'])       , 'CCD gain (channel 1)')
+        hdu_g.header['GAIN2'] =( np.float32(self.information['g_gain2'])       , 'CCD gain (channel 2)')
+        hdu_g.header['GAIN3'] =( np.float32(self.information['g_gain3'])       , 'CCD gain (channel 3)')
+        hdu_g.header['GAIN4'] =( np.float32(self.information['g_gain4'])       , 'CCD gain (channel 4)')
+        hdu_g.header['GAIN5'] =( np.float32(self.information['g_gain5'])       , 'CCD gain (channel 5)')
+        hdu_g.header['GAIN6'] =( np.float32(self.information['g_gain6'])       , 'CCD gain (channel 6)')
+        hdu_g.header['GAIN7'] =( np.float32(self.information['g_gain7'])       , 'CCD gain (channel 7)')
+        hdu_g.header['GAIN8'] =( np.float32(self.information['g_gain8'])       , 'CCD gain (channel 8)')
+        hdu_g.header['GAIN9'] =( np.float32(self.information['g_gain9'])       , 'CCD gain (channel 9)')
+        hdu_g.header['GAIN10']=( np.float32(self.information['g_gain10'])       , 'CCD gain (channel 10)')
+        hdu_g.header['GAIN11']=( np.float32(self.information['g_gain11'])       , 'CCD gain (channel 11)')
+        hdu_g.header['GAIN12']=( np.float32(self.information['g_gain12'])       , 'CCD gain (channel 12)')
+        hdu_g.header['GAIN13']=( np.float32(self.information['g_gain13'])       , 'CCD gain (channel 13)')
+        hdu_g.header['GAIN14']=( np.float32(self.information['g_gain14'])       , 'CCD gain (channel 14)')
+        hdu_g.header['GAIN15']=( np.float32(self.information['g_gain15'])       , 'CCD gain (channel 15)')
+        hdu_g.header['GAIN16']=( np.float32(self.information['g_gain16'])       , 'CCD gain (channel 16)')
+        #######
+        hdu_g.header['RDNOIS1'] =( np.float32(self.information['g_rdnois1'])       , 'read noise (channel 1)')
+        hdu_g.header['RDNOIS2'] =( np.float32(self.information['g_rdnois2'])       , 'read noise (channel 2)')
+        hdu_g.header['RDNOIS3'] =( np.float32(self.information['g_rdnois3'])       , 'read noise (channel 3)')
+        hdu_g.header['RDNOIS4'] =( np.float32(self.information['g_rdnois4'])       , 'read noise (channel 4)')
+        hdu_g.header['RDNOIS5'] =( np.float32(self.information['g_rdnois5'])       , 'read noise (channel 5)')
+        hdu_g.header['RDNOIS6'] =( np.float32(self.information['g_rdnois6'])       , 'read noise (channel 6)')
+        hdu_g.header['RDNOIS7'] =( np.float32(self.information['g_rdnois7'])       , 'read noise (channel 7)')
+        hdu_g.header['RDNOIS8'] =( np.float32(self.information['g_rdnois8'])       , 'read noise (channel 8)')
+        hdu_g.header['RDNOIS9'] =( np.float32(self.information['g_rdnois9'])       , 'read noise (channel 9)')
+        hdu_g.header['RDNOIS10']=( np.float32(self.information['g_rdnois10'])       , 'read noise (channel 10)')
+        hdu_g.header['RDNOIS11']=( np.float32(self.information['g_rdnois11'])       , 'read noise (channel 11)')
+        hdu_g.header['RDNOIS12']=( np.float32(self.information['g_rdnois12'])       , 'read noise (channel 12)')
+        hdu_g.header['RDNOIS13']=( np.float32(self.information['g_rdnois13'])       , 'read noise (channel 13)')
+        hdu_g.header['RDNOIS14']=( np.float32(self.information['g_rdnois14'])       , 'read noise (channel 14)')
+        hdu_g.header['RDNOIS15']=( np.float32(self.information['g_rdnois15'])       , 'read noise (channel 15)')
+        hdu_g.header['RDNOIS16']=( np.float32(self.information['g_rdnois16'])       , 'read noise (channel 16)')
+        
+        #######
+        hdu_g.header['DETBIA1']=( np.float32(self.information['g_detbia1'])      , 'amplifier bias voltage (channel 1)')
+        hdu_g.header['DETBIA2']=( np.float32(self.information['g_detbia2'])      , 'amplifier bias voltage (channel 2)')
+        hdu_g.header['DETBIA3']=( np.float32(self.information['g_detbia3'])      , 'amplifier bias voltage (channel 3)')
+        hdu_g.header['DETBIA4']=( np.float32(self.information['g_detbia4'])      , 'amplifier bias voltage (channel 4)')
+        hdu_g.header['DETBIA5']=( np.float32(self.information['g_detbia5'])      , 'amplifier bias voltage (channel 5)')
+        hdu_g.header['DETBIA6']=( np.float32(self.information['g_detbia6'])      , 'amplifier bias voltage (channel 6)')
+        hdu_g.header['DETBIA7']=( np.float32(self.information['g_detbia7'])      , 'amplifier bias voltage (channel 7)')
+        hdu_g.header['DETBIA8']=( np.float32(self.information['g_detbia8'])      , 'amplifier bias voltage (channel 8)')
+        hdu_g.header['DETBIA9']=( np.float32(self.information['g_detbia9'])      , 'amplifier bias voltage (channel 9)')
+        hdu_g.header['DETBIA10']=( np.float32(self.information['g_detbia10'])      , 'amplifier bias voltage (channel 10)')
+        hdu_g.header['DETBIA11']=( np.float32(self.information['g_detbia11'])      , 'amplifier bias voltage (channel 11)')
+        hdu_g.header['DETBIA12']=( np.float32(self.information['g_detbia12'])      , 'amplifier bias voltage (channel 12)')
+        hdu_g.header['DETBIA13']=( np.float32(self.information['g_detbia13'])      , 'amplifier bias voltage (channel 13)')
+        hdu_g.header['DETBIA14']=( np.float32(self.information['g_detbia14'])      , 'amplifier bias voltage (channel 14)')
+        hdu_g.header['DETBIA15']=( np.float32(self.information['g_detbia15'])      , 'amplifier bias voltage (channel 15)')
+        hdu_g.header['DETBIA16']=( np.float32(self.information['g_detbia16'])      , 'amplifier bias voltage (channel 16)')
+        
+        ######
+  
+        hdu_g.header['READT0']  =( time2jd(t2)       , 'read start time (UTC)')
+        hdu_g.header['READT1']  =( time2jd(t3)         , 'read end time (UTC)')
+        hdu_g.header['DETTEMP0']=( np.float32(0)       , 'detector temperature at READT0')
+        hdu_g.header['DETTEMP1']=( np.float32(0)       , 'detector temperature at READT1')
+        hdu_g.header['RDSPEED'] =( np.float32(3)       , 'read speed (in MHz)')
+        
+        #####  Shutter information #####
+        hdu_g.header['SHUTHWV'] =( 'XXXXXXXXXXXXXXXX'  , 'shutter hardware version')
+        hdu_g.header['SHUTSWV'] =( 'XXXXXXXXXXXXXXXX'  , 'shutter software')
+        hdu_g.header['SHUTSTAT']=( bool(True)          , 'shutter status,T: open, F: close')
+        hdu_g.header['SHTOPEN0']=( np.float64(0)           , 'shutter open time (begin)')
+        hdu_g.header['SHTOPEN1']=( np.float64(0)           , 'shutter open time (end)')
+        hdu_g.header['SHTCLOS0']=( np.float64(0)           , 'shutter close time (begin)')
+        hdu_g.header['SHTCLOS1']=( np.float64(0)           , 'shutter close time (end)')
+        ################ finish MCI header for image layer
+        
+
+        ##############################################################################################333   
+        
+        hdulist_g=fits.HDUList([ofd_g, hdu_g])
+        
+        if self.source=='DARK' or self.source=='FLAT' or self.source=='BIAS':            
+            file_g=self.result_path+'/cali_Data/'+filename_g
+        else:
+            file_g=self.result_path+'/sky_Data/'+filename_g 
+            
+        hdulist_g.writeto(file_g, overwrite=True)
+
+        ###########################################################################
+        ### r band   
+    
+        #create a new FITS file, using HDUList instance
+        ofd_r = fits.PrimaryHDU()
+        
+        filename_r='CSST_MCI_C2_'+self.source+'_'+exp_starttime+'_'+exp_endtime+'_'+str(obsid)+'_'+ \
+                   str(self.filterP[self.filter_r]['number'])+'_L0_'+self.information['ver']+'.fits'
+        
+        #########    add primary header   ##########
+        ofd_r.header['GROUPS']=( bool(False),   'always F')
+        ofd_r.header['DATE']  =( ofd_g.header['DATE'],   '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']=( 'MCI',   '   ')
+        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  #############             
+        
+        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']=(ofd_g.header['DATE-OBS'] ,   'date of the observation (yyyy-mm-dd hh:mm:ss)')
+        
+        ofd_r.header['EXPSTART']=(ofd_g.header['EXPSTART'],   '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['pos_x']) ,   'the orbital position of CSST in X direction at EXPSTART')
+        ofd_r.header['POSI0_Y'] =(np.float64(self.information['pos_y']) ,   'the orbital position of CSST in Y direction at EXPSTART')
+        ofd_r.header['POSI0_Z'] =(np.float64(self.information['pos_z']) ,   'the orbital position of CSST in Z direction at EXPSTART')
+        ofd_r.header['VELO0_X']=( np.float64(self.information['velocity_x']) ,   'the orbital velocity of CSST in X direction at EXPSTART')
+        ofd_r.header['VELO0_Y']=( np.float64(self.information['velocity_y']) ,   'the orbital velocity of CSST in Y direction at EXPSTART')
+        ofd_r.header['VELO0_Z']=( np.float64(self.information['velocity_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']  =(ofd_g.header['EXPEND'] ,       '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(p1x) ,   'the orbital position of CSST in X direction at EXPEND')
+        ofd_r.header['POSI1_Y'] =(np.float64(p1y) ,   'the orbital position of CSST in Y direction at EXPEND')
+        ofd_r.header['POSI1_Z'] =(np.float64(p1z) ,   'the orbital position of CSST in Z direction at EXPEND')
+        
+        ofd_r.header['VELO1_X']=(np.float64(p1vx) ,   'the orbital velocity of CSST in X direction at EXPEND')
+        ofd_r.header['VELO1_Y']=(np.float64(p1vy) ,   'the orbital velocity of CSST in Y direction at EXPEND')
+        ofd_r.header['VELO1_Z']=(np.float64(p1vz) ,   '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(0.0),    'RA of the pointing (degrees) at EXPEND in deg')
+        ofd_r.header['DEC_PNT1']=(np.float64(0.0),   '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 header for 0 layer  
+        #new image HDU
+        hdu_r = fits.ImageHDU(data=np.uint16(self.image_r))
+
+        
+        ###  MCI  header for immage layer ###
+        #### Instrument information (at EXPSTART time)
+        hdu_r.header['PMIRRPOS']=( bool(True) ,  'FSM pointing,T: point to MCI, F: not point to MCI   ')
+        hdu_r.header['DIFFUSER']=( bool(True) ,  'insert diffuser for flat calibration,T: insert diffuser, F: not insert diffuser   ')
+        hdu_r.header['FLAMP']   =( np.int16(0)  ,   'status of flat lamp,0: off, 1: 1, 2: 2, 3: 3')
+        hdu_r.header['MCITEMP'] =( np.float32(0.0), 'MCI components temperature')
+        hdu_r.header['MCISTAT'] =( np.int16(0)  ,   'MCI components status parameters')
+        
+        #####  Filter information  #####
+        hdu_r.header['FILTER'] =(self.filter_r , 'filter band name')
+        
+        #####  Detector information ####
+        hdu_r.header['DETNAM'] =('C2'     , 'detector name for channle 2')
+        hdu_r.header['DETTEMP'] =(np.float32(0)     , 'detector temperature')
+        hdu_r.header['DETSIZE'] =(str(self.image_r.shape[0])+'*'+str(self.image_r.shape[1]), 'detector size')
+        hdu_r.header['DATASEC'] =(str(self.information['ysize'])+'*'+str( self.information['xsize']), 'data section')
+        hdu_r.header['PIXSCAL1']=(np.float32(0.05)     , 'pixel scale for axis 1')
+        hdu_r.header['PIXSCAL2']=(np.float32(0.05)     , 'pixel scale for axis 2')
+        hdu_r.header['PIXSIZE1']=(np.float32(10.0)     , 'pixel size for axis 1')
+        hdu_r.header['PIXSIZE2']=(np.float32(10.0)     , 'pixel size for axis 2')
+        hdu_r.header['NCHAN']   =(np.int16(16)     , 'number of readout channels')
+        hdu_r.header['NCHAN1']  =(np.int16(8)       , 'number of horizontal channels')
+        hdu_r.header['NCHAN2']  =(np.int16(2)       , 'number of vertical channels')
+        hdu_r.header['PSCAN1']  =(np.int32(27)      , 'horizontal prescan width, per readout channel')
+        hdu_r.header['PSCAN2']  =(np.int32(8)       , 'vertical prescan height, per readout channel')
+        hdu_r.header['OSCAN1']  =(np.int32(16)      , 'horizontal overscan width, per readout channel')
+        hdu_r.header['OSCAN2']  =(np.int32(8)       , 'vertical overscan width, per readout channel')
+        hdu_r.header['BIN_X']   =(np.int16(0)       , 'bin number in X (wavelength)')
+        hdu_r.header['BIN_Y']   =(np.int16(0)       , 'bin number in Y (spatial)')
+        
+        ####   World coordinate system and related parameters  #####
+        
+        hdu_r.header['WCSAXES']=( np.int16(2)        ,    'number of World Coordinate System axes')
+        hdu_r.header['CRPIX1'] =( np.float64(self.information['CRPIX1'])       , 'x-coordinate of reference pixel')
+        hdu_r.header['CRPIX2'] =( np.float64(self.information['CRPIX2'])       , 'y-coordinate of reference pixel')
+        hdu_r.header['CRVAL1'] =( np.float64(self.information['CRVAL1'])       , 'first axis value at reference pixel')
+        hdu_r.header['CRVAL2'] =( np.float64(self.information['CRVAL2'])       , 'second axis value at reference pixel')
+        hdu_r.header['CTYPE1'] =( 'RA---TAN'                               , 'the coordinate type for the first axis')
+        hdu_r.header['CTYPE2'] =( 'DEC---TAN'                              , 'the coordinate type for the second axis')
+        hdu_r.header['CD1_1']  =( np.float64(self.information['CD1_1'])        , 'partial of first axis coordinate w.r.t. x')
+        hdu_r.header['CD1_2']  =( np.float64(self.information['CD1_2'])        , 'partial of first axis coordinate w.r.t. y')
+        hdu_r.header['CD2_1']  =( np.float64(self.information['CD2_1'])        , 'partial of second axis coordinate w.r.t. x')
+        hdu_r.header['CD2_2']  =( np.float64(self.information['CD2_2'])             , 'partial of second axis coordinate w.r.t. y')
+        hdu_r.header['others'] =( 0.0       ,         'other specific parameters')
+        
+        ######  Readout information  #######
+        hdu_r.header['GAIN1'] =( np.float32(self.information['r_gain1'])       , 'CCD gain (channel 1)')
+        hdu_r.header['GAIN2'] =( np.float32(self.information['r_gain2'])       , 'CCD gain (channel 2)')
+        hdu_r.header['GAIN3'] =( np.float32(self.information['r_gain3'])       , 'CCD gain (channel 3)')
+        hdu_r.header['GAIN4'] =( np.float32(self.information['r_gain4'])       , 'CCD gain (channel 4)')
+        hdu_r.header['GAIN5'] =( np.float32(self.information['r_gain5'])       , 'CCD gain (channel 5)')
+        hdu_r.header['GAIN6'] =( np.float32(self.information['r_gain6'])       , 'CCD gain (channel 6)')
+        hdu_r.header['GAIN7'] =( np.float32(self.information['r_gain7'])       , 'CCD gain (channel 7)')
+        hdu_r.header['GAIN8'] =( np.float32(self.information['r_gain8'])       , 'CCD gain (channel 8)')
+        hdu_r.header['GAIN9'] =( np.float32(self.information['r_gain9'])       , 'CCD gain (channel 9)')
+        hdu_r.header['GAIN10']=( np.float32(self.information['r_gain10'])       , 'CCD gain (channel 10)')
+        hdu_r.header['GAIN11']=( np.float32(self.information['r_gain11'])       , 'CCD gain (channel 11)')
+        hdu_r.header['GAIN12']=( np.float32(self.information['r_gain12'])       , 'CCD gain (channel 12)')
+        hdu_r.header['GAIN13']=( np.float32(self.information['r_gain13'])       , 'CCD gain (channel 13)')
+        hdu_r.header['GAIN14']=( np.float32(self.information['r_gain14'])       , 'CCD gain (channel 14)')
+        hdu_r.header['GAIN15']=( np.float32(self.information['r_gain15'])       , 'CCD gain (channel 15)')
+        hdu_r.header['GAIN16']=( np.float32(self.information['r_gain16'])       , 'CCD gain (channel 16)')
+        #######
+        hdu_r.header['RDNOIS1'] =( np.float32(self.information['r_rdnois1'])       , 'read noise (channel 1)')
+        hdu_r.header['RDNOIS2'] =( np.float32(self.information['r_rdnois2'])       , 'read noise (channel 2)')
+        hdu_r.header['RDNOIS3'] =( np.float32(self.information['r_rdnois3'])       , 'read noise (channel 3)')
+        hdu_r.header['RDNOIS4'] =( np.float32(self.information['r_rdnois4'])       , 'read noise (channel 4)')
+        hdu_r.header['RDNOIS5'] =( np.float32(self.information['r_rdnois5'])       , 'read noise (channel 5)')
+        hdu_r.header['RDNOIS6'] =( np.float32(self.information['r_rdnois6'])       , 'read noise (channel 6)')
+        hdu_r.header['RDNOIS7'] =( np.float32(self.information['r_rdnois7'])       , 'read noise (channel 7)')
+        hdu_r.header['RDNOIS8'] =( np.float32(self.information['r_rdnois8'])       , 'read noise (channel 8)')
+        hdu_r.header['RDNOIS9'] =( np.float32(self.information['r_rdnois9'])       , 'read noise (channel 9)')
+        hdu_r.header['RDNOIS10']=( np.float32(self.information['r_rdnois10'])       , 'read noise (channel 10)')
+        hdu_r.header['RDNOIS11']=( np.float32(self.information['r_rdnois11'])       , 'read noise (channel 11)')
+        hdu_r.header['RDNOIS12']=( np.float32(self.information['r_rdnois12'])       , 'read noise (channel 12)')
+        hdu_r.header['RDNOIS13']=( np.float32(self.information['r_rdnois13'])       , 'read noise (channel 13)')
+        hdu_r.header['RDNOIS14']=( np.float32(self.information['r_rdnois14'])       , 'read noise (channel 14)')
+        hdu_r.header['RDNOIS15']=( np.float32(self.information['r_rdnois15'])       , 'read noise (channel 15)')
+        hdu_r.header['RDNOIS16']=( np.float32(self.information['r_rdnois16'])       , 'read noise (channel 16)')
+        
+        #######
+        hdu_r.header['DETBIA1']=( np.float32(self.information['r_detbia1'])      , 'amplifier bias voltage (channel 1)')
+        hdu_r.header['DETBIA2']=( np.float32(self.information['r_detbia2'])      , 'amplifier bias voltage (channel 2)')
+        hdu_r.header['DETBIA3']=( np.float32(self.information['r_detbia3'])      , 'amplifier bias voltage (channel 3)')
+        hdu_r.header['DETBIA4']=( np.float32(self.information['r_detbia4'])      , 'amplifier bias voltage (channel 4)')
+        hdu_r.header['DETBIA5']=( np.float32(self.information['r_detbia5'])      , 'amplifier bias voltage (channel 5)')
+        hdu_r.header['DETBIA6']=( np.float32(self.information['r_detbia6'])      , 'amplifier bias voltage (channel 6)')
+        hdu_r.header['DETBIA7']=( np.float32(self.information['r_detbia7'])      , 'amplifier bias voltage (channel 7)')
+        hdu_r.header['DETBIA8']=( np.float32(self.information['r_detbia8'])      , 'amplifier bias voltage (channel 8)')
+        hdu_r.header['DETBIA9']=( np.float32(self.information['r_detbia9'])      , 'amplifier bias voltage (channel 9)')
+        hdu_r.header['DETBIA10']=( np.float32(self.information['r_detbia10'])      , 'amplifier bias voltage (channel 10)')
+        hdu_r.header['DETBIA11']=( np.float32(self.information['r_detbia11'])      , 'amplifier bias voltage (channel 11)')
+        hdu_r.header['DETBIA12']=( np.float32(self.information['r_detbia12'])      , 'amplifier bias voltage (channel 12)')
+        hdu_r.header['DETBIA13']=( np.float32(self.information['r_detbia13'])      , 'amplifier bias voltage (channel 13)')
+        hdu_r.header['DETBIA14']=( np.float32(self.information['r_detbia14'])      , 'amplifier bias voltage (channel 14)')
+        hdu_r.header['DETBIA15']=( np.float32(self.information['r_detbia15'])      , 'amplifier bias voltage (channel 15)')
+        hdu_r.header['DETBIA16']=( np.float32(self.information['r_detbia16'])      , 'amplifier bias voltage (channel 16)')
+        
+        ######
+        ######
+        hdu_r.header['READT0']  =( hdu_g.header['READT0']       , 'read start time (UTC)')
+        hdu_r.header['READT1']  =( hdu_g.header['READT1']       , 'read end time (UTC)')
+        hdu_r.header['DETTEMP0']=( np.float32(0)       , 'detector temperature at READT0')
+        hdu_r.header['DETTEMP1']=( np.float32(0)       , 'detector temperature at READT1')
+        hdu_r.header['RDSPEED'] =( np.float32(3)       , 'read speed (in MHz)')
+        
+        #####  Shutter information #####
+        hdu_r.header['SHUTHWV'] =( 'XXXXXXXXXXXXXXXX'  , 'shutter hardware version')
+        hdu_r.header['SHUTSWV'] =( 'XXXXXXXXXXXXXXXX'  , 'shutter software')
+        hdu_r.header['SHUTSTAT']=( bool(True)          , 'shutter status,T: open, F: close')
+        hdu_r.header['SHTOPEN0']=( np.float64(0)           , 'shutter open time (begin)')
+        hdu_r.header['SHTOPEN1']=( np.float64(0)           , 'shutter open time (end)')
+        hdu_r.header['SHTCLOS0']=( np.float64(0)           , 'shutter close time (begin)')
+        hdu_r.header['SHTCLOS1']=( np.float64(0)           , 'shutter close time (end)')
+        ################ finish MCI header for image layer
+        
+        #write the actual file
+        hdulist_r=fits.HDUList([ofd_r, hdu_r])
+        
+        #file_r=self.result_path+'/sky_Data/'+filename_r
+        
+        if self.source=='DARK' or self.source=='FLAT' or self.source=='BIAS':            
+            file_r=self.result_path+'/cali_Data/'+filename_r
+        else:
+            file_r=self.result_path+'/sky_Data/'+filename_r 
+            
+        hdulist_r.writeto(file_r, overwrite=True)
+        
+        
+        #################################################################################
+        ### i band
+        
+        #create a new FITS file, using HDUList instance
+        ofd_i = fits.PrimaryHDU()
+        
+        filename_i='CSST_MCI_C3_'+self.source+'_'+exp_starttime+'_'+exp_endtime+'_'+str(obsid)+'_'+ \
+                   str(self.filterP[self.filter_i]['number'])+'_L0_'+self.information['ver']+'.fits'
+        
+        #### add primary header
+        ofd_i.header['GROUPS']=( bool(False),   'always F')
+        ofd_i.header['DATE']  =( ofd_g.header['DATE'],   'date this file was written' ) 
+        
+        ofd_i.header['FILENAME']=(filename_i,   ' file name C48  ')        
+        ofd_i.header['OBSTYPE'] =( self.source,   'observation type raw,flt, mask, bias, dark, sci')        
+        ofd_i.header['TELESCOP']=('CSST',   'always CSST')
+        ofd_i.header['INSTRUME']=( 'MCI',   '   ')
+        ofd_i.header['RADECSYS']=('ICRS',   ' always ICRS ')
+        ofd_i.header['EQUINOX'] =( float(2000.0),   'always 2000.0')
+        ofd_i.header['FITSCREA']=( '4.2.1',   'FITS create software version')
+        
+        #########   Object information  #############             
+        
+        ofd_i.header['OBJECT']=( self.information['name_obj'],   'object name')
+        ofd_i.header['TARGET']=( self.information['target'],     'target name, hhmmss+ddmmss')
+        ofd_i.header['OBJ_RA'] =(np.float64(self.information['ra_obj']) ,   'RA of the object in deg')
+        ofd_i.header['OBJ_DEC']=(np.float64(self.information['dec_obj']) ,   'DEC of the object in deg')
+        
+        ofd_i.header['RA_PNT0'] =( np.float64(self.information['ra_pnt0']) ,     'RA of the pointing (degrees) at EXPSTART')
+        ofd_i.header['DEC_PNT0']=(np.float64(self.information['dec_pnt0']) ,    'DEC of the pointing (degrees) at EXPSTART')
+    
+
+        
+        ##############
+        ofd_i.header['OBSID'] =(str(obsid) ,   'observation ID, 3+8bit')
+        
+        ########  Telescope information ###############
+        # ofd_i.header['COMMENT'] ='=========================================================='
+        # ofd_i.header['COMMENT'] ='Telescope information'
+        # ofd_i.header['COMMENT'] ='=========================================================='       
+        
+        ofd_i.header['REFFRAME']=('CSSTGSC-1.0' ,   'guide star catalog version')
+        ofd_i.header['DATE-OBS']=(ofd_g.header['DATE-OBS'] ,   'date of the observation (yyyy-mm-dd hh:mm:ss)')
+        
+        ofd_i.header['EXPSTART']=(ofd_g.header['EXPSTART'],   'exposure start time')
+        ofd_i.header['SUNANGL0']=(np.float32(0.0) ,   'angle between sun and optical axis at EXPSTART')
+        ofd_i.header['MOONANG0']=(np.float32(0.0) ,   'angle between moon and optical axis at EXPSTART')
+        ofd_i.header['POS_ANG0']=(np.float64(0.0),   'angle between optical axis and the North Pole at EXPSTART in arcsec')
+        ofd_i.header['TEL_ALT0']=(np.float64(0.0),   'angle between optical axis and the ground- piston at EXPSTART in deg')
+        
+        ofd_i.header['HOODSTA0']=(np.float32(0.0) ,   'lens hood altitude at EXPSTART')
+        ofd_i.header['HOODANG0']=(np.float32(0.0),    'lens hood azimuth at EXPSTART')
+        ofd_i.header['POSI0_X'] =(np.float64(self.information['pos_x']) ,   'the orbital position of CSST in X direction at EXPSTART')
+        ofd_i.header['POSI0_Y'] =(np.float64(self.information['pos_y']) ,   'the orbital position of CSST in Y direction at EXPSTART')
+        ofd_i.header['POSI0_Z'] =(np.float64(self.information['pos_z']) ,   'the orbital position of CSST in Z direction at EXPSTART')
+        ofd_i.header['VELO0_X']=( np.float64(self.information['velocity_x']) ,   'the orbital velocity of CSST in X direction at EXPSTART')
+        ofd_i.header['VELO0_Y']=( np.float64(self.information['velocity_y']) ,   'the orbital velocity of CSST in Y direction at EXPSTART')
+        ofd_i.header['VELO0_Z']=( np.float64(self.information['velocity_z']) ,   'the orbital velocity of CSST in Z direction at EXPSTART')
+        
+        ofd_i.header['Euler0_1']=( np.float64(0.0),   'Euler angle 1 at EXPSTART')
+        ofd_i.header['Euler0_2']=( np.float64(0.0),   'Euler angle 2 at EXPSTART')
+        ofd_i.header['Euler0_3']=( np.float64(0.0),   'Euler angle 3 at EXPSTART')        
+        
+        ofd_i.header['EXPEND']  =(ofd_g.header['EXPEND'] ,       'exposure end time')
+        
+        ofd_i.header['SUNANGL1']=(np.float32(0.0),   'angle between sun and optical axis at EXPEND')
+        ofd_i.header['MOONANG1']=(np.float32(0.0) ,   'angle between moon and optical axis at EXPEND   ')
+        ofd_i.header['POS_ANG1']=(np.float64(0.0) ,   'angle between optical axis and the North Pole at EXPEND in arcsec')
+        ofd_i.header['TEL_ALT1']=(np.float64(0.0) ,   'angle between optical axis and the ground- piston at EXPEND in deg  ')
+        ofd_i.header['HOODSTA1']=(np.float32(0.0),   'lens hood altitude at EXPEND   ')
+        ofd_i.header['HOODANG1']=(np.float32(0.0),   'lens hood azimuth at EXPEND   ')
+        ofd_i.header['POSI1_X'] =(np.float64(p1x) ,   'the orbital position of CSST in X direction at EXPEND')
+        ofd_i.header['POSI1_Y'] =(np.float64(p1y) ,   'the orbital position of CSST in Y direction at EXPEND')
+        ofd_i.header['POSI1_Z'] =(np.float64(p1z) ,   'the orbital position of CSST in Z direction at EXPEND')
+        
+        ofd_i.header['VELO1_X']=(np.float64(p1vx) ,   'the orbital velocity of CSST in X direction at EXPEND')
+        ofd_i.header['VELO1_Y']=(np.float64(p1vy) ,   'the orbital velocity of CSST in Y direction at EXPEND')
+        ofd_i.header['VELO1_Z']=(np.float64(p1vz) ,   'the orbital velocity of CSST in Z direction at EXPEND')
+        
+        ofd_i.header['Euler1_1']=( np.float64(0.0),   'Euler angle 1 at EXPEND')
+        ofd_i.header['Euler1_2']=( np.float64(0.0),   'Euler angle 2 at EXPEND')
+        ofd_i.header['Euler1_3']=( np.float64(0.0),   'Euler angle 3 at EXPEND')
+        
+        ofd_i.header['RA_PNT1']=(np.float64(0.0),    'RA of the pointing (degrees) at EXPEND in deg')
+        ofd_i.header['DEC_PNT1']=(np.float64(0.0),   'DEC of the pointing (degrees) at EXPEND in deg')
+        
+        ofd_i.header['EXPTIME']=(self.information['exptime'],   'exposure duration')
+        ofd_i.header['EPOCH']  =(np.float32(0.0),   'coordinate epoch')
+        ofd_i.header['CHECKSUM']=(   0           ,   'hdu-checksum')  
+        ########## finish header for 0 layer  
+
+        #new image HDU
+        hdu_i = fits.ImageHDU(data=np.uint16(self.image_i))
+        
+        
+        
+        ###  MCI  header for immage layer ###
+        #### Instrument information (at EXPSTART time)
+        hdu_i.header['PMIRRPOS']=( bool(True) ,  'FSM pointing,T: point to MCI, F: not point to MCI   ')
+        hdu_i.header['DIFFUSER']=( bool(True) ,  'insert diffuser for flat calibration,T: insert diffuser, F: not insert diffuser   ')
+        hdu_i.header['FLAMP']   =( np.int16(0)  ,   'status of flat lamp,0: off, 1: 1, 2: 2, 3: 3')
+        hdu_i.header['MCITEMP'] =( np.float32(0.0), 'MCI components temperature')
+        hdu_i.header['MCISTAT'] =( np.int16(0)  ,   'MCI components status parameters')
+        
+        #####  Filter information  #####
+        hdu_i.header['FILTER'] =(self.filter_i , 'filter band name')
+       
+        #####  Detector information ####
+        hdu_i.header['DETNAM'] =('C3'     , 'detector name for channle 3')
+        hdu_i.header['DETTEMP'] =(np.float32(0)     , 'detector temperature')
+        hdu_i.header['DETSIZE'] =(str(self.image_i.shape[0])+'*'+str(self.image_i.shape[1]), 'detector size')
+        hdu_i.header['DATASEC'] =(str(self.information['ysize'])+'*'+str( self.information['xsize']), 'data section')
+        hdu_i.header['PIXSCAL1']=(np.float32(0.05)     , 'pixel scale for axis 1')
+        hdu_i.header['PIXSCAL2']=(np.float32(0.05)     , 'pixel scale for axis 2')
+        hdu_i.header['PIXSIZE1']=(np.float32(10.0)     , 'pixel size for axis 1')
+        hdu_i.header['PIXSIZE2']=(np.float32(10.0)     , 'pixel size for axis 2')
+        hdu_i.header['NCHAN']   =(np.int16(16)     , 'number of readout channels')
+        hdu_i.header['NCHAN1']  =(np.int16(8)       , 'number of horizontal channels')
+        hdu_i.header['NCHAN2']  =(np.int16(2)       , 'number of vertical channels')
+        hdu_i.header['PSCAN1']  =(np.int32(27)      , 'horizontal prescan width, per readout channel')
+        hdu_i.header['PSCAN2']  =(np.int32(8)       , 'vertical prescan height, per readout channel')
+        hdu_i.header['OSCAN1']  =(np.int32(16)      , 'horizontal overscan width, per readout channel')
+        hdu_i.header['OSCAN2']  =(np.int32(8)       , 'vertical overscan width, per readout channel')
+        hdu_i.header['BIN_X']   =(np.int16(0)       , 'bin number in X (wavelength)')
+        hdu_i.header['BIN_Y']   =(np.int16(0)       , 'bin number in Y (spatial)')
+        
+        ####   World coordinate system and related parameters  #####
+        
+        hdu_i.header['WCSAXES']=( np.int16(2)        ,    'number of World Coordinate System axes')
+        hdu_i.header['CRPIX1'] =( np.float64(self.information['CRPIX1'])       , 'x-coordinate of reference pixel')
+        hdu_i.header['CRPIX2'] =( np.float64(self.information['CRPIX2'])       , 'y-coordinate of reference pixel')
+        hdu_i.header['CRVAL1'] =( np.float64(self.information['CRVAL1'])       , 'first axis value at reference pixel')
+        hdu_i.header['CRVAL2'] =( np.float64(self.information['CRVAL2'])       , 'second axis value at reference pixel')
+        hdu_i.header['CTYPE1'] =( 'RA---TAN'                               , 'the coordinate type for the first axis')
+        hdu_i.header['CTYPE2'] =( 'DEC---TAN'                              , 'the coordinate type for the second axis')
+        hdu_i.header['CD1_1']  =( np.float64(self.information['CD1_1'])        , 'partial of first axis coordinate w.r.t. x')
+        hdu_i.header['CD1_2']  =( np.float64(self.information['CD1_2'])        , 'partial of first axis coordinate w.r.t. y')
+        hdu_i.header['CD2_1']  =( np.float64(self.information['CD2_1'])        , 'partial of second axis coordinate w.r.t. x')
+        hdu_i.header['CD2_2']  =( np.float64(self.information['CD2_2'])             , 'partial of second axis coordinate w.r.t. y')
+        hdu_i.header['others'] =( 0.0       ,         'other specific parameters')
+        
+        ######  Readout information  #######
+        hdu_i.header['GAIN1'] =( np.float32(self.information['i_gain1'])       , 'CCD gain (channel 1)')
+        hdu_i.header['GAIN2'] =( np.float32(self.information['i_gain2'])       , 'CCD gain (channel 2)')
+        hdu_i.header['GAIN3'] =( np.float32(self.information['i_gain3'])       , 'CCD gain (channel 3)')
+        hdu_i.header['GAIN4'] =( np.float32(self.information['i_gain4'])       , 'CCD gain (channel 4)')
+        hdu_i.header['GAIN5'] =( np.float32(self.information['i_gain5'])       , 'CCD gain (channel 5)')
+        hdu_i.header['GAIN6'] =( np.float32(self.information['i_gain6'])       , 'CCD gain (channel 6)')
+        hdu_i.header['GAIN7'] =( np.float32(self.information['i_gain7'])       , 'CCD gain (channel 7)')
+        hdu_i.header['GAIN8'] =( np.float32(self.information['i_gain8'])       , 'CCD gain (channel 8)')
+        hdu_i.header['GAIN9'] =( np.float32(self.information['i_gain9'])       , 'CCD gain (channel 9)')
+        hdu_i.header['GAIN10']=( np.float32(self.information['i_gain10'])       , 'CCD gain (channel 10)')
+        hdu_i.header['GAIN11']=( np.float32(self.information['i_gain11'])       , 'CCD gain (channel 11)')
+        hdu_i.header['GAIN12']=( np.float32(self.information['i_gain12'])       , 'CCD gain (channel 12)')
+        hdu_i.header['GAIN13']=( np.float32(self.information['i_gain13'])       , 'CCD gain (channel 13)')
+        hdu_i.header['GAIN14']=( np.float32(self.information['i_gain14'])       , 'CCD gain (channel 14)')
+        hdu_i.header['GAIN15']=( np.float32(self.information['i_gain15'])       , 'CCD gain (channel 15)')
+        hdu_i.header['GAIN16']=( np.float32(self.information['i_gain16'])       , 'CCD gain (channel 16)')
+        #######
+        hdu_i.header['RDNOIS1'] =( np.float32(self.information['i_rdnois1'])       , 'read noise (channel 1)')
+        hdu_i.header['RDNOIS2'] =( np.float32(self.information['i_rdnois2'])       , 'read noise (channel 2)')
+        hdu_i.header['RDNOIS3'] =( np.float32(self.information['i_rdnois3'])       , 'read noise (channel 3)')
+        hdu_i.header['RDNOIS4'] =( np.float32(self.information['i_rdnois4'])       , 'read noise (channel 4)')
+        hdu_i.header['RDNOIS5'] =( np.float32(self.information['i_rdnois5'])       , 'read noise (channel 5)')
+        hdu_i.header['RDNOIS6'] =( np.float32(self.information['i_rdnois6'])       , 'read noise (channel 6)')
+        hdu_i.header['RDNOIS7'] =( np.float32(self.information['i_rdnois7'])       , 'read noise (channel 7)')
+        hdu_i.header['RDNOIS8'] =( np.float32(self.information['i_rdnois8'])       , 'read noise (channel 8)')
+        hdu_i.header['RDNOIS9'] =( np.float32(self.information['i_rdnois9'])       , 'read noise (channel 9)')
+        hdu_i.header['RDNOIS10']=( np.float32(self.information['i_rdnois10'])       , 'read noise (channel 10)')
+        hdu_i.header['RDNOIS11']=( np.float32(self.information['i_rdnois11'])       , 'read noise (channel 11)')
+        hdu_i.header['RDNOIS12']=( np.float32(self.information['i_rdnois12'])       , 'read noise (channel 12)')
+        hdu_i.header['RDNOIS13']=( np.float32(self.information['i_rdnois13'])       , 'read noise (channel 13)')
+        hdu_i.header['RDNOIS14']=( np.float32(self.information['i_rdnois14'])       , 'read noise (channel 14)')
+        hdu_i.header['RDNOIS15']=( np.float32(self.information['i_rdnois15'])       , 'read noise (channel 15)')
+        hdu_i.header['RDNOIS16']=( np.float32(self.information['i_rdnois16'])       , 'read noise (channel 16)')
+        
+        #######
+        hdu_i.header['DETBIA1']=( np.float32(self.information['i_detbia1'])      , 'amplifier bias voltage (channel 1)')
+        hdu_i.header['DETBIA2']=( np.float32(self.information['i_detbia2'])      , 'amplifier bias voltage (channel 2)')
+        hdu_i.header['DETBIA3']=( np.float32(self.information['i_detbia3'])      , 'amplifier bias voltage (channel 3)')
+        hdu_i.header['DETBIA4']=( np.float32(self.information['i_detbia4'])      , 'amplifier bias voltage (channel 4)')
+        hdu_i.header['DETBIA5']=( np.float32(self.information['i_detbia5'])      , 'amplifier bias voltage (channel 5)')
+        hdu_i.header['DETBIA6']=( np.float32(self.information['i_detbia6'])      , 'amplifier bias voltage (channel 6)')
+        hdu_i.header['DETBIA7']=( np.float32(self.information['i_detbia7'])      , 'amplifier bias voltage (channel 7)')
+        hdu_i.header['DETBIA8']=( np.float32(self.information['i_detbia8'])      , 'amplifier bias voltage (channel 8)')
+        hdu_i.header['DETBIA9']=( np.float32(self.information['i_detbia9'])      , 'amplifier bias voltage (channel 9)')
+        hdu_i.header['DETBIA10']=( np.float32(self.information['i_detbia10'])      , 'amplifier bias voltage (channel 10)')
+        hdu_i.header['DETBIA11']=( np.float32(self.information['i_detbia11'])      , 'amplifier bias voltage (channel 11)')
+        hdu_i.header['DETBIA12']=( np.float32(self.information['i_detbia12'])      , 'amplifier bias voltage (channel 12)')
+        hdu_i.header['DETBIA13']=( np.float32(self.information['i_detbia13'])      , 'amplifier bias voltage (channel 13)')
+        hdu_i.header['DETBIA14']=( np.float32(self.information['i_detbia14'])      , 'amplifier bias voltage (channel 14)')
+        hdu_i.header['DETBIA15']=( np.float32(self.information['i_detbia15'])      , 'amplifier bias voltage (channel 15)')
+        hdu_i.header['DETBIA16']=( np.float32(self.information['i_detbia16'])      , 'amplifier bias voltage (channel 16)')
+        
+        ######
+        ######
+        hdu_i.header['READT0']  =( hdu_g.header['READT0']       , 'read start time (UTC)')
+        hdu_i.header['READT1']  =( hdu_g.header['READT1']       , 'read end time (UTC)')
+        hdu_i.header['DETTEMP0']=( np.float32(0)       , 'detector temperature at READT0')
+        hdu_i.header['DETTEMP1']=( np.float32(0)       , 'detector temperature at READT1')
+        hdu_i.header['RDSPEED'] =( np.float32(3)       , 'read speed (in MHz)')
+        
+        #####  Shutter information #####
+        hdu_i.header['SHUTHWV'] =( 'XXXXXXXXXXXXXXXX'  , 'shutter hardware version')
+        hdu_i.header['SHUTSWV'] =( 'XXXXXXXXXXXXXXXX'  , 'shutter software')
+        hdu_i.header['SHUTSTAT']=( bool(True)          , 'shutter status,T: open, F: close')
+        hdu_i.header['SHTOPEN0']=( np.float64(0)           , 'shutter open time (begin)')
+        hdu_i.header['SHTOPEN1']=( np.float64(0)           , 'shutter open time (end)')
+        hdu_i.header['SHTCLOS0']=( np.float64(0)           , 'shutter close time (begin)')
+        hdu_i.header['SHTCLOS1']=( np.float64(0)           , 'shutter close time (end)')
+        ################ finish MCI header for image layer
+
+        #write the actual file
+        hdulist_i=fits.HDUList([ofd_i, hdu_i])
+        
+        #file_i=self.result_path+'/sky_Data/'+filename_i
+        if self.source=='DARK' or self.source=='FLAT' or self.source=='BIAS':            
+            file_i=self.result_path+'/cali_Data/'+filename_i
+        else:
+            file_i=self.result_path+'/sky_Data/'+filename_i 
+        
+        hdulist_i.writeto(file_i, overwrite=True)
+        
+#################################################################################################        
+
+
+    def simulate(self, simnumber):
+        """
+        Create a single simulated image of a quadrant defined by the configuration file.
+        Will do all steps defined in the config file sequentially.
+
+        :return: None
+        """
+        
+        ##############################################################################################
+        ##############################################################################################
+        sourcelist=['XDF','CS','PI','TO']
+        
+        self.source=self.information['sourcein']
+        
+        self.information['ver']='1.0'   
+        
+        self._loadGhostModel()    
+         
+        self.load_filter_PSF()   
+        
+        self.information['simnumber']=simnumber
+        
+        
+        ############################################################################
+        flag=0
+        for k in range(1,50,1):         
+                        
+            fn=self.information['indata_path']+'/TianCe/orbit20160925/'+str(k)+'.txt';
+            d=np.loadtxt(fn);
+            self.dt_num=int(self.information['simnumber']*(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
+            
+        self.orbit_pars=d
+        self.orbit_file_num=k
+        self.orbit_exp_num =kk
+        
+        self.information['pos_x']=float(self.orbit_pars[self.orbit_exp_num,1])
+        self.information['pos_y']=float(self.orbit_pars[self.orbit_exp_num,2])
+        self.information['pos_z']=float(self.orbit_pars[self.orbit_exp_num,3])
+        self.information['velocity_x']=float(self.orbit_pars[self.orbit_exp_num,4])
+        self.information['velocity_y']=float(self.orbit_pars[self.orbit_exp_num,5])
+        self.information['velocity_z']=float(self.orbit_pars[self.orbit_exp_num,6])           
+        
+        exptime_start_jd=d[kk,0]  #### jd time, utc format                
+        
+        self.dt=julian.from_jd(exptime_start_jd, fmt='jd')
+        
+        #############################################################################
+        #############################################################################
+        
+
+        
+        if self.source in sourcelist:
+            
+            self.earthshine_theta=30.0       # in degree
+            
+
+
+            self.TianCe_day=str(self.dt.year)+self.dt.strftime("-%m-%d")  ####str(self.dt.year)+'-'+str(self.dt.month)+'-'+str(self.dt.day)
+            
+            self.TianCe_exp_start=dt2hmd(self.dt)
+            
+           ###################################
+            if self.sky_shift_rot:                
+                
+                np.random.seed()
+                ud = np.random.random()        # Choose a random pointing ra in degree
+                dis_ra  = 2*(ud-1)*10/3600.0
+                
+                np.random.seed()
+                ud = np.random.random()        # Choose a random pointing dec in degree
+                dis_dec = 2*(ud-1)*10/3600.0   #  
+                
+                np.random.seed()
+                ud = np.random.random()        # Choose a random telescope rotation in degree 
+                rotTelPos = 2*(ud-1)*5.0*galsim.degrees  # the image will rotate in clockwise rotaion between -10-10 degree
+                rotSkyPos = 0.0*galsim.degrees
+                
+                ################################
+                   
+            else:
+            
+                dis_ra  = 0     # 
+                dis_dec = 0     # T       
+                
+                rotTelPos = 0.0*galsim.degrees  # the image will rotate in clockwise rotaion between -10-10 degree
+                rotSkyPos = 0.0*galsim.degrees
+            
+            #####################################################
+            
+            self.information['T_disRa']  =dis_ra    # telescope pointing Ra  displacement  
+            self.information['T_disDec'] =dis_dec   # telescope pointing Dec displacement
+            
+            self.information['rotTelPos']=rotTelPos   # telescope rotation in degree
+            self.information['rotSkyPos']=rotSkyPos   # sky rotation in degree
+            
+            self.information['name_obj'] ='MCI_obj'     
+                            
+            theta = self.information['rotTelPos'] - self.information['rotSkyPos'] 
+            
+            self.information['rotangle']=theta.deg 
+           
+
+            self.log.info('dis_ra(in pixel)=%f, dis_dec(in pixel)=%f, sky_rot(in deg)=%f' % (dis_ra*3600/0.05, dis_dec*3600/0.05, theta.deg))
+           #######################################################################################################            
+   
+           #############################################################################
+           ##########   simulate star images from star SED data       ##################
+           
+            if  self.sim_star :
+                ############################
+        
+                starimg=self.cal_StarSED_img()
+                
+
+            ##########################################################################
+            # # ########################################################################
+            # ######   generate galaxy SED  image without Lensing effect    ########## 
+            
+            if self.sim_galaxy :
+                
+                losimg=self.cal_Lensing_galaxy_img()
+             
+
+            if self.sim_star ==False and self.sim_galaxy == False:
+                print('error: No image will be simulated...... ')
+                sys.exit(1)
+                
+            #### stack the image 
+               
+            if 'losimg' in dir():                    
+                self.image_g+=losimg['g']
+                self.image_r+=losimg['r']
+                self.image_i+=losimg['i']
+                
+            if 'starimg' in dir():                     
+                self.image_g+=starimg['g']
+                self.image_r+=starimg['r']
+                self.image_i+=starimg['i']     
+                
+                
+            ###########################################################################################################
+    
+            #########################################################################################################          
+            
+        ####### end if source =='SKY'    
+            
+        elif self.source=='FLAT':   
+            
+            theta=0
+            
+            self.information['name_obj']='FLAT'   
+                            
+            self.information['rotangle']=0    
+                
+            self.information['CRVAL1']=0
+            self.information['CRVAL2']=0            
+            
+            self.information['CRPIX1']=self.information['xsize']/2-0.5 #### 
+            self.information['CRPIX2']=self.information['ysize']/2-0.5 #### 
+            
+            self.information['tel_Pra']=0
+            self.information['tel_Pdec']=0
+                        
+            self.information['CD1_1']=-np.cos(theta)*self.information['pixel_size'] #### 
+            self.information['CD1_2']= np.sin(theta)*self.information['pixel_size'] #### 
+            
+            self.information['CD2_1']= np.sin(theta)*self.information['pixel_size'] #### 
+            self.information['CD2_2']= np.cos(theta)*self.information['pixel_size'] ####  
+            
+            self.information['ra_obj']  =0
+            self.information['dec_obj'] =0
+            self.information['ra_pnt0'] =0
+            self.information['dec_pnt0']=0
+            self.information['target']  =0                     
+                
+            self.image_g +=100*self.information['exptime']
+            self.image_r +=100*self.information['exptime']
+            self.image_i +=100*self.information['exptime']              
+                        
+        elif self.source=='DARK':  
+            
+            theta=0
+            
+                        
+            self.information['name_obj']='DARK'   
+                            
+            self.information['rotangle']=0    
+                
+            self.information['CRVAL1']=0
+            self.information['CRVAL2']=0
+            
+            self.information['tel_Pra']=0
+            self.information['tel_Pdec']=0
+            
+            self.information['CRPIX1']=self.information['xsize']/2-0.5 #### 
+            self.information['CRPIX2']=self.information['ysize']/2-0.5 #### 
+            
+            self.information['CD1_1']=-np.cos(theta)*self.information['pixel_size'] #### 
+            self.information['CD1_2']= np.sin(theta)*self.information['pixel_size'] #### 
+            
+            self.information['CD2_1']= np.sin(theta)*self.information['pixel_size'] #### 
+            self.information['CD2_2']= np.cos(theta)*self.information['pixel_size'] ####  
+            
+            self.information['ra_obj']  =0
+            self.information['dec_obj'] =0
+            self.information['ra_pnt0'] =0
+            self.information['dec_pnt0']=0
+            self.information['target']  =0                     
+                
+            self.image_g=np.zeros((self.information['ysize'], self.information['xsize']), dtype=np.float64)
+            self.image_r=np.zeros((self.information['ysize'], self.information['xsize']), dtype=np.float64)
+            self.image_i=np.zeros((self.information['ysize'], self.information['xsize']), dtype=np.float64)          
+             
+        elif self.source=='BIAS':
+            
+            theta=0
+    
+            self.information['exptime']=0.0
+  
+            self.information['name_obj']='BIAS'   
+                            
+            self.information['rotangle']=0    
+                
+            self.information['CRVAL1']=0
+            self.information['CRVAL2']=0
+            
+            self.information['tel_Pra']=0
+            self.information['tel_Pdec']=0
+            
+            self.information['CRPIX1']=self.information['xsize']/2-0.5 #### 
+            self.information['CRPIX2']=self.information['ysize']/2-0.5 #### 
+            
+            self.information['CD1_1']=-np.cos(theta)*self.information['pixel_size'] #### 
+            self.information['CD1_2']= np.sin(theta)*self.information['pixel_size'] #### 
+            
+            self.information['CD2_1']= np.sin(theta)*self.information['pixel_size'] #### 
+            self.information['CD2_2']= np.cos(theta)*self.information['pixel_size'] ####  
+            
+            self.information['ra_obj']  =0
+            self.information['dec_obj'] =0
+            self.information['ra_pnt0'] =0
+            self.information['dec_pnt0']=0
+            self.information['target']  =0  
+            
+            self.image_g=np.zeros((self.information['ysize'], self.information['xsize']), dtype=np.float64)
+            self.image_r=np.zeros((self.information['ysize'], self.information['xsize']), dtype=np.float64)
+            self.image_i=np.zeros((self.information['ysize'], self.information['xsize']), dtype=np.float64)
+            
+        else:
+            print('Souce is not correct and programe will return')  
+            sys.exit(1)
+            
+            
+   ########################################################################################################################
+
+   ########################################################################################################################         
+      
+        #Apply flat-field large scale structure for one chip
+        if self.flatfieldM:
+            self.applyFlat()
+            #print('applyFlat')
+        ########################################################
+        if self.PRNUeffect:
+            self.applyPRNUeffect()   
+            #print('applyPRNUeffect')
+
+        ####################################################################
+        if self.source in sourcelist:
+            if self.cosmicRays:
+                self.addCosmicRays()
+                #print('addCosmicRays finisth')
+                     
+        ##################################################
+        if self.skyback and self.source in sourcelist:
+            self.applyPoissonSkyNoise()
+            #print('apply Poisson Sky Noise')
+        #################################################    
+            
+        if self.darknoise:  
+            self.applyDarkCurrent()
+            #print('applyDarkCurrent finisth')
+               
+        #################################################
+
+        if self.shutterEffect:
+            self.addshuttereffect()
+            #print('add shutter effect ')
+        #################################################    
+            
+        if self.source not in sourcelist: 
+            self.applyPoissonNoise()
+         
+        ###################################################
+        if self.nonlinearity:
+            self.applyNonlinearity()
+            #print('applyNonlinearity')
+          
+        # ################################################
+        if self.bleeding:
+            self.image_g=self.applyBleeding(self.image_g.copy())
+            self.image_r=self.applyBleeding(self.image_r.copy())
+            self.image_i=self.applyBleeding(self.image_i.copy())
+            #print('apply bleeding effect finisth')
+     
+        ################################################
+
+        #### read CCD image to output matrix as defined         
+        
+        if self.overscans and self.information['xsize']==9216 and self.information['ysize']==9232 :           
+           
+            self.image_g=self.addScanEffect(self.image_g.copy())
+            self.image_r=self.addScanEffect(self.image_r.copy())
+            self.image_i=self.addScanEffect(self.image_i.copy())
+            #print('apply overscan zone 1-16')
+            
+        #################################################################
+        
+        if self.radiationDamage:
+            self.applyRadiationDamage()
+            #print('applyRadiationDamage()')
+        ######################################################################
+        ###### apply readoutNoise   ######
+    
+        if self.readoutNoise:
+            self.applyReadoutNoise()
+            #print('apply readout noise')
+                        
+        ######################################################################
+        if self.information['xsize']==9216 and self.information['ysize']==9232 :
+            self.electrons2ADU()    ### apply gain
+            #print('apply 1-16 zone by different gain ')
+            
+        else:
+            self.image_g=self.image_g/1.5
+            self.image_r=self.image_r/1.5
+            self.image_i=self.image_i/1.5
+            #print('apply gain finish')
+        ######################################################################
+        #size of the output image array, xsize is column, ysize is row, xsize = 9216,ysize = 9232
+
+        if self.information['xsize']==9216 and self.information['ysize']==9232 :
+            self.applyBias()  ### 
+            #print('apply bias finish')
+        else:
+            self.image_g=self.image_g+500.0
+            self.image_r=self.image_r+500.0
+            self.image_i=self.image_i+500.0
+            
+        ###################################################################    
+        if self.cosmetics and self.information['xsize']==9216 and self.information['ysize']==9232 :
+            self.applyCosmetics()
+            #print('applyCosmetics')
+        
+        ##################################################
+        if self.intscale:
+            self.discretise()                  
+        ##################################################
+        self.writeOutputs(simnumber)   
+        
+        self.log.info('Finished the simulation.') 
+        
+        #print('The iLoop= % d simlaiton finished. ' )
+
+       
+################################################################################################
+
+def mcisim(iLoop,configfile):
+    
+    opts, args = processArgs()                
+    opts.configfile=configfile
+    simulate = MCIsimulator(opts)  
+    simulate.configure(1)    
+
+    simulate.simulate(1)   
+       
+    return
+
+################################################################################################   
+################################################################################################
+################################################################################################
+################################### main function  #############################################
+
+#############################################################################
+
+if __name__ == "__main__":
+    
+    
+    
+    if  len(sys.argv[:]) <2:
+        
+        configfile='./mci_data/mci_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) 
+            
+            
+     
+
+    mcisim(1, configfile)
+    
+    print('MCI simulation is successful!')
+    
+   
+
+
+        
\ No newline at end of file

csst_mci_sim/mci_data/mci_sim_example.config

+[TEST]
+
+### user should define file path below
+inData_path =/home/yan/MCI_sim_Fabu/MCI_inputData
+
+result_path =/home/yan/MCI_sim_Fabu/mci_sim_result
+
+
+#size of the output image array, xsize is column, ysize is row, xsize = 9216,ysize = 9232
+xsize = 1000
+ysize = 1000 
+
+##prescan and overscan , do not change!!!
+prescan  = 27
+overscan = 320
+
+###sourcein = XDF, DARK, FLAT, BIAS
+
+sourcein=XDF
+
+##################################################
+####                                          ####    
+####Control flags (can be yes/no) ####
+##################################################
+
+cosmicRays     = yes
+
+darknoise      = yes 
+
+cosmetics      = yes
+
+radiationDamage= no
+
+bleeding       = yes
+
+overscans      = yes
+
+nonlinearity   = yes                           
+
+readoutNoise   = yes                             
+
+skyback        = yes
+
+TianceEffect   = yes
+
+intscale       = yes 
+
+ghosts         = yes
+
+shutterEffect      = yes                              
+
+flatfieldM         = yes 
+
+PRNUeffect         = yes 
+
+appFatt            = no 
+
+sky_shift_rot      = yes
+
+distortion         = no
+
+sim_star           = yes
+
+sim_galaxy         = yes   
+
+##############################################                             
+##############################################                            
+
+#CCD properties
+fullwellcapacity = 90000
+
+#dark noise in electrons per second
+dark = 0.001
+
+#exposure to simulate, exposure time
+exptime = 300.0
+
+###PNRU matrix sigma
+prnu_sigma=0.001
+
+### cosmicray coveringFraction
+coveringFraction=1.0
+
+
+#offset from the object, note that at the moment this is fixed, but in reality a focal plane position dependent.
+ghostOffsetX = 50
+ghostOffsetY = 50
+ghostRatio  = 1e-04
+
+
+
+####treering effect,bright fatter effect and difusion effect
+treering=0.1
+
+fatter=1.0
+
+difusion=0.1
+    
+### the choosen Filters in three CCDs   
+filter_g=u
+filter_r=F555W
+filter_i=F814W
+
+##['G_filters']=["F275W",  "F280N","NUV",   "WU",   "CBU",   "F343N",  "u",      "F373N", "F395N"," F336W"]
+##['R_filters']=["F487N", "F502N", "CBV",   "r",    "F656N", "F658N",  "F467M",  "F555W", "F606W",  "F673N"]        
+##['I_filters']=["z",      "y",    "F815N", "CBI",  "F925N", "F960M",  "F968N",  "F845M" ,"F850LP" ,"F814W"]
+
+####################
+g_gain1=1.53
+g_gain2=1.54
+g_gain3=1.55
+g_gain4=1.53
+g_gain5=1.51
+g_gain6=1.56
+g_gain7=1.58
+g_gain8=1.53
+g_gain9=1.54
+g_gain10=1.57
+g_gain11=1.51
+g_gain12=1.53
+g_gain13=1.55
+g_gain14=1.57
+g_gain15=1.53
+g_gain16=1.52
+
+#######################
+g_rdnois1=4.53
+g_rdnois2=4.54
+g_rdnois3=4.55
+g_rdnois4=4.53
+g_rdnois5=4.51
+g_rdnois6=4.56
+g_rdnois7=4.58
+g_rdnois8=4.53
+g_rdnois9=4.54
+g_rdnois10=4.57
+g_rdnois11=4.51
+g_rdnois12=4.53
+g_rdnois13=4.55
+g_rdnois14=4.57
+g_rdnois15=4.53
+g_rdnois16=4.52
+
+############################
+g_detbia1=500
+g_detbia2=510
+g_detbia3=514
+g_detbia4=520
+g_detbia5=524
+g_detbia6=540
+g_detbia7=530
+g_detbia8=532
+g_detbia9=534
+g_detbia10=526
+g_detbia11=532
+g_detbia12=516
+g_detbia13=540
+g_detbia14=560
+g_detbia15=536
+g_detbia16=528
+
+############################
+r_gain1=1.52
+r_gain2=1.54
+r_gain3=1.52
+r_gain4=1.55
+r_gain5=1.56
+r_gain6=1.52
+r_gain7=1.54
+r_gain8=1.6
+r_gain9=1.54
+r_gain10=1.52
+r_gain11=1.54
+r_gain12=1.57
+r_gain13=1.52
+r_gain14=1.55
+r_gain15=1.52
+r_gain16=1.55
+
+################################
+r_rdnois1=4.23
+r_rdnois2=4.24
+r_rdnois3=4.25
+r_rdnois4=4.23
+r_rdnois5=4.21
+r_rdnois6=4.26
+r_rdnois7=4.28
+r_rdnois8=4.23
+r_rdnois9=4.24
+r_rdnois10=4.27
+r_rdnois11=4.21
+r_rdnois12=4.23
+r_rdnois13=4.25
+r_rdnois14=4.27
+r_rdnois15=4.23
+r_rdnois16=4.22
+
+###################################
+r_detbia1=600
+r_detbia2=610
+r_detbia3=614
+r_detbia4=620
+r_detbia5=624
+r_detbia6=640
+r_detbia7=630
+r_detbia8=632
+r_detbia9=634
+r_detbia10=626
+r_detbia11=632
+r_detbia12=616
+r_detbia13=640
+r_detbia14=660
+r_detbia15=636
+r_detbia16=628
+
+###################################
+i_gain1=1.62
+i_gain2=1.64
+i_gain3=1.62
+i_gain4=1.65
+i_gain5=1.66
+i_gain6=1.62
+i_gain7=1.64
+i_gain8=1.6
+i_gain9=1.64
+i_gain10=1.62
+i_gain11=1.64
+i_gain12=1.67
+i_gain13=1.62
+i_gain14=1.65
+i_gain15=1.62
+i_gain16=1.65
+
+###################################
+i_rdnois1=4.63
+i_rdnois2=4.64
+i_rdnois3=4.65
+i_rdnois4=4.63
+i_rdnois5=4.61
+i_rdnois6=4.66
+i_rdnois7=4.68
+i_rdnois8=4.63
+i_rdnois9=4.64
+i_rdnois10=4.67
+i_rdnois11=4.61
+i_rdnois12=4.63
+i_rdnois13=4.65
+i_rdnois14=4.67
+i_rdnois15=4.63
+i_rdnois16=4.62
+
+#######################################
+i_detbia1=400
+i_detbia2=410
+i_detbia3=414
+i_detbia4=420
+i_detbia5=424
+i_detbia6=440
+i_detbia7=430
+i_detbia8=432
+i_detbia9=434
+i_detbia10=426
+i_detbia11=432
+i_detbia12=416
+i_detbia13=440
+i_detbia14=460
+i_detbia15=436
+i_detbia16=428
+
+####### end #######################
+
+
+
+
+
+        
+        
+
+
+
+
+

csst_mci_sim/support/MCIinstrumentModel.py

+"""
+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 matplotlib
+# import matplotlib.pyplot as plt
+# import datetime, math
+# import numpy as np
+# import numexpr as ne
+
+
+def MCIinformation():
+    """
+    Returns a dictionary describing MCI CCD. The following information is provided (id: value - reference)::
+
+        
+         dob: 0 - CDM03 (Short et al. 2010)
+         fwc: 90000 - CCD spec EUCL-EST-RS-6-002 (for CDM03)      
+         rdose: 30000000000.0 - derived (above the PLM requirement)
+         sfwc: 730000.0 - CDM03 (Short et al. 2010), see also the CCD spec EUCL-EST-RS-6-002       
+         st: 5e-06 - CDM03 (Short et al. 2010)        
+         svg: 1e-10 - CDM03 (Short et al. 2010)
+         t: 0.01024 - CDM03 (Short et al. 2010)
+         trapfile: cdm_euclid.dat - CDM03 (derived, refitted to CCD204 data)
+         vg: 6e-11 - CDM03 (Short et al. 2010)
+         vth: 11680000.0 - CDM03 (Short et al. 2010)
+
+
+
+    :return: instrument model parameters
+    :rtype: dict
+    """
+
+    #########################################################################################################
+    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()

csst_mci_sim/support/cosmicrays.py

+"""
+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__":
+    import sys
+    sys.path.append('/home/yan/csst-master/')
+    
+    
+    from support import logger as lg
+    from support import files as fileIO
+    from scipy import ndimage
+    from astropy.io import fits
+
+    #set up logger
+    log = lg.setUpLogger('VISsim.log')
+
+    #test section
+    crImage = np.zeros((2066, 2048), dtype=np.float64)
+
+    #cosmic ray instance
+    cosmics = cosmicrays(log, crImage)
+
+    #add cosmic rays up to the covering fraction
+    CCD_cr = cosmics.addUpToFraction(1.4, limit=None, verbose=True)
+    
+    print(CCD_cr)
+    print(type(CCD_cr))
+    
+    effected = np.count_nonzero(CCD_cr)
+    print (effected, effected*100./(CCD_cr.shape[0]*CCD_cr.shape[1]))
+
+    #save to FITS
+    fits.writeto(CCD_cr, 'cosmicrayTest.fits', overwrite=True)
+
+    #smooth with a charge diffusion kernel
+    smooth = ndimage.filters.gaussian_filter(CCD_cr, (0.32, 0.32))
+
+    fits.writeto(smooth, 'cosmicrayTestSmoothed.fits')

csst_mci_sim/support/logger.py

+"""
+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
+
+
+class SimpleLogger(object):
+    """
+    A simple class to create a log file or print the information on screen.
+    """
+
+    def __init__(self, filename, verbose=False):
+        self.file = open(filename, 'w')
+        self.verbose = verbose
+
+    def write(self, text):
+        """
+        Writes text either to file or screen.
+        """
+        print >> self.file, text
+        if self.verbose: print( text)

csst_mci_sim/support/sed.py

+#!/usr/bin/env python
+# -*- coding:utf-8 -*-
+
+# @Author: Shuai Feng (hebtu.edu.cn)
+# @Time:   2022-09-25
+
+import numpy as np
+from scipy.interpolate import interp1d
+from astropy.io import fits
+from astropy.io import ascii
+from astropy import units as u
+
+# ----------------
+# Magnitude Module
+
+from scipy.interpolate import interp1d
+
+def Calzetti_Law(wave, Rv = 4.05):
+    """Dust Extinction Curve by Calzetti et al. (2000)
+
+    Args:
+        wave (float): Wavelength
+        Rv (float, optional): Extinction curve. Defaults to 4.05.
+
+    Returns:
+        float: Extinction value E(B-V)
+    """
+    
+    wave_number = 1./(wave * 1e-4)
+    reddening_curve = np.zeros(len(wave))
+    
+    idx = np.logical_and(wave >= 1200, wave <= 6300)
+    reddening_curve[idx] = 2.659 * ( -2.156 + 1.509 * wave_number[idx] - 0.198 * \
+                    (wave_number[idx] ** 2)) + 0.011 * (wave_number[idx] **3 ) + Rv
+                                    
+    idx = np.logical_and(wave >= 6300, wave <= 22000)
+    reddening_curve[idx] = 2.659 * ( -1.857 + 1.040 * wave_number[idx]) + Rv
+    return reddening_curve
+
+def reddening(wave, flux, ebv = 0.0, law = 'calzetti', Rv = 4.05):
+    """
+    Reddening an input spectra through a given reddening curve.
+
+    Args:
+        wave (float): Wavelength of input spectra
+        flux (float): Flux of input spectra
+        ebv (float, optional): Extinction value. Defaults to 0.
+        law (str, optional): Extinction curve. Defaults to 'calzetti'.
+        Rv (float, optional): _description_. Defaults to 4.05.
+
+    Returns:
+        float: Flux of spectra after reddening.
+    """
+    if law == 'calzetti':
+        curve = Calzetti_Law(wave, Rv = Rv)
+        fluxNew = flux / (10. ** (0.4 * ebv * curve))
+    return fluxNew
+
+def flux_to_mag(wave, flux, filepath, band='GAIA_bp'):
+    """Convert flux of given spectra to magnitude
+
+    Args:
+        wave (float): Wavelength
+        flux (float): Flux of spectra
+        band (str, optional): Filter band name. Defaults to 'GAIA_bp'.
+
+    Returns:
+        float: value of magnitude
+    """
+    # /home/yan/MCI_sim/MCI_input/SED_Code/data
+    import os
+    
+    #parent = os.path.dirname(os.path.realpath(__file__))
+    
+    parent=filepath
+    
+    band = ascii.read(parent+'/SED_Code/seddata/' + band + '.dat')
+    
+    wave0= band['col1']
+    curv0= band['col2']
+    
+    # Setting the response
+    func = interp1d(wave0, curv0)
+    response = np.copy(wave)
+    ind_extra = (wave > max(wave0)) | (wave < min(wave0))
+    response[ind_extra] = 0
+    ind_inside = (wave < max(wave0)) & (wave > min(wave0))
+    response[ind_inside] = func(wave[ind_inside])
+
+    # Total Flux
+    Tflux = np.trapz(flux * response, wave) / np.trapz(response, wave)
+    
+    return -2.5 * np.log10(Tflux)
+
+def calibrate(wave, flux, mag, filepath,band='GAIA_bp'):
+    """
+    Calibrate the spectra according to the magnitude.
+
+    Args:
+        wave (float): Wavelength
+        flux (float): Flux of spectra
+        mag (float): Input magnitude.
+        band (str, optional): Filter band name. Defaults to 'GAIA_bp'.
+
+    Returns:
+        float: Flux of calibrated spectra. Units: 1e-17 erg/s/A/cm^2
+    """
+    inst_mag = flux_to_mag(wave, flux,filepath ,band = band)
+    instflux = 10 ** (-0.4 * inst_mag)
+    realflux = (mag * u.STmag).to(u.erg/u.s/u.cm**2/u.AA).value
+
+    # Normalization
+    flux_ratio = realflux / instflux
+    flux_calibrate = flux * flux_ratio * 1e17                        # Units: 10^-17 erg/s/A/cm^2
+    
+    return flux_calibrate
+
+# ------------
+# SED Template
+
+class Gal_Temp():
+    """
+    Template of Galaxy SED
+    """
+    
+    def __init__(self,parent):
+        import os
+        
+        self.parent = parent
+        
+        hdulist = fits.open(self.parent+'/SED_Code/seddata/galaxy_temp.fits')
+        self.wave = hdulist[1].data['wave']
+        self.flux = hdulist[2].data
+        self.age_grid = hdulist[3].data['logAge']
+        self.feh_grid = hdulist[3].data['FeH']
+        
+    def toMag(self, redshift = 0):
+        """Calculating magnitude
+
+        Args:
+            redshift (float, optional): redshift of spectra. Defaults to 0.
+        """
+        wave = self.wave * (1 + redshift)
+        self.umag = flux_to_mag(wave, self.flux, self.parent,band='SDSS_u')
+        self.gmag = flux_to_mag(wave, self.flux, self.parent,band='SDSS_g')
+        self.rmag = flux_to_mag(wave, self.flux, self.parent,band='SDSS_r')
+        self.imag = flux_to_mag(wave, self.flux, self.parent,band='SDSS_i')
+        self.zmag = flux_to_mag(wave, self.flux, self.parent,band='SDSS_z')
+        
+class Star_Temp():
+    """
+    Template of Stellar SED
+    """
+    
+    def __init__(self, parent):
+        
+        #import os
+        
+        #parent = os.path.dirname(os.path.realpath(__file__))
+
+        ##/home/yan/MCI_sim_Fabu/MCI_inputData/SED_Code
+        
+        self.parent = parent
+        
+        hdulist = fits.open(parent+'/SED_Code/seddata/stellar_temp.fits')
+        
+        self.wave = hdulist[1].data['wave']
+        self.flux = hdulist[2].data
+        self.Teff_grid = hdulist[3].data['Teff']
+        self.FeH_grid = hdulist[3].data['FeH']
+        
+    def toMag(self):
+        wave = self.wave
+        self.bpmag = flux_to_mag(wave, self.flux, self.parent,band='GAIA_bp')
+        self.rpmag = flux_to_mag(wave, self.flux, self.parent,band='GAIA_rp')
+        
+# -------------
+# SED Modelling
+
+def Model_Stellar_SED(wave, bp, rp, temp, filepath):
+    """Modelling stellar SED based on bp, rp magnitude
+
+    Args:
+        wave (float): Wavelength
+        bp (float): Magnitude of GAIA BP band
+        rp (float): Magnitude of GAIA RP band
+        temp (class): Class of stellar template
+
+    Returns:
+        float array: Spectral energy distribution of stellar SED, 
+        which have the same length to the input wave
+    """
+
+    color0 = bp - rp
+    colors = temp.bpmag - temp.rpmag
+    
+    idx = np.argmin(np.abs(colors - color0))
+    flux0 = temp.flux[idx]
+    flux1 = np.interp(wave, temp.wave, flux0)
+    flux  = calibrate(wave, flux1, rp, filepath, band = 'GAIA_rp')
+    return flux
+
+def Model_Galaxy_SED(wave, ugriz, z, temp):
+    """Modelling galaxy SED based on u,g,r,i,z magnitude
+
+    Args:
+        wave (float): Wavelength
+        ugriz (float, array): The array of magnitude of SDSS ugriz band
+        z (float): Redshift
+        temp (class): Class of gaalxy template
+
+    Returns:
+        float array: Spectral energy distribution of stellar SED, 
+        which have the same length to the input wave
+    """
+    sed = 10. ** (-0.4 * ugriz)
+    sed = sed / sed[2]
+    ntemp = len(temp.rmag)
+    dmag  = np.zeros(ntemp)
+    for j in range(ntemp):
+        ugriz0 = np.array([temp.umag[j], temp.gmag[j], temp.rmag[j], 
+                           temp.imag[j], temp.zmag[j]])
+        sed0 = 10. ** (-0.4 * ugriz0)
+        sed0 = sed0 / sed0[2]
+        dmag[j] = np.sum(np.abs(sed - sed0))
+    
+    idx = np.argmin(dmag)
+    flux0 = temp.flux[idx]
+    
+    # Effect of E(B-V)
+    ri0 = ugriz[2] - ugriz[3]
+    ri  = temp.rmag - temp.imag
+    dri = ri0 - ri[idx]
+    Alambda = Calzetti_Law(np.array([6213 / (1 + z), 7625 / (1 + z)]))
+    eri0 = (Alambda[0] - Alambda[1])
+    ebv = dri/eri0
+    if ebv<0:
+        ebv=0
+    if ebv>0.5:
+        ebv=0.5
+    flux1 = reddening(temp.wave, flux0, ebv = ebv)
+    
+    flux2 = np.interp(wave, temp.wave * (1 + z), flux1)
+    flux  = calibrate(wave, flux2, ugriz[2], band = 'SDSS_r')
+    return flux
\ No newline at end of file

csst_mci_sim/support/shao.py

+from ctypes import *
+
+def checkInputList(input_list, n):
+    if type(input_list) != type([1, 2, 3]):
+        raise TypeError("Input type is not list!", input_list)
+    for i in input_list:
+        if type(i) != type(1.1):
+            if type(i) != type(1):
+                raise TypeError("Input list's element is not float or int!", input_list)
+    if len(input_list) != n:
+        raise RuntimeError("Length of input list is not equal to stars' number!", input_list)
+            
+def onOrbitObsPosition(input_ra_list, input_dec_list, input_pmra_list, input_pmdec_list, input_rv_list, \
+                       input_parallax_list, input_nstars, input_x, input_y, input_z, input_vx, input_vy, \
+                       input_vz, input_epoch, input_date_str, input_time_str):
+    #Check input parameters
+    if type(input_nstars) != type(1):
+        raise TypeError("Parameter 7 is not int!", input_nstars)
+    
+    checkInputList(input_ra_list, input_nstars)
+    checkInputList(input_dec_list, input_nstars)
+    checkInputList(input_pmra_list, input_nstars)
+    checkInputList(input_pmdec_list, input_nstars)
+    checkInputList(input_rv_list, input_nstars)
+    checkInputList(input_parallax_list, input_nstars)
+    
+    if type(input_x) != type(1.1):
+        raise TypeError("Parameter 8 is not double!", input_x)
+    if type(input_y) != type(1.1):
+        raise TypeError("Parameter 9 is not double!", input_y)
+    if type(input_z) != type(1.1):
+        raise TypeError("Parameter 10 is not double!", input_z)
+    if type(input_vx) != type(1.1):
+        raise TypeError("Parameter 11 is not double!", input_vx)
+    if type(input_vy) != type(1.1):
+        raise TypeError("Parameter 12 is not double!", input_vy)
+    if type(input_vz) != type(1.1):
+        raise TypeError("Parameter 13 is not double!", input_vz)
+    
+    #Convert km -> m
+    input_x = input_x*1000.0 
+    input_y = input_y*1000.0 
+    input_z = input_z*1000.0
+    input_vx = input_vx*1000.0
+    input_vy = input_vy*1000.0
+    input_vz  = input_vz*1000.0
+    
+    if type(input_date_str) != type("2025-03-05"):
+        raise TypeError("Parameter 15 is not string!", input_date_str)
+    else:
+        input_date_str = input_date_str.strip()
+        if not (input_date_str[4]=="-" and input_date_str[7]=="-"):
+            raise TypeError("Parameter 15 format error (1)!", input_date_str)
+        else:
+            tmp = input_date_str.split("-")
+            if len(tmp) != 3:
+                raise TypeError("Parameter 15 format error (2)!", input_date_str)
+            input_year = int(tmp[0])
+            input_month = int(tmp[1])
+            input_day = int(tmp[2])
+            if not (input_year>=1900 and input_year<=2100):
+                raise TypeError("Parameter 15 year range error [1900 ~ 2100]!", input_year)
+            if not (input_month>=1 and input_month<=12):
+                raise TypeError("Parameter 15 month range error [1 ~ 12]!", input_month)
+            if not (input_day>=1 and input_day<=31):
+                raise TypeError("Parameter 15 day range error [1 ~ 31]!", input_day)
+    
+    if type(input_time_str) != type("20:15:15.15"):
+        raise TypeError("Parameter 16 is not string!", input_time_str)
+    else:
+        input_time_str = input_time_str.strip()
+        if not (input_time_str[2]==":" and input_time_str[5]==":"):
+            raise TypeError("Parameter 16 format error (1)!", input_time_str)
+        else:
+            tmp = input_time_str.split(":")
+            if len(tmp) != 3:
+                raise TypeError("Parameter 16 format error (2)!", input_time_str)
+            input_hour = int(tmp[0])
+            input_minute = int(tmp[1])
+            input_second = float(tmp[2])
+            if not (input_hour>=0 and input_hour<=23):
+                raise TypeError("Parameter 16 hour range error [0 ~ 23]!", input_hour)
+            if not (input_minute>=0 and input_minute<=59):
+                raise TypeError("Parameter 16 minute range error [0 ~ 59]!", input_minute)
+            if not (input_second>=0 and input_second<60.0):
+                raise TypeError("Parameter 16 second range error [0 ~ 60)!", input_second)
+    #Inital dynamic lib
+   
+    # import os
+    # currfile=os.getcwd()
+    # print(currfile)
+    
+    shao = cdll.LoadLibrary('./mci_so/libshao.so')
+    
+    shao.onOrbitObs.restype = c_int
+    
+    d3 = c_double * 3
+    shao.onOrbitObs.argtypes = [c_double, c_double, c_double, c_double, c_double, c_double, \
+                                c_int, c_int, c_int, c_int, c_int, c_double, \
+                                c_double, POINTER(d3), POINTER(d3), \
+                                c_int, c_int, c_int, c_int, c_int, c_double, \
+                                POINTER(c_double), POINTER(c_double) ]
+    output_ra_list = list()
+    output_dec_list = list()
+    for i in range(input_nstars):
+        input_ra = c_double(input_ra_list[i])
+        input_dec = c_double(input_dec_list[i])
+        input_pmra = c_double(input_pmra_list[i])
+        input_pmdec = c_double(input_pmdec_list[i])
+        input_rv = c_double(input_rv_list[i])
+        input_parallax = c_double(input_parallax_list[i])
+        p3 = d3(input_x, input_y, input_z)
+        v3 = d3(input_vx, input_vy, input_vz)
+        input_year_c = c_int(input_year)
+        input_month_c = c_int(input_month)
+        input_day_c = c_int(input_day) 
+        input_hour_c = c_int(input_hour)
+        input_minute_c = c_int(input_minute)
+        input_second_c  = c_double(input_second)
+        DAT = c_double(37.0)
+        output_ra = c_double(0.0)
+        output_dec = c_double(0.0)
+        rs = shao.onOrbitObs(input_ra, input_dec, input_parallax, input_pmra, input_pmdec, input_rv,  \
+                             input_year_c, input_month_c, input_day_c, input_hour_c, input_minute_c, input_second_c, \
+                             DAT, byref(p3), byref(v3), \
+                             input_year_c, input_month_c, input_day_c, input_hour_c, input_minute_c, input_second_c, \
+                             byref(output_ra), byref(output_dec))
+        if rs != 0:
+            raise RuntimeError("Calculate error!")
+        output_ra_list.append(output_ra.value)
+        output_dec_list.append(output_dec.value)
+    
+    return output_ra_list, output_dec_list
+    
\ No newline at end of file

setup.py

+import setuptools
+
+with open("README.md", "r") as fh:
+    long_description = fh.read()
+
+setuptools.setup(
+    name='csst_mci_sim',
+    version='2.0.0-MCI1.0.0',
+    author='CSST Team',
+    author_email='zhaojunyan@shao.ac.cn',
+    description='The CSST - mci - 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_mci_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_mci_sim': 'csst_mci_sim'},
+    # include_package_data=True,
+    package_data={"": ["LICENSE", "README.md"],
+                  "csst_mci_sim": ["mci_so/*", "mci_data/*"]},
+    # install_requires=['sphinx',
+    #                   'numpy',
+    #                   'scipy', 'matplotlib',
+    #                   'astropy', 'healpy', 'ccdproc', 'deepCR', 'photutils'],
+    python_requires='>=3.8',
+)