camera.py

import math
import numpy as np
import scipy.ndimage as nd
from astropy.io import fits
import matplotlib.pyplot as plt

from .config import config, S
from .utils import region_replace, random_seed_select
from .io import log
from .optics import filter_throughput

cpism_refdata = config['cpism_refdata']
MAG_SYSTEM = config['mag_system']
solar_spectrum = S.FileSpectrum(config['solar_spectrum'])
solar_spectrum.convert('photlam')
def sky_frame_maker(band, skybg, platescale, shape):
    """
    generate a sky background frame.

    Parameters
    ----------
    band : str
        The band of the sky background.
    skybg : str
        The sky background file name. 
    platescale : float
        The platescale of the camera in arcsec/pixel.
    shape : tuple
        The shape of the output frame. (y, x)

    Returns
    -------
    sky_bkg_frame : numpy.ndarray
        The sky background frame.
    """
    filter = filter_throughput(band)
    sk_spec = solar_spectrum.renorm(skybg, MAG_SYSTEM, filter)
    sky_bkg_frame = np.zeros(shape)
    sky_bkg_frame += (sk_spec * filter).integrate() * platescale**2
    return sky_bkg_frame


class CRobj(object):
    """
    Cosmic ray object.

    Attributes
    ----------
    flux : float
        The flux of the cosmic ray.
    angle : float
        The angle of the cosmic ray.
    sigma : float
        The width of the cosmic ray.
    length : int
        The length of the cosmic ray.
    """

    def __init__(self, flux, angle, sigma, length) -> None:
        self.flux = flux
        self.angle = angle
        self.sigma = sigma
        self.length = length


class CosmicRayFrameMaker(object):
    """
    Cosmic ray frame maker.

    Parameters
    ----------
    depth : float
        The depth of the camera pixel in um.
    pitch : float
        The pitch of the camera pixel in um.
    cr_rate : float
        The cosmic ray rate per second per cm2.

    """

    def __init__(self) -> None:
        self.tmp_size = [7, 101]
        self.freq_power = -0.9
        self.trail_std = 0.1
        self.depth = 10  # um
        self.pitch = 13  # um
        self.cr_rate = 1  # particle per s per cm2 from Miles et al. 2021

    def make_CR(self, length, sigma, seed=-1):
        """
        make a image of cosmic ray with given length and sigma.

        Parameters
        ----------
        length : int
            The length of the cosmic ray in pixel.
        sigma : float
            The width of the cosmic ray in pixel.

        Returns
        -------
        output : numpy.ndarray
            The image of cosmic ray.
        """
        h = self.tmp_size[0]
        w = self.tmp_size[1]

        freq = ((w-1)/2-np.abs(np.arange(w)-(w-1)/2)+1)**(self.freq_power)

        x = np.arange(w) - (w-1)/2
        hl = (length-1)/2
        x_wing = np.exp(-(np.abs(x)-hl)**2/sigma/sigma/2)
        x_wing[np.abs(x) < hl] = 1

        cr = np.zeros([h, w])
        center = (h-1)/2

        for i in range(h):
            phase = np.random.rand(w)*2*np.pi
            trail0 = abs(np.fft.fft(freq*np.sin(phase) + 1j*x*np.cos(phase)))
            # TODO maybe somthing wrong
            trail_norm = (trail0 - trail0.mean())/trail0.std()
            cr[i, :] = np.exp(-(i - center)**2/sigma/sigma/2) \
                * (trail_norm * self.trail_std + 1) * x_wing

        output = np.zeros([w, w])
        d = (w-h)//2
        output[d:d+h, :] = cr
        return output

    def _length_rand(self, N, seed=-1):
        """
        randomly generate N cosmic ray length.
        """
        len_out = []
        seed = random_seed_select(seed=seed)
        log.debug(f"cr length seed: {seed}")
        for i in range(N):
            x2y2 = 2
            while x2y2 > 1:
                lx = 1 - 2 * np.random.rand()
                ly = 1 - 2 * np.random.rand()
                x2y2 = lx * lx + ly * ly

            z = 1 - 2 * x2y2
            r = 2 * np.sqrt(x2y2 * (1 - x2y2))
            length = abs(r / z * self.depth)
            pitch = self.pitch

            len_out.append(int(length/pitch))
        return np.array(len_out)

    def _number_rand(self, expt, pixsize, random=False, seed=-1):
        """
        randomly generate the number of cosmic rays.
        """
        area = (self.pitch / 1e4)**2 * pixsize[0] * pixsize[1]
        ncr = area * expt * self.cr_rate
        if random:
            seed = random_seed_select(seed=seed)
            log.debug(f"cr count: {seed}")
            ncr = np.random.poisson(ncr)
        else:
            ncr = int(ncr)
        self.ncr = ncr
        return ncr

    def _sigma_rand(self, N, seed=-1):
        """
        randomly generate N cosmic ray sigma.
        """
        sig_sig = 0.5  # asuming the sigma of size of cosmic ray is 0.5
        seed = random_seed_select(seed=seed)
        log.debug(f"cr width seed: {seed}")
        sig = abs(np.random.randn(N))*sig_sig + 1/np.sqrt(12) * 1.2
        # assume sigma is 1.2 times of pictch sig
        return sig

    def _flux_rand(self, N, seed=-1):
        """
        randomly generate N cosmic ray flux.
        """
        seed = random_seed_select(seed=seed)
        log.debug(f"cr flux seed: {seed}")
        u = np.random.rand(N)
        S0 = 800
        lam = 0.57
        S = (-np.log(1-u)/lam + S0**0.25)**4
        return S

    def random_CR_parameter(self, expt, pixsize):
        """
        randomly generate cosmic ray parameters, including number, length, flux, sigma and angle.

        Parameters
        ----------
        expt : float
            The exposure time in second.
        pixsize : list
            The size of the image in pixel.

        Returns
        -------
        CRs : list
            A list of cosmic ray objects.

        """
        N = self._number_rand(expt, pixsize)
        log.debug(f"cr count: {N}")
        length = self._length_rand(N)
        if N > 0:
            log.debug(f"cr length, max: {length.max()}, min: {length.min()}")
            flux = self._flux_rand(N)
            log.debug(f"cr flux, max: {flux.max()}, min: {flux.min()}")
            sig = self._sigma_rand(N)
            log.debug(f"cr width, max: {sig.max()}, min: {sig.min()}")
            seed = random_seed_select(seed=-1)
            log.debug(f"cr angle seed: {seed}")
            angle = np.random.rand(N) * 180

        CRs = []
        for i in range(N):
            CRs.append(CRobj(flux[i], angle[i], sig[i], length[i]))
        return CRs

    def make_cr_frame(self, shape, expt, seed=-1):
        """
        make a cosmic ray frame.

        Parameters
        ----------
        shape : list
            The size of the image in pixel.
        expt : float
            The exposure time in second.
        seed : int, optional
            The random seed. The default is -1. If seed is -1, the seed will be randomly selected.

        Returns
        -------
        image : numpy.ndarray
            The cosmic ray frame.
        """
        image = np.zeros(shape)
        sz = shape
        cr_array = self.random_CR_parameter(expt, shape)
        cr_center = (self.tmp_size[1] - 1)/2
        seed = random_seed_select(seed=seed)
        log.debug(f"cr position seed: {seed}")

        for i in range(len(cr_array)):
            cr = cr_array[i]
            x = np.random.rand() * sz[1]
            y = np.random.rand() * sz[0]
            cr_img = self.make_CR(cr.length, cr.sigma)
            cr_img *= cr.flux
            cr_img = abs(nd.rotate(cr_img, cr.angle, reshape=False))

            if i == 0:
                pdin = False
            else:
                pdin = True

            if i == len(cr_array) - 1:
                pdout = False
            else:
                pdout = True

            image = region_replace(
                image, cr_img,
                [y-cr_center, x-cr_center],
                padded_in=pdin,
                padded_out=pdout, subpix=True
            )
            image = np.maximum(image, 0)
        log.debug(f"cr image max: {image.max()}, min: {image.min()}")
        return image


class CpicVisEmccd(object):
    """
    The class for Cpic EMCCD camera.

    Attributes
    -----------
    switch : dict
        A dictionary to switch on/off the camera effects.
    """
    def __init__(self, config_dict=None):
        """Initialize the camera.
        
        Parameters
        ----------
        config_dict : dict, optional
            The configuration dictionary.

        Returns
        -------
        None
        """
        self._defaut_config()
        if config_dict is not None:
            old_switch = self.switch
            self.__dict__.update(config_dict) #not safe, be careful
            old_switch.update(self.switch)
            self.switch = old_switch
        self.config_init()

    def _defaut_config(self):
        """set up defaut config for the camera."""
        self.plszx = 1024
        self.plszy = 1024
        self.pscan1 = 8
        self.pscan2 = 0
        self.oscan1 = 16
        self.oscan2 = 18
        self.udark = 6
        self.bdark = 2
        self.ldark = 16
        self.rdark = 16

        self.max_adu = 16_383

        self.switch = {
            'bias_vp': True,
            'bias_hp': True,
            'bias_ci': True,
            'bias_shift': True,
            'cic': True,
            'dark': True,
            'flat': True,
            'badcolumn': True,
            'shutter': True,
            'cte': True,
            'nonlinear': True,
            'cosmicray': True,
            'blooming': True,
            'em_blooming': True,
            }

        self.dark_file = cpism_refdata + '/camera/emccd_dark_current.fits'
        self.flat_file = cpism_refdata + '/camera/emccd_flat_field.fits'
        self.cic_file = cpism_refdata + '/camera/emccd_cic2.fits'
        self.bad_col_file = cpism_refdata + '/camera/emccd_bad_columns.fits'
        self.cic = None
        self.dark = None
        self.flat = None
        
        self.fullwell = 80_000
        self.em_fullwell = 500_000 #780_000
        self.em_cte = 0.9996
        self.emreg_cal_num = 10 # 用来加速计算
        self.emreg_num = 604
        self.readout_speed = 6.25e6 # Hz
        self.readout_time = 0.365 # s

        self.heat_speed = 1 / 1000 # voltage / 1000 degree per frame
        self.temper_speed = 0.05 # degree per second
        self.cooler_temp = -80

        self.readout_noise = 160
        self.ph_per_adu = 59
        self.bias_level = 200
        self.vertical_param = [0, 14]
        self.horizontal1_param = [5.3/2, 12, 5/2, 50]
        self.horizontal2_param = [5.3/4, 16.17, 5/4, 76.65]
        self.bias_hp_resd = 2
        self.cooler_interfence = 20
        self.bias_level_std = 3
        self.bias_shift_per_volt = -1.34
        self.shift_time = 1 / 1000

        self.nonlinear_coefficient = -0.1

        self.detector_name = 'EMCCD'
        self.ccd_label= 'CCD201-20'
        self.pitch_size = 13

        
    
    def config_init(self):
        """initialize the camera. 
            If the config is set, call this function to update the config.
        """

        self._img_index = 0
        self._ccd_temp = self.cooler_temp
        self._vertical_i0 = np.random.randint(0, 2)
        self._frame_read_time = 2080 * 1055 / 6.25e6

        self.flat_shape = [self.plszy, self.plszx]

        darksz_x = self.plszx + self.rdark + self.ldark
        darksz_y = self.plszy + self.udark + self.bdark
        self.dark_shape = [darksz_y, darksz_x]

        biassz_x = darksz_x + self.pscan1 + self.oscan1
        biassz_y = darksz_y + self.pscan2 + self.oscan2
        self.bias_shape = [biassz_y, biassz_x]

        if self.flat is None:
            self.flat = fits.getdata(self.flat_file)
        if self.cic is None:
            self.cic = fits.getdata(self.cic_file)
        if self.dark is None:
            self.dark = fits.getdata(self.dark_file)

        self.bad_col = fits.getdata(self.bad_col_file)
        self.ccd_temp = self.cooler_temp
        self.system_time = 0
        self.time_syn(0, initial=True)
        self.emgain_set(1024, -80)

        def CTE_cal(n, N_p, CTE, S_0):
            '''
            CTE_cal(order of trail pixels, number of pixel transfers, CTE, initial intensity of target pixel)
            '''
            CTI = 1 - CTE
            S_0_n = S_0 * ((N_p * CTI) ** n) / math.factorial(n) * math.exp(-N_p * CTI)
            return S_0_n
        
        def cte_201(cte, start=0, length=10):
            N_p = 604
            S_0 = 1
            res = [0]*start
            for n in range(length):
                s = CTE_cal(n, N_p, cte, S_0)
                res.append(s)
            return np.array(res)

        cti_trail = cte_201(self.em_cte, start=0, length=10)
        self.cti_trail = cti_trail / cti_trail.sum()


    def em_fix_fuc_fit(self, emgain):
        """Calculate the emgain fix coeficient to fix the gamma distribution.
        The coeficient is from fixing of ideal emgain distribution.

        Parameters
        ----------
        emgain : float
            The emgain.
        
        Returns
        -------
        float
            The coeficient.
        """
        emgain = np.array([emgain]).flatten()
        p = [0.01014486, -0.00712984, -0.17163414,  0.09523666, -0.53926089]
        def kernel(em):
            log_em = np.log10(em)
            loglog_g = np.log10(log_em)
            loglog_t = np.polyval(p, loglog_g)
            log_t = 10**loglog_t
            t = 10**log_t - 1
            return t
        output = []
        for em in emgain:
            if em <= 1:
                output.append(0)
            elif em > 80:
                output.append(kernel(80))
            else:
                output.append(kernel(em))
        return np.array(output)


    def bias_frame(self):
        """Generate bias frame
        The bias frame contains vertical, horizontal, peper-salt noise, bias drift effect. 
        Can be configurable using self.switch.

        Returns
        -------
        np.ndarray
            bias frame
        """

        shape = self.bias_shape
        TPI = np.pi * 2
        
        # vertical pattern
        # 使用一维的曲线描述竖条纹的截面
        vp_1d = np.zeros(shape[1])
        # 以下代码用于模拟垂直间隔的竖条纹在奇偶幅图像时的不同表现
        # 后续相机更新过程中，已经将改特性进行修改，因此不再使用此代码
        # vp_1d[0::2] = self.vertical_param[self._vertical_i0]
        # self._vertical_i0 = 1 - self._vertical_i0
        # vp_1d[1::2] = self.vertical_param[self._vertical_i0]
        vp_1d[0::2] = self.vertical_param[self._vertical_i0]
        vp_1d[1::2] = self.vertical_param[1 - self._vertical_i0]

        vp_frame = np.zeros(shape)
        if self.switch['bias_vp']:
            vp_frame += vp_1d

        # if show: # pragma: no cover
        #     plt.figure(figsize=(10, 3))
        #     plt.plot(vp_1d)
        #     plt.xlim([0, 100])
        #     plt.xlabel('x-axis')
        #     plt.ylabel('ADU')
        #     plt.title('vertical pattern')

        # horizontal pattern
        # 图像上的横条纹是梳状，分为两个部分，左边大约77个像素是周期小一点，其余的会大一点
        boundary = 77 # boundary between left and width
        boundary_width = 5 # 左右需要平滑过度一下

        y = np.arange(self.bias_shape[0])
        hp_left_param = self.horizontal1_param # 实测数据拟合得到的
        hp_left_1d = hp_left_param[0] * np.sin(TPI * (y / hp_left_param[1] + np.random.rand()))
        hp_left_1d += hp_left_param[2] * np.sin(TPI * (y / hp_left_param[3] + np.random.rand()))
        hp_left_frame = np.broadcast_to(hp_left_1d, [boundary+boundary_width, len(hp_left_1d),]).T

        hp_right_param = self.horizontal2_param
        hp_right_1d = hp_right_param[0] * np.sin(TPI * (y / hp_right_param[1] + np.random.rand()))
        hp_right_1d += hp_right_param[2] * np.sin(TPI * (y / hp_right_param[3] + np.random.rand()))
        hp_right_frame = np.broadcast_to(hp_right_1d, [shape[1] - boundary, len(hp_right_1d)]).T
        
        combine_profile_left = np.ones(boundary + boundary_width)
        combine_profile_left[-boundary_width:] = (boundary_width - np.arange(boundary_width) - 1) / boundary_width
        
        combine_profile_right = np.ones(shape[1] - boundary)
        combine_profile_right[:boundary_width] = np.arange(boundary_width)/boundary_width

        hp_frame = np.zeros(shape)
        hp_frame[:, :boundary+boundary_width] = hp_left_frame * combine_profile_left
        hp_frame[:, boundary:] = hp_right_frame * combine_profile_right

        # residual frame 横条纹外，还有一个垂直方向的梯度，根据测试数据，使用直线加指数函数的方式来拟合
        exp_a, exp_b, lin_a, lin_b = (-1.92377463e+00, 1.32698365e-01, 8.39509583e-04, 4.25384480e-01)
        y_trend = exp_a * np.exp(-(y + 1) * exp_b) + y * lin_a - lin_b

        # random horizontal pattern generated in frequency domain
        # 除了规则横条纹外，还有随机的横条纹，随机的横条纹在频域空间生成，具有相同的频率谱和随机的相位
        rsd_freq_len = len(y_trend) * 4
        red_freq = np.arange(rsd_freq_len)
        red_freq = red_freq - (len(red_freq) - 1)/2
        red_freq = np.exp(-red_freq**2 / 230**2) * 240 + np.random.randn(rsd_freq_len)*30
        red_freq = np.fft.fftshift(red_freq)

        phase = np.random.rand(rsd_freq_len) * TPI
        hp_rsd_1d = np.fft.ifft(np.exp(1j * phase) * red_freq)

        hp_rsd_1d = np.abs(hp_rsd_1d) * self.bias_hp_resd
        hp_rsd_1d = hp_rsd_1d[:rsd_freq_len//4]
        hp_rsd_1d = hp_rsd_1d - np.mean(hp_rsd_1d)

        hp_rsd_1d = y_trend + hp_rsd_1d
        hp_frame = (hp_frame.T + hp_rsd_1d).T

        if not self.switch['bias_hp']:
            hp_frame *= 0

        # if show:  # pragma: no cover
        #     plt.figure(figsize=(10, 3))
        #     # plt.plot(hp_right_1d)
        #     plt.plot(hp_rsd_1d)

        #     # plt.xlim([0, 200])
        #     plt.xlabel('y-axis')
        #     plt.ylabel('ADU')
        #     plt.title('vertical pattern')

        # 接上制冷机后，会有亮暗点
        #cooler interfence effect
        ci_position = 10
        ci_sub_struct = 80
        ci_sub_exp = 2.5
        ci_x_shft = 3
        ci_interval = 250 # 6.25MHz readout / 2.5KHz interfence
        ci_dn = self.cooler_interfence

        npix = shape[0] * shape[1]
        n_ci_event = npix // ci_interval 
        ci_align = np.zeros((n_ci_event, ci_interval))
        ci_align[:, ci_position] = np.random.randn(n_ci_event) * ci_dn
        ci_align[:, ci_position+1] = np.random.randn(n_ci_event) * ci_dn
        yi0 = np.random.randint(0, ci_sub_struct)
        xs0 = (ci_interval - ci_position) / (ci_sub_struct / 2)**ci_sub_exp
        
        for yi in range(n_ci_event):
            sub_yi = (yi - yi0) % ci_sub_struct
            sub_yi = abs(sub_yi - ci_sub_struct / 2)
            shiftx = int(sub_yi**ci_sub_exp * xs0)
            ci_align[yi, :] = np.roll(ci_align[yi, :], shiftx)

        ci_align = np.pad(ci_align.flatten(), (0, npix-n_ci_event*ci_interval))
        ci_frame = ci_align.reshape(shape[0], shape[1])
        
        for yi in range(shape[0]):
            ci_frame[yi, :] = np.roll(ci_frame[yi, :], yi * ci_x_shft)

        if not self.switch['bias_ci']:
            ci_frame *= 0

        bias_shift = 0
        if self.switch['bias_shift']:
            bias_shift = (self.volt - 25) * self.bias_shift_per_volt

        # 混合在一起
        rn_adu = self.readout_noise / self.ph_per_adu
        bias_level = self.bias_level + np.random.randn() * self.bias_level_std
        bias_frame = vp_frame + ci_frame + hp_frame + bias_level
        bias_frame += rn_adu * np.random.randn(shape[0], shape[1]) + bias_shift

        return bias_frame
    
    def nonlinear_effect(self, image):
        """
        nonlinear effect
        """
        fullwell = self.fullwell
        nonlinear_coefficient = self.nonlinear_coefficient
        log.debug(
            f"nonlinear effect added with coefficient {nonlinear_coefficient}")
        image += (image / fullwell)**2 * nonlinear_coefficient * fullwell

        return image
        
    def time_syn(self, t, readout=True, initial=False):
        """
        time synchronization and update the system time and ccd temperature
        
        Parameters
        ----------
        t : float
            relative time
        readout : bool, optional
            If the camera is readout before time synchronization, set readout to True
            if True, the ccd temperature will increase, otherwise, it will decrease
        initial : bool, optional
            If inital is True, the ccd will be intialized to the cooler temperature
        
        """
        if initial:
            self.ccd_temp = self.cooler_temp
            self.system_time = t
            return
        
        dt = np.maximum(t, 0)
        heat = 0
        if readout:
            heat = self.volt * self.heat_speed

        self.ccd_temp = heat + self.cooler_temp + (self.ccd_temp - self.cooler_temp) * np.exp(-dt * self.temper_speed)
        if self.ccd_temp < self.cooler_temp: # 
            self.ccd_temp = self.cooler_temp

        self.system_time += dt

    # def em_cte(self, img):
        
    #     i_shift = 1
    #     cte_coe = 0.1
    #     img_shift_i = np.zeros_like(img)
    #     img_shift_i = np.random.poisson(img * cte_coe)
    #     pass

    def emgain_set(self, em_set, ccd_temp=None, self_update=True):
        """Set emgain from em set value.

        Parameters
        ----------
        em_set : int
            em set value. 3FF is about 1.17×, 200 is about 1000×.
        ccd_temp : float, optional
            CCD temperature. If not given, use the current ccd temperature.
        self_update : bool, optional
            if True, update the emgain and emset. Default is True.
            if False, only return the emgain.
        """
        
        if ccd_temp is None:
            ccd_temp = self.ccd_temp

        volt_coe_a = -0.01828
        volt_coe_b = 43.61

        volt_func = lambda es: volt_coe_a * es + volt_coe_b

        self.volt = volt_func(em_set)

        volt_3ff = volt_func(int('3ff', 16))
        volt_190 = volt_func(int('190', 16))

        em_coe_c = 0.24
        # using the expression of em_b = ln(g199) / constant to make fitting easier
        constant = (np.exp(em_coe_c * volt_190) - np.exp(em_coe_c * volt_3ff))
        # fitting from the ccd test result
        ln_g190 = (-ccd_temp - 7) * 0.0325
        em_coe_b = ln_g190 / constant
        
        emgain = np.exp(em_coe_b * np.exp(em_coe_c * self.volt))
        emgain = np.maximum(1, emgain)
        # print(emgain, em_coe_b, em_coe_c * self.volt, self.volt, np.exp(em_coe_c * self.volt))
        if self_update:
           self.emgain = emgain
           self.emset = em_set

        return emgain
    
    def vertical_blooming(self, image):
        """
        vertical blooming effect
        """
        fullwell = self.fullwell
        line = np.arange(image.shape[0])
        yp, xp = np.where(image > fullwell)
        n_saturated = len(xp)
        log.debug(f"{len(xp)} pixels are saturated!")
        if n_saturated > 5000:
            log.warning(f"More than 5000({len(xp)}) pixels are saturated!")
        img0 = image.copy()
        for x, y in zip(xp, yp):
            image[:, x] += np.exp(-(line-y)**2/20**2) * img0[y, x] * 0.2
        return np.minimum(image, fullwell)

    def emregester_blooming(self, image, max_iteration=5):
        """
        emregester blooming effect
        """
        line = image.flatten().copy()

        curve_x = np.arange(1300)+2
        curve_y = np.exp(11*curve_x**(-0.19)-11)
        curve_y[0] = 0
        curve_y /= curve_y.sum()

        over_limit_coe = 0.999

        saturated = image > self.em_fullwell
        n_saturated = saturated.sum()
        if n_saturated > 0:
            log.debug(f"{n_saturated} pixels are saturated during EM process.")

        if n_saturated > 2000:
            log.warning(
                f"More than 2000 ({n_saturated}) pixels are saturated during EM process!")

        for index in range(max_iteration):
            over_limit = np.maximum(
                line - self.em_fullwell * over_limit_coe, 0)
            line = np.minimum(line, self.em_fullwell * over_limit_coe)
            blooming = np.convolve(over_limit, curve_y, mode='full')[
                :len(line)]
            line = line + blooming
            n_over = (line > self.em_fullwell).sum()
            if n_over <= 0:
                break

            log.debug(
                f'{index}/{max_iteration} loop: saturated pixel number: {n_over}')
        line = np.minimum(line, self.em_fullwell)
        return line.reshape(image.shape)

    # def add_em_cte_conv(self, image):
    #     shape = image.shape
    #     img_line = np.convolve(image.flatten(), self.cti_trail, mode='full')
    #     return img_line[:shape[0]*shape[1]].reshape(shape)

    def readout(self, image_focal, em_set, expt_set, image_cosmic_ray=False, emgain=None):
        """From focal planet image to ccd output.
        Interface function for emccd. Simulate the readout process.
        
        Parameters
        ----------
        image_focal : np.ndarray
            focal planet image. Unit: photon-electron/s/pixel
        em_set : int
            emgain set value.
        expt_set: float
            exposure time set value. Unit: s (0 mains 1ms)
        image_cosmic_ray : np.ndarray, optional
            cosmic ray image. Unit: photon-electron/pixel
        emgain: float, optional
            if not None, use the given emgain. Else, calculate the emgain using em_set

        Returns
        -------
        np.ndarray
            ccd output image. Unit: ADU
        """

        expt = expt_set
        if expt_set == 0:
            expt = 0.001
        
        dt = self.readout_time + expt
        self.time_syn(dt, readout=True)

        if emgain is None:
            self.emgain_set(em_set, self.ccd_temp)
        else:
            self.emgain = emgain
            self.emset = 0

        emgain = self.emgain

        image = image_focal.astype(float)
        log.debug(f"image total photon: {image.sum()}")
        if self.switch['flat']:
            image = image * self.flat

        img_dark = np.zeros(self.dark_shape)
        img_dark[
            self.bdark:self.plszy+self.bdark,
            self.ldark:self.ldark+self.plszx
        ] = image
        image = img_dark

        if self.switch['dark']:
            image = image + self.dark

        image *= expt

        if image_cosmic_ray is not None and self.switch['cosmicray']:
            image += image_cosmic_ray

        if self.switch['nonlinear']:
            image = self.nonlinear_effect(image)

        if self.switch['blooming']:
            image = self.vertical_blooming(image)

        if self.switch['badcolumn']:
            for i in range(self.bad_col.shape[1]):
                deadpix_x = self.bad_col[0, i]
                deadpix_y = self.bad_col[1, i]
                image[deadpix_y:, deadpix_x] = 0

        img_bias = np.zeros(self.bias_shape, dtype=int)
        
        img_bias[
            self.pscan2: -self.oscan2,
            self.pscan1: -self.oscan1
        ] = np.random.poisson(image)
        image = img_bias

        if self.switch['shutter']:
            img_line = image_focal.sum(axis=0)
            image_shutter = np.broadcast_to(img_line, [self.bias_shape[0], self.flat_shape[1]])
            image_shutter = image_shutter / self.flat_shape[1] * self.shift_time
            image_shutter = np.random.poisson(image_shutter)
            image[:, self.pscan1+self.ldark:-self.oscan1-self.rdark] += image_shutter
        
        
        if self.switch['cic']:
            cic_frame = np.zeros((self.dark_shape[0], self.bias_shape[1])) + self.cic
            image[self.pscan2:-self.oscan2, :] += np.random.poisson(cic_frame)

        # em_fix found to fitting the gamma distribution to theoritical one.
        # here is the theoritical distribution. See Robbins and Hadwen 2003 for more details.
        # >>> pEM = np.exp(np.log(emgain)/self.emreg_num) - 1
        # >>> for _ in range(self.emreg_num):
        # >>>     image += np.random.binomial(image, pEM)
        # This code is too slow, so we used a modified gamma
        
        em_fix = self.em_fix_fuc_fit(emgain) * emgain
        image = np.random.gamma(image, em_fix) + image * (emgain - em_fix)

        
        if self.switch['em_blooming']:
            image = self.emregester_blooming(image)

        image = np.maximum(image, 0)
        image = image.astype(int)

        if self.switch['cte']:
            big_cte = self.em_cte ** (self.emreg_num / self.emreg_cal_num)
            for _ in range(self.emreg_cal_num):
                residual = np.random.binomial(image, 1-big_cte)
                image[:, 1:] += residual[:, :-1] - residual[:, 1:]

        bias = self.bias_frame()

        image = image / self.ph_per_adu + bias
        image = np.minimum(image, self.max_adu)
        image = np.maximum(image, 0)

        log.debug(
            f"emccd paramters: \
            emset: {em_set} \
            emgain: {emgain} \
            expt: {expt} \
            readout time: {dt}"
        )

        return image.astype(dtype=np.uint16)


# if __name__ == '__main__':
    # import matplotlib.pyplot as plt
    # emccd = EMCCD()
    # image_focal = np.zeros((emccd.plszy, emccd.plszx)) + 1000
    # image_focal[100:105, 100:105] = 10_000_000
    # after_cte = emccd.emregester_blooming(image_focal, max_iteration=100)
    # print(after_cte.sum(), image_focal.sum())

    # fits.writeto('after_cte.fits', after_cte, overwrite=True)

#     # darksz_x = emccd.plszx + emccd.rdark + emccd.ldark
#     # darksz_y = emccd.plszy + emccd.udark + emccd.bdark
#     # iamge_cosmic_ray = np.zeros((darksz_y, darksz_x))
#     # emgain = 10
#     # expt = 10
#     # image = emccd.readout(image_focal, emgain, expt, iamge_cosmic_ray)
#     # fits.writeto('test.fits', image, overwrite=True)

    # image = np.zeros((1000, 1000))
    # make_cosmic_ray_frame = CosmicRayFrameMaker()
    # crimage = make_cosmic_ray_frame(image.shape, 3000)
    # fits.writeto('crimage.fits', crimage, overwrite=True)

#     # emccd.add_stripe_effect(image)