diff --git a/CpicImgSim/camera.py b/CpicImgSim/camera.py
new file mode 100644
index 0000000000000000000000000000000000000000..d5ada923b052d6d7dbcbec2e76a7e97f9ad4d82e
--- /dev/null
+++ b/CpicImgSim/camera.py
@@ -0,0 +1,633 @@
+import os
+import yaml
+import numpy as np
+import scipy.ndimage as nd
+from astropy.io import fits
+
+from .config import cpism_refdata, solar_spectrum, MAG_SYSTEM
+from .utils import region_replace, random_seed_select
+from .io import log
+from .optics import filter_throughput
+
+
+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 EMCCD(object):
+    """
+    EMCCD camera class
+
+    Parameters
+    ----------
+    config_file : str
+        config file name
+
+    Attributes
+    ----------
+    switch : dict
+        switch for each camera effects, including:
+            - 'flat': bool,
+            - 'dark': bool,
+            - 'stripe': bool,
+            - 'cic': bool,
+            - 'cte': bool,
+            - 'badcolumn': bool,
+            - 'nonlinear': bool,
+            - 'cosmicray': bool,
+            - 'blooming': bool,
+
+    """
+
+    def __init__(self, config_file="emccd_config.yaml"):
+        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.fullwell = 80_000
+        self.em_fullwell = 780_000
+
+        # if config file exists, load it, otherwise use default values
+        config_file = cpism_refdata + '/camera/' + config_file
+        log.debug(f"emccd config file: {config_file}")
+        if os.path.exists(config_file):
+            self.load_config(config_file)
+
+        else:  # pragma: no cover
+            # set default values for EMCCD
+            # note: these values are default values, you can change them by load_config()
+            # ↓↓↓↓↓↓↓start default values setting↓↓↓↓↓↓
+
+            self.readout_noise = 40
+            self.ph_per_adu = 8
+            self.bias_level = 30
+            self.max_adu = 16_383
+
+            self.switch = {
+                'flat': True,
+                'dark': True,
+                'stripe': True,
+                'cic': False,
+                'cte': False,
+                'badcolumn': True,
+                'nonlinear': False,
+                'cosmicray': True,
+                'blooming': False,
+            }
+
+            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/emcid_cic.fits'
+            self.bad_col_file = cpism_refdata + '/camera/emccd_bad_columns.fits'
+            # ↑↑↑↑↑↑↑end default values setting↑↑↑↑↑↑
+            # note: these values are default values, you can change them by load_config()
+
+        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.image_shape = [biassz_y, biassz_x]
+
+        self.flat = fits.getdata(self.flat_file)
+        self.cic = fits.getdata(self.cic_file)
+        self.dark = fits.getdata(self.dark_file)
+        self.bad_col = fits.getdata(self.bad_col_file)
+
+    def load_config(self, config_file):
+        """
+        load config file. Only for internal use.
+        """
+
+        with open(config_file, 'r') as f:
+            config = yaml.load(f, Loader=yaml.FullLoader)
+
+        self.switch = config['switch']
+
+        self.readout_noise = config['readout_noise']
+        self.ph_per_adu = config['ph_per_adu']
+        self.bias_level = config['bias_level']
+        self.max_adu = config['max_adu']
+
+        self.dark_file = cpism_refdata + "/camera/" + config['dark_file']
+        self.flat_file = cpism_refdata + "/camera/" + config['flat_file']
+        self.cic_file = cpism_refdata + "/camera/" + config['cic_file']
+        self.bad_col_file = cpism_refdata + \
+            "/camera/" + config['bad_col_file']
+
+    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 nonlinear_effect(self, image):
+        """
+        nonlinear effect
+        """
+        fullwell = self.fullwell
+        nonlinear_coefficient = 0.1
+        log.debug(
+            f"nonlinear effect added with coefficient {nonlinear_coefficient}")
+        image += (image / fullwell)**2 * nonlinear_coefficient * fullwell
+
+        return image
+
+    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}')
+
+        return line.reshape(image.shape)
+
+    def cte(self, image):
+        """
+        cte effect
+        """
+        image = self.emregester_blooming(image, max_iteration=5)
+        return image
+
+    def readout(self, image_focal, emgain, expt, image_cosmic_ray=None):
+        """
+        emccd readout. Generate a image with emccd readout effect.
+
+        Parameters
+        ----------
+        image_focal : numpy.ndarray
+            image at focal plane. Unit: electron/second
+        emgain : float
+            emgain of emccd
+        expt : float
+            exposure time. Unit: second
+        image_cosmic_ray : numpy.ndarray, optional
+            cosmic ray image. Unit: electron/second, by default None
+
+        Returns
+        -------
+        numpy.ndarray
+            image with emccd readout effect. Unit: ADU
+
+        Notes
+        -----
+        1. effects include: dark, flat, cte, blooming, nonlinear, etc. Can be turned on/off by switch.
+        2. size of input image_focal must be 1024x1024
+        3. size of output image is 1080x1056 (including overscan and dark reference region)
+        4. Q.E is not included in this function. It should be included in image_focal. See optics.py for details.
+        """
+        log.debug(
+            fr"EMCCD readout: {emgain=:}, {expt=:}, image_comic_ray:{'None' if image_cosmic_ray is None else 'Not None'}")
+
+        log.debug(
+            f"camera effects  switch={self.switch}"
+        )
+        image = image_focal * expt
+        if self.switch['flat']:
+            image = image * self.flat
+
+        if self.switch['nonlinear']:
+            image = self.nonlinear_effect(image)
+
+        darksz_x = self.plszx + self.rdark + self.ldark
+        darksz_y = self.plszy + self.udark + self.bdark
+        img_dark = np.zeros((darksz_y, darksz_x))
+        img_dark[
+            self.bdark:self.plszy+self.bdark,
+            self.ldark:self.ldark+self.plszx
+        ] = image
+        image = img_dark
+
+        if self.switch['dark']:
+            image += self.dark * expt
+
+        if self.switch['cic']:
+            image += self.cic
+
+        if image_cosmic_ray is not None and self.switch['cosmicray']:
+            image += image_cosmic_ray
+
+        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
+
+        biassz_x = darksz_x + self.pscan1 + self.oscan1
+        biassz_y = darksz_y + self.pscan2 + self.oscan2
+        img_bias = np.zeros((biassz_y, biassz_x), dtype=int)
+
+        seed = random_seed_select()
+        log.debug(f"photon noise seed: {seed}")
+        img_bias[
+            self.pscan2:self.pscan2+darksz_y,
+            self.pscan1:self.pscan1+darksz_x
+        ] = np.random.poisson(image)
+        image = img_bias
+
+        if self.switch['cte']:
+            image = self.cte(image * emgain) / emgain
+
+        seed = random_seed_select()
+        log.debug(f"gamma noise seed: {seed}")
+        if emgain != 1:
+            image = np.random.gamma(image, emgain)
+
+        image = np.minimum(image, self.em_fullwell)
+
+        seed = random_seed_select()
+        log.debug(f"readout noise seed: {seed}")
+        image += np.random.randn(biassz_y, biassz_x) * self.readout_noise
+        image = image / self.ph_per_adu + self.bias_level
+
+        if self.switch['stripe']:
+            image += self.add_stripe_effect(image)
+
+        image = np.minimum(image, self.max_adu)
+        image = np.maximum(image, 0)
+
+        return image.astype(np.uint16)
+
+    def add_stripe_effect(self, image):
+        """
+        add stripe effect
+        """
+        shape = image.shape
+
+        v_width = 1
+        v_amplitude = 30
+        v_limit = 0.01
+        v_base = 10
+
+        h_width = 20
+        h_amplitude = 10
+        h_limit = 0.9
+        h_base = 20
+
+        index = np.linspace(0, np.pi, shape[0] * shape[1])
+
+        def stripe(width, limit, amplitude, base, axis=0):
+            seed = random_seed_select()
+            log.debug(f"stripe noise seed: {seed}")
+            dim_axis = shape[axis]
+            dim_other = shape[0] * shape[1] // shape[axis]
+            value = np.sin(index / width * dim_axis + np.pi *
+                           dim_axis / width * np.random.randint(1024))
+
+            value = np.maximum(value, -limit)
+            value = np.minimum(value, limit)
+            value = (value / limit + limit) / 2 * amplitude + base
+            value = value.reshape(dim_axis, dim_other)
+
+            if axis == 1:
+                value = value.T
+
+            return value
+
+        output = stripe(v_width, v_limit, v_amplitude, v_base, axis=1)
+        output += stripe(h_width, h_limit, h_amplitude, h_base, axis=0)
+        return output
+
+        # # plt.plot(horizontal_index, horizontal_value)
+        # # # plt.xlim([0, 6.28])
+        # # plt.show()
+
+        # fits.writeto('horizontal_value.fits', output, overwrite=True)
+
+
+# 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)