csst_mci_sim.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Oct 18 09:38:35 2021


@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 object catalog:



      #. 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 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.
         
       # add galaxy lensing effect to self.information['lensing'], set this parameter in .config file; update@24.10.16

.. 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
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.table import Table
from  astropy.io import fits
from astropy import units as u
import os, sys, math 
import configparser as ConfigParser
#from matplotlib import pyplot as plt 

from scipy import ndimage
######################
try:
    from support import logger as lg
    from support import cosmicrays
    from support import shao
    from support import sed 
    from support import MCIinstrumentModel
    from mci_so import cdm03bidir
    
except:    
    from .support import logger as lg
    from .support import cosmicrays
    from .support import shao
    from .support import sed 
    from .support import MCIinstrumentModel
    from .mci_so import cdm03bidir

####################################################
#from joblib import Parallel, delayed
from astropy.coordinates import SkyCoord
from scipy import interpolate
from scipy.signal import fftconvolve
import pandas
import ctypes
import astropy.coordinates as coord
from scipy.interpolate import interp1d

########################### functions  #########################

"""
Charge Transfer Inefficiency
============================

This file contains a simple class to run a CDM03 CTI model developed by Alex Short (ESA).

This now contains both the official CDM03 and a new version that allows different trap
parameters in parallel and serial direction.

:requires: NumPy
:requires: CDM03 (FORTRAN code, f2py -c -m cdm03bidir cdm03bidir.f90)

:version: 0.35
"""
import numpy as np


#CDM03bidir
class CDM03bidir():
    """
    Class to run CDM03 CTI model, class Fortran routine to perform the actual CDM03 calculations.

     :param settings: input parameters
     :type settings: dict
     :param data: input data to be radiated
     :type data: ndarray
     :param log: instance to Python logging
     :type log: logging instance
    """
    def __init__(self, settings, data, log=None):
        """
        Class constructor.

        :param settings: input parameters
        :type settings: dict
        :param data: input data to be radiated
        :type data: ndarray
        :param log: instance to Python logging
        :type log: logging instance
        """
        self.data = data
        self.values = dict(quads=(0,1,2,3), xsize=2048, ysize=2066, dob=0.0, rdose=8.0e9)
        self.values.update(settings)
        self.log = log
        self._setupLogger()

        #default CDM03 settings
        self.params = dict(beta_p=0.6, beta_s=0.6, fwc=200000., vth=1.168e7, vg=6.e-11, t=20.48e-3,
                           sfwc=730000., svg=1.0e-10, st=5.0e-6, parallel=1., serial=0.0)
        #update with inputs
        self.params.update(self.values)

        #read in trap information
        trapdata = np.loadtxt(self.values['dir_path']+self.values['paralleltrapfile'])
        if trapdata.ndim > 1:
            self.nt_p = trapdata[:, 0]
            self.sigma_p = trapdata[:, 1]
            self.taur_p = trapdata[:, 2]
        else:
            #only one trap species
            self.nt_p = [trapdata[0],]
            self.sigma_p = [trapdata[1],]
            self.taur_p = [trapdata[2],]

        trapdata = np.loadtxt(self.values['dir_path']+self.values['serialtrapfile'])
        if trapdata.ndim > 1:
            self.nt_s = trapdata[:, 0]
            self.sigma_s = trapdata[:, 1]
            self.taur_s = trapdata[:, 2]
        else:
            #only one trap species
            self.nt_s = [trapdata[0],]
            self.sigma_s = [trapdata[1],]
            self.taur_s = [trapdata[2],]

        #scale thibaut's values
        if 'thibaut' in self.values['parallelTrapfile']:
            self.nt_p /= 0.576  #thibaut's values traps / pixel
            self.sigma_p *= 1.e4 #thibaut's values in m**2
        if 'thibaut' in self.values['serialTrapfile']:
            self.nt_s *= 0.576 #thibaut's values traps / pixel  #should be division?
            self.sigma_s *= 1.e4 #thibaut's values in m**2


    def _setupLogger(self):
        """
        Set up the logger.
        """
        self.logger = True
        # if self.log is None:
        #     self.logger = False


    def applyRadiationDamage(self, data, iquadrant=0):
        """
        Apply radian damage based on FORTRAN CDM03 model. The method assumes that
        input data covers only a single quadrant defined by the iquadrant integer.

        :param data: imaging data to which the CDM03 model will be applied to.
        :type data: ndarray

        :param iquandrant: number of the quadrant to process
        :type iquandrant: int

        cdm03 - Function signature::

              sout = cdm03(sinp,iflip,jflip,dob,rdose,in_nt,in_sigma,in_tr,[xdim,ydim,zdim])
            Required arguments:
              sinp : input rank-2 array('d') with bounds (xdim,ydim)
              iflip : input int
              jflip : input int
              dob : input float
              rdose : input float
              in_nt : input rank-1 array('d') with bounds (zdim)
              in_sigma : input rank-1 array('d') with bounds (zdim)
              in_tr : input rank-1 array('d') with bounds (zdim)
            Optional arguments:
              xdim := shape(sinp,0) input int
              ydim := shape(sinp,1) input int
              zdim := len(in_nt) input int
            Return objects:
              sout : rank-2 array('d') with bounds (xdim,ydim)

        .. Note:: Because Python/NumPy arrays are different row/column based, one needs
                  to be extra careful here. NumPy.asfortranarray will be called to get
                  an array laid out in Fortran order in memory. Before returning the
                  array will be laid out in memory in C-style (row-major order).

        :return: image that has been run through the CDM03 model
        :rtype: ndarray   
        """""
        #return data
    
        iflip = iquadrant / 2
        jflip = iquadrant % 2

        params = [self.params['beta_p'], self.params['beta_s'], self.params['fwc'], self.params['vth'],
                  self.params['vg'], self.params['t'], self.params['sfwc'], self.params['svg'], self.params['st'],
                  self.params['parallel'], self.params['serial']]

        if self.logger:
            self.log.info('nt_p=' + str(self.nt_p))
            self.log.info('nt_s=' + str(self.nt_s))
            self.log.info('sigma_p= ' + str(self.sigma_p))
            self.log.info('sigma_s= ' + str(self.sigma_s))
            self.log.info('taur_p= ' + str(self.taur_p))
            self.log.info('taur_s= ' + str(self.taur_s))
            self.log.info('dob=%f' % self.values['dob'])
            self.log.info('rdose=%e' % self.values['rdose'])
            self.log.info('xsize=%i' % data.shape[1])
            self.log.info('ysize=%i' % data.shape[0])
            self.log.info('quadrant=%i' % iquadrant)
            self.log.info('iflip=%i' % iflip)
            self.log.info('jflip=%i' % jflip)

#################################################################################

        
        CTIed = cdm03bidir.cdm03(np.asfortranarray(data),
                                  jflip, iflip,
                                  self.values['dob'], self.values['rdose'],
                                  self.nt_p, self.sigma_p, self.taur_p,
                                  self.nt_s, self.sigma_s, self.taur_s,
                                  params,
                                  [data.shape[0], data.shape[1], len(self.nt_p), len(self.nt_s), len(self.params)])
       
        
        return np.asanyarray(CTIed)
       
#################################################################################################################
#################################################################################################################


def transRaDec2D(ra, dec):
    # radec转为竞天程序里的ob, 赤道坐标系下的笛卡尔三维坐标xyz.
    x1 = np.cos(dec / 57.2957795) * np.cos(ra / 57.2957795)
    y1 = np.cos(dec / 57.2957795) * np.sin(ra / 57.2957795)
    z1 = np.sin(dec / 57.2957795)
    return np.array([x1, y1, z1])
###############################################################################

################################################################
def ill2flux(E,path):

    # use template from sky_bkg (background_spec_hst.dat)
    filename = path+'MCI_inputData/refs/background_spec_hst.dat'
    cat_spec = pd.read_csv(filename, sep='\s+', header=None, comment='#')
    wave0 = cat_spec[0].values       # A
    spec0 = cat_spec[2].values      # erg/s/cm^2/A/arcsec^2

    # convert erg/s/cm^2/A/arcsec^2 to erg/s/cm^2/A/sr
    flux1 = spec0 / (1/4.25452e10)
    # convert erg/s/cm^2/A/sr to W/m^2/sr/um
    flux2 = flux1 * 10

    # 对接收面积积分，输出单位 W/m^2/nm
    D = 2   # meter
    f = 28  # meter, 焦距，转换关系来源于王维notes.
    flux3 = flux2 * np.pi * D**2 / 4 / f**2 / 10**3

    f = interp1d(wave0, flux3)
    wave_range = np.arange(3800, 7800)
    flux3_mean = f(wave_range)
    delta_lamba = 0.1   # nm
    E0 = np.sum(flux3_mean * delta_lamba)

    factor = E / E0
    spec_scaled = factor * spec0

    return wave0, spec_scaled

##############################################################

def earth_angle(time_jd, x_sat, y_sat, z_sat, ra_obj, dec_obj):

    ra_sat = np.arctan2(y_sat, x_sat) / np.pi * 180
    dec_sat = np.arctan2(z_sat, np.sqrt(x_sat**2+y_sat**2)) / np.pi * 180
    radec_sat = SkyCoord(ra=ra_sat*u.degree, dec=dec_sat*u.degree, frame='gcrs')
    lb_sat = radec_sat.transform_to('geocentrictrueecliptic')

    # get the obj location
    radec_obj = SkyCoord(ra=ra_obj*u.degree, dec=dec_obj*u.degree, frame='gcrs')
    lb_obj = radec_obj.transform_to('geocentrictrueecliptic')

    # calculate the angle between sub-satellite point and the earth side
    earth_radius = 6371     # km
    sat_height = np.sqrt(x_sat**2 + y_sat**2 + z_sat**2)
    angle_a = np.arcsin(earth_radius/sat_height) / np.pi * 180

    # calculate the angle between satellite position and the target position
    angle_b = lb_sat.separation(lb_obj)

    # calculat the earth angle
    angle = 180 - angle_a - angle_b.degree

    return angle
#######################################

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
#############################################################33333

class StrayLight(object):

    def __init__(self, path, jtime = 2460843., sat = np.array([0,0,0]), radec = np.array([0,0])):
        self.path = path
        self.jtime = jtime
        self.sat = sat
        self.equator = coord.SkyCoord(radec[0]*u.degree, radec[1]*u.degree,frame='icrs')
        self.ecliptic = self.equator.transform_to('barycentrictrueecliptic')
        self.pointing = transRaDec2D(radec[0], radec[1])
        self.slcdll = ctypes.CDLL(self.path+'MCI_inputData/refs/libstraylight.so') #dylib

        self.slcdll.Calculate.argtypes = [ctypes.c_double, ctypes.POINTER(ctypes.c_double),
                                         ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_double),
                                         ctypes.POINTER(ctypes.c_double), ctypes.c_char_p]

        self.slcdll.PointSource.argtypes = [ctypes.c_double ,ctypes.POINTER(ctypes.c_double),
                                           ctypes.POINTER(ctypes.c_double),ctypes.POINTER(ctypes.c_double),
                                           ctypes.POINTER(ctypes.c_double), ctypes.c_char_p]

        self.slcdll.EarthShine.argtypes = [ctypes.c_double ,ctypes.POINTER(ctypes.c_double),
                                          ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_double),
                                          ctypes.POINTER(ctypes.c_double)]

        self.slcdll.Zodiacal.argtypes = [ctypes.c_double ,ctypes.POINTER(ctypes.c_double),
                                        ctypes.POINTER(ctypes.c_double)]
        self.slcdll.ComposeY.argtypes = [ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_double),
                                       ctypes.POINTER(ctypes.c_double)]
        self.slcdll.Init.argtypes=[ctypes.c_char_p,ctypes.c_char_p,ctypes.c_char_p,ctypes.c_char_p]
        self.deFn = self.path+"MCI_inputData/refs/DE405"
        self.PSTFn = self.path+"MCI_inputData/refs/PST"
        self.RFn = self.path+"MCI_inputData/refs/R"
        self.ZolFn =self.path+"MCI_inputData/refs/Zodiacal"
        self.brightStarTabFn = self.path+"MCI_inputData/refs/BrightGaia_with_csst_mag"

        self.slcdll.Init(str.encode(self.deFn),str.encode(self.PSTFn),str.encode(self.RFn),str.encode(self.ZolFn))


    def caculateStarLightFilter(self, filter='i'):
        
        filterIndex = {'nuv':0,'u':1,'g':2,'r':3, 'i':4, 'z':5,'y':6}

        sat = (ctypes.c_double*3)()
        sat[:] = self.sat
        ob = (ctypes.c_double*3)()
        ob[:] = self.pointing

        py1 = (ctypes.c_double*3)()
        py2 = (ctypes.c_double*3)()
        self.slcdll.ComposeY(ob, py1, py2)

        star_e1 = (ctypes.c_double*7)()
        self.slcdll.PointSource(self.jtime, sat, ob, py1, star_e1, str.encode(self.brightStarTabFn))
        star_e2 = (ctypes.c_double*7)()
        self.slcdll.PointSource(self.jtime, sat, ob, py2, star_e2, str.encode(self.brightStarTabFn))

        band_star_e1 = star_e1[:][filterIndex[filter]]
        band_star_e2 = star_e2[:][filterIndex[filter]]

        return max(band_star_e1, band_star_e2)

    def caculateEarthShineFilter(self, filter='i'):
        
        filterIndex = {'nuv':0,'u':1,'g':2,'r':3, 'i':4, 'z':5,'y':6}
        sat = (ctypes.c_double*3)()
        sat[:] = self.sat
        ob = (ctypes.c_double*3)()
        ob[:] = self.pointing

        py1 = (ctypes.c_double*3)()
        py2 = (ctypes.c_double*3)()
        self.slcdll.ComposeY(ob, py1, py2)

        earth_e1 = (ctypes.c_double*7)()
        self.slcdll.EarthShine(self.jtime, sat, ob, py1, earth_e1)          # e[7]代表7个波段的照度
        
        earth_e2 = (ctypes.c_double*7)()
        self.slcdll.EarthShine(self.jtime, sat, ob, py2, earth_e2)

        band_earth_e1 = earth_e1[:][filterIndex[filter]]
        band_earth_e2 = earth_e2[:][filterIndex[filter]]

        return max(band_earth_e1, band_earth_e2)
    
###############################################################################  
###############################################################################
###############################################################################

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

##############################################################################    

##############################################################################
    
#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(ra0, dec0):
    '''convert deg to ra's HMS and  dec's DMS'''
    c = SkyCoord(ra=ra0*u.degree, dec=dec0*u.degree)
    ss=c.to_string('hmsdms')
    
    ## print(ss)
    
    for k in range(5,len(ss)-1):

        if ss[k]=='s':
            s_idx=k
            
        if ss[k]=='d':
            d_idx=k
        
        if ss[k]=='m':
            m_idx=k
            
    if c.ra.hms.s<10 and  c.ra.hms.s>0.1:
        temp=str(c.ra.hms.s)
        ra_ss='0'+temp[:3]
    if c.ra.hms.s>10:
        temp=str(c.ra.hms.s)
        ra_ss=temp[:4]
        
    if c.ra.hms.s<0.1:
        temp=str(c.ra.hms.s)
        ra_ss='00.0'
        
    dms_d=ss[s_idx+2:d_idx]
    dms_m=ss[d_idx+1:m_idx]
    dms_s=ss[m_idx+1:m_idx+3]
    
    hhmmss=ss[0:2]+ss[3:5]+ra_ss
    return hhmmss+dms_d+dms_m+dms_s

###############################################################################

###############################################################################
#############################################################################
#### 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, dn=5, IDWindex=3, 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
    
    ####################################
    #######   my code  ######
    psfWeight = np.zeros([dn])

    if OnlyNeighbors == True:
        
        neigh = findNeighbors(px, py, cen_col, cen_row, dn)
        
    #######################################################################
    for ipsf in range(len(neigh)):

        idx=neigh[ipsf]
        dist = np.sqrt((ref_col - cen_col[idx])**2 + (ref_row - cen_row[idx])**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
        if IDWindex == 5:
            psfWeight[ipsf] = dist**5
            
        psfWeight[ipsf] = max(psfWeight[ipsf], minimum_psf_weight)
        psfWeight[ipsf] = 1./psfWeight[ipsf]
        
    ##################################3    
        
    psfWeight /= np.sum(psfWeight)

    psfMaker  = np.zeros([ngy, ngx, 7 ], dtype=np.float32)
    
    #########################################################
    for waven in range(7):
    
        for ipsf in range(len(neigh)):
    
            idx=neigh[ipsf]
    
            iPSFMat = PSFMat[:, :, waven,idx].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, configfile):
        """
        Class Constructor.

        :param opts: OptionParser instance
        :type opts: OptionParser instance
        """
        self.configfile = configfile
        self.section = 'TEST'       #####simulation 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,
                                     pixel_size=0.05,
                                     dark=0.001,
                                     exptime=300.0,
                                     readouttime=4.,
                                     rdose=8.0e9,                                  
                                     ghostCutoff=22.0,
                                     ghostratio=5.e-5,
                                     coveringfraction=0.05,    #CR: 
                                     # cosmicraylengths =self.information['dir_path']+'MCI_inputData/data/cdf_cr_length.dat',
                                     # cosmicraydistance=self.information['dir_path']+'MCI_inputData/data/cdf_cr_total.dat',                                
                                     # parallelTrapfile =self.information['dir_path']+'MCI_inputData/data/cdm_euclid_parallel.dat',
                                     # serialTrapfile   =self.information['dir_path']+'MCI_inputData/data/cdm_euclid_serial.dat',
                                     mode='same'))
    
####################################################

###############################################################################
    def readConfigs(self,simnumber,source):
        """
        Reads the config file information using configParser and sets up a logger.
        """
        self.config = ConfigParser.RawConfigParser()
        
        if simnumber==1:
            print('beging config file name : ')
            print(self.configfile)        
        
        #self.config.readfp(open(self.configfile))
        self.config.read_file(open(self.configfile))
        

       
###########################################################################

    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...
        self.information['ghostratio'] = float(self.config.get(self.section, 'ghostRatio'))

    ###################################################################################################
    ##################################################################################################
    
        #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')
        
        self.save_starpsf         =self.config.getboolean(self.section, 'save_starpsf')
        self.save_cosmicrays      =self.config.getboolean(self.section, 'save_cosmicrays')
        self.lensing              =self.config.getboolean(self.section, 'lensing')
        
        ###############################################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,
                             save_starpsf       =self.save_starpsf,
                             save_cosmicrays    =self.save_cosmicrays ,
                             lensing            =self.lensing )        
        #####################################################################  


            
        return
            
 ##############################################################################           
###############################################################################

    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=float)
        self.image_r=np.zeros((self.information['ysize'], self.information['xsize']), dtype=float)
        self.image_i=np.zeros((self.information['ysize'], self.information['xsize']), dtype=float)
        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']