V2.0 developing initial

.gitignore

+# Build and Release Folders
+bin-debug/
+bin-release/
+[Oo]bj/
+[Bb]in/
+build/
+dist/
+starmodel/
+cpism_refdata/
+*.egg-info
+
+# Other files and folders
+.settings/
+_*/
+
+~*
+.VSCodeCounter
+.vscode
+*.log
+output/
+*.log.*
+/*.fits
+
+docs/notebooks/image_files/_*
+
+# unitest output
+tests/.coverage
+tests/htmlcov/
+
+# Executables
+*.swf
+*.air
+*.ipa
+*.apk
+
+# Project files, i.e. `.project`, `.actionScriptProperties` and `.flexProperties`
+# should NOT be excluded as they contain compiler settings and other important
+# information for Eclipse / Flash Builder.

CpicImgSim/__init__.py

+# from .main import quick_run, observation_simulation
+from .optics import make_focus_image, focal_mask
+from .target import star_photlam, planet_contrast, extract_target_x_y, spectrum_generator
+from .camera import EMCCD, CosmicRayFrameMaker, sky_frame_maker
+from .config import __version__
+
+__all__ = [
+    "EMCCD",
+    "CosmicRayFrameMaker",
+    "sky_frame_maker",
+    "star_photlam",
+    "planet_contrast",
+    "extract_target_x_y",
+    "spectrum_generator",
+    "make_focus_image",
+    "focal_mask",
+    # "quick_run",
+    # "observation_simulation"
+]
\ No newline at end of file

CpicImgSim/camera.py

+import os
+import yaml
+import numpy as np
+import scipy.ndimage as nd
+from astropy.io import fits
+import matplotlib.pyplot as plt
+import math
+
+from CpicImgSim.config import cpism_refdata, solar_spectrum, MAG_SYSTEM
+from CpicImgSim.utils import region_replace, random_seed_select
+from CpicImgSim.io import log
+from CpicImgSim.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 CpicVisEmccd(object):
+    def __init__(self):
+        self._defaut_config()
+        self.config_init()
+
+    @classmethod
+    def from_config(cls, config):
+        obj = cls.__new__()
+        obj._defaut_config()
+        obj.__dict__.update(config)
+        obj._refresh_config()
+        return obj
+
+    def _defaut_config(self):
+        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 = {
+            'flat': True,
+            'dark': True,
+            'stripe': True,
+            'cic': True,
+            'cte': True,
+            'badcolumn': True,
+            'nonlinear': True,
+            'cosmicray': True,
+            'blooming': False,
+            'bias_vp': True,
+            'bias_hp': True,
+            'bias_ci': True,
+            'bias_shift': 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_cic.fits'
+        self.bad_col_file = cpism_refdata + '/camera/emccd_bad_columns.fits'
+        self.cic = 0.2
+        
+        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 = 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
+
+
+    def shutter_effect(self, image, expt):
+        image_mean = image.mean(axis=0)
+        image_shutter = np.broadcast_to(image_mean, image.shape)
+        image_shutter /= image_shutter.mean() 
+        image_shutter *= image.mean() * self.shift_time / expt
+        return image_shutter.astype(int)
+        
+    
+    def config_init(self):
+        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.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)
+
+        self.ccd_temp = self.cooler_temp
+        self.system_time = 0
+        self.time_syn(0, initial=True)
+        self.emgain_set(1024, -30)
+
+        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 bias_frame(self, show=False):
+        shape = self.image_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.image_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: 
+            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 = 0.1
+        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):
+        if initial:
+            self.ccd_temp = self.cooler_temp
+            self.system_time = t
+            return
+        
+        dt = t
+        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):
+        
+        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 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):
+        expt = expt_set
+        if expt_set == 0:
+            expt = 0.001
+        
+        dt = self.readout_time + expt
+        self.time_syn(dt, readout=True)
+
+        emgain = self.emgain_set(em_set, self.ccd_temp)
+
+        image = image_focal
+        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
+
+        image *= expt
+
+        if self.switch['nonlinear']:
+            image = self.nonlinear_effect(image)
+
+        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)
+
+        img_bias[
+            self.pscan2:self.pscan2+darksz_y,
+            self.pscan1:self.pscan1+darksz_x
+        ] = np.random.poisson(image)
+
+        sub_img = image[100:-100, 100:-100]
+        print(np.random.poisson(sub_img).std(), sub_img.mean())
+
+        image = img_bias
+
+        lower_limit_n = int(np.log(emgain)/np.log(1 + 0.2)) # to make pEM < 0.2 
+        emregister_num = np.maximum(lower_limit_n, self.emreg_cal_num)  
+        emregister_num = np.minimum(emregister_num, self.emreg_num)
+        
+        pEM = np.exp(np.log(emgain)/emregister_num) - 1
+        big_cte = self.em_cte ** (self.emreg_num / emregister_num)
+
+        image += np.random.poisson(image ) 
+
+        # for _ in range(emregister_num):
+        #     # image += np.random.binomial(image, pEM)
+        #     residual = np.random.binomial(image, 1-big_cte)
+        #     image[:, 1:] += residual[:, :-1] - residual[:, 1:]
+
+
+        sub_img = image[100:-100, 100:-100] / emgain
+        print(sub_img.std(), sub_img.mean(), emgain)
+
+        bias = self.bias_frame()
+
+        image = image / self.ph_per_adu + bias
+        image = np.minimum(image, self.max_adu)
+        image = np.maximum(image, 0)
+
+        return image#.astype(dtype=np.uint16)
+        
+
+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="ccd201_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:
+            raise(ValueError('config_file error'))
+        
+        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, 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 
+        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
+        image *= 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}")
+        print(image.mean())
+        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)

CpicImgSim/config.py

+import os, yaml
+import warnings
+
+cpism_refdata = os.environ.get('cpism_refdata', './cpism_refdata')
+
+if not os.path.exists(cpism_refdata):  # pragma: no cover
+    raise Exception(
+        "Can not find CPISM reference data.")
+
+
+if not os.path.exists(os.environ.get('PYSYN_CDBS', './trd')):  # pragma: no cover
+    raise Exception(
+        "Can not find PYSYN Stellar reference data.")
+
+# we need to ignore the warning from pysynphot, because we only use the star models.
+with warnings.catch_warnings():  # pragma: no cover
+    warnings.filterwarnings("ignore")
+    import pysynphot as S
+
+solar_spectrum = S.FileSpectrum(
+    f"{os.environ['PYSYN_CDBS']}/grid/solsys/solar_spec.fits")
+solar_spectrum.convert('photlam')
+
+config_file = cpism_refdata + '/optics/optics_config.yaml'
+if not os.path.exists(config_file):  # pragma: no cover
+    raise FileNotFoundError(f"光学配置文件不存在({config_file})")
+with open(config_file, 'r') as f:
+    optics_config = yaml.load(f, Loader=yaml.FullLoader)
+
+MAG_SYSTEM = 'abmag'
+
+__version__ = '2.0.0'
+
+def which_focalplane(band):
+    """
+    Return the name of the focalplane which the band belongs to.
+
+    Parameters
+    -----------
+    band: str
+        The name of the band.
+        from ['f565', 'f661', 'f743', 'f883', 'f940', 'f1265', 'f1425', 'f1542', 'wfs']
+
+    Returns
+    --------
+    str
+        The name of the focalplane.
+        'vis' or 'nir' or 'wfs'
+
+    Raises
+    -------
+    ValueError
+        If the band is not in ['f565', 'f661', 'f743', 'f883', 'f940', 'f1265', 'f1425', 'f1542', 'wfs']
+    """
+    band = band.lower()
+    if band in ['f565', 'f661', 'f743', 'f883']:
+        return 'vis'
+    if band in ['f940', 'f1265', 'f1425', 'f1542']:
+        return 'nir'
+    if band in ['wfs']:
+        return 'wfs'
+
+    raise ValueError(f"未知的波段{band}")
\ No newline at end of file

CpicImgSim/io.py

+import numpy as np
+from astropy.io import fits
+import yaml
+from CpicImgSim.config import cpism_refdata, __version__, which_focalplane
+from CpicImgSim.utils import Logger
+
+import os
+from datetime import datetime, timedelta
+from astropy.coordinates import SkyCoord
+import re
+import json
+import pandas as pd
+
+
+config_file = f'{cpism_refdata}/cpism_config.yaml'
+with open(config_file, 'r') as fid:
+    config = yaml.load(fid, Loader=yaml.FullLoader)
+
+output_dir = config['output_dir']
+if config['relative_path']:
+    ref_dir_base = os.path.dirname(cpism_refdata)
+    output_dir = f'{ref_dir_base}/{output_dir}'
+
+
+log_dir = output_dir + '/LOG'
+tmp_dir = config['tmp_dir']
+log_level = config['log_level']
+header_check = config['check_fits_header']
+
+for dir in ['', 'TMP', 'CAL', 'SCI', 'LOG']:
+    sub_dir = f"{output_dir}/{dir}"
+    if not os.path.exists(sub_dir):
+        os.makedirs(sub_dir)
+
+tmp_folder_path = '.'
+if tmp_dir == 'TMP':
+    tmp_folder_path = output_dir + '/TMP'
+
+log = Logger(log_dir+'/cpism_pack.log', log_level).logger
+
+
+def check_and_update_fits_header(header):
+    """
+    Check the header keywords and update the description according to the data model.
+
+    Parameters
+    -----------
+    header: astropy.io.fits.header.Header
+        The header to be checked.
+
+    Returns
+    --------
+    None
+    """
+
+    hdu = 'image'
+    if 'FILETYPE' in header.keys():
+        hdu = 'primary'
+
+    model_file = f"{cpism_refdata}/io/image_header.json"
+    if hdu == 'primary':
+        model_file = f"{cpism_refdata}/io/primary_header.json"
+
+    with open(model_file, 'r', encoding='utf-8') as fid:
+        data_model = json.load(fid)
+
+    dm_comment = {}
+
+    def print_warning(info):
+        if header_check:
+            log.warning(info)
+
+    # check existance and format of keyword in fits header
+    for keyword, comment, format, _, _ in data_model:
+        if pd.isnull(comment):
+            comment = ''
+
+        if len(comment) > 46:
+            comment = comment[:46]
+            print_warning(
+                f"Keyword {keyword} has a comment longer than 46 characters. It will be truncated to 46 characters.")
+
+        dm_comment[keyword] = comment
+
+        if keyword not in header.keys():
+            print_warning(f"Keyword {keyword} not found in [{hdu}] header.")
+
+        elif not pd.isnull(format):
+            value = header[keyword]
+            # check the type of the value, I for int, R for float, C for str
+            if isinstance(value, str):
+                type = 'C'
+            elif isinstance(value, float):
+                type = 'R'
+            elif isinstance(value, bool):
+                type = 'L'
+            elif isinstance(value, int):
+                type = 'I'
+            else:
+                type = 'U'
+
+            if type != format[0]:
+                print_warning(
+                    f"Keyword {keyword} has wrong type in [{hdu}]. {format[0]} expected, {type} found.")
+
+    # check if there are extral keyword in header, and update the comment
+    for keyword in header.keys():
+        # print(keyword)
+        if keyword not in dm_comment.keys():
+            print_warning(
+                f"Keyword {keyword} not found in the [{hdu}] data model.")
+        else:
+            header[keyword] = (header[keyword], dm_comment[keyword])
+
+    return header
+
+
+def obsid_parser(
+        obsid: int):
+    """
+    Parse the obsid to get the obstype.
+
+    Parameters
+    ----------
+    obsid: str
+        The obsid of the observation.
+        Obsid must be 11 digits and start with 5 for CPIC.
+
+    Returns
+    -------
+    str
+        The obstype of the observation.
+
+    Raises
+    ------
+    ValueError
+        If the obsid is not 11 digits or does not start with 5.
+    """
+    obsid = str(obsid)
+    if len(obsid) != 11:
+        raise ValueError('Obsid must be 11 digits.')
+
+    if obsid[0] != '5':
+        raise ValueError('Obsid must start with 5 for CPIC')
+
+    obstype_dict = {
+        '00': 'BIAS',
+        '01': 'DARK',
+        '02': 'FLAT',
+        '03': 'BKGD',
+        '04': 'LASR',
+        '10': 'SCIE',
+        '11': 'DENF',
+        '12': 'CALS',
+        '15': 'TEMP'
+    }
+    obstype = obstype_dict.get(obsid[1:3], 'DEFT')
+    return obstype
+
+
+def datetime_obj_to_mjd(time_obj):
+    """
+    transfer datetime object to mean julian date (MJD).
+
+    Parameters
+    ----------
+    time_obj: datetime.datetime
+        The datetime object.
+
+    Returns
+    -------
+    float
+        The mean julian date (MJD).
+    """
+    return (time_obj - datetime(1858, 11, 17)).total_seconds() / 86400
+
+
+def primary_hdu(
+        obs_info: dict,
+        gnc_info: dict,
+        filename_output=False):
+    """
+    Generate the primary hdu of the fits file.
+
+    Parameters
+    ----------
+    obs_info: dict
+        The parameters of the observation. See `save_fits` function.
+    gnc_info: dict
+        The gnc information of the observation.
+    filename_output: bool (optional)
+        If True, the folder and the filename will be returned.
+
+    Returns
+    -------
+    fits.PrimaryHDU
+        The primary hdu of the fits file.
+    str, str
+        The folder and filename of the fits file. Only returned if filename_output is True.
+
+    Notes
+    -----
+    The gnc_info dict should contain the information of orbit and observation.
+    these informations are used to genrated the hdu header. Refer to the data model for more information.
+    """
+
+    camera_config, _ = load_camera_and_optics_config(obs_info['band'])
+
+    obsid = obs_info['obsid']
+
+    exp_start = gnc_info.get(
+        'EXPSTART', datetime.now().isoformat(timespec='seconds'))
+    exp_start = datetime.fromisoformat(exp_start)
+
+    duartion = (obs_info['expt'] +
+                camera_config['readout_time']) * obs_info['nframe']
+    default_end = exp_start + timedelta(seconds=duartion)
+
+    exp_end = gnc_info.get('EXPEND', default_end.isoformat(timespec='seconds'))
+    exp_end = datetime.fromisoformat(exp_end)
+
+    filename = "CSST_CPIC"
+    filename += "_" + which_focalplane(obs_info['band']).upper()
+    filename += "_" + obsid_parser(obsid)
+    filename += "_" + exp_start.strftime("%Y%m%d%H%M%S")
+    filename += "_" + exp_end.strftime("%Y%m%d%H%M%S")
+    filename += f"_{obsid}_X_L0_V01.fits"
+
+    type_dir = 'CAL'
+    if str(obsid)[1] == '1':
+        type_dir = 'SCI'
+
+    mjd_dir = f"{datetime_obj_to_mjd(exp_start):.0f}"
+    folder = f"{type_dir}/{mjd_dir}"
+
+    header = fits.Header()
+    # General keywords
+    header['SIMPLE'] = True
+    header['BITPIX'] = 8
+    header['NAXIS'] = 0
+    header['EXTEND'] = True
+    header['NEXTEND'] = 1  # + parameters['nframe']
+    header['GROUPS'] = False
+
+    header['DATE'] = datetime.now().isoformat(timespec='seconds')
+    heaer_filename = filename[:-4]
+    if len(heaer_filename) > 68:
+        heaer_filename = heaer_filename[:68]
+    header['FILENAME'] = heaer_filename
+    header['FILETYPE'] = obsid_parser(obsid)
+    header['TELESCOP'] = 'CSST'
+    header['INSTRUME'] = 'CPIC'
+    header['RADECSYS'] = 'ICRS'
+    header['EQUINOX'] = 2000.0
+    header['FITSCREA'] = f'CPISM V{__version__}'
+
+    cstar = {'ra': '0d', 'dec': '0d'}
+    if obs_info['target'] != {}:
+        cstar = obs_info['target']['cstar']
+
+    radec = SkyCoord(cstar['ra'], cstar['dec'])
+    target_name = radec.to_string('hmsdms')
+    target_name = re.sub(R'[hdms\s]', '', target_name)
+    header['OBJECT'] = cstar.get('name', target_name)
+    header['TARGET'] = target_name
+    header['OBSID'] = str(obsid)
+    header['OBJ_RA'] = radec.ra.degree
+    header['OBJ_DEC'] = radec.dec.degree
+
+    # telescope information
+    header['REFFRAME'] = 'CSSTGSC-1.0'
+    header['DATE-OBS'] = exp_start.isoformat(timespec='seconds')
+    header['SATESWV'] = 'SIMULATION'
+
+    header['EXPSTART'] = datetime_obj_to_mjd(exp_start)
+    header['CABSTART'] = gnc_info.get('CABSTART', header['EXPSTART'])
+    header['SUNANGL0'] = gnc_info.get('SUNANGL0', -1.0)
+    header['MOONANG0'] = gnc_info.get('MOONANG0', -1.0)
+    header['TEL_ALT0'] = gnc_info.get('TEL_ALT0', -1.0)
+    header['POS_ANG0'] = gnc_info.get(
+        'POS_ANG0', float(obs_info['rotation']))
+    header['POSI0_X'] = gnc_info.get('POSI0_X', -1.0)
+    header['POSI0_Y'] = gnc_info.get('POSI0_Y', -1.0)
+    header['POSI0_Z'] = gnc_info.get('POSI0_Z', -1.0)
+    header['VELO0_X'] = gnc_info.get('VELO0_X', -1.0)
+    header['VELO0_Y'] = gnc_info.get('VELO0_Y', -1.0)
+    header['VELO0_Z'] = gnc_info.get('VELO0_Z', -1.0)
+    header['EULER0_1'] = gnc_info.get('EULER0_1', -1.0)
+    header['EULER0_2'] = gnc_info.get('EULER0_2', -1.0)
+    header['EULER0_3'] = gnc_info.get('EULER0_3', -1.0)
+    header['RA_PNT0'] = gnc_info.get('RA_PNT0', header['OBJ_RA'])
+    header['DEC_PNT0'] = gnc_info.get('DEC_PNT0', header['OBJ_DEC'])
+
+    header['EXPEND'] = datetime_obj_to_mjd(exp_end)
+    header['CABEND'] = gnc_info.get('CABEDN', header['EXPEND'])
+    header['SUNANGL1'] = gnc_info.get('SUNANGL1', header['SUNANGL0'])
+    header['MOONANG1'] = gnc_info.get('MOONANG1', header['MOONANG0'])
+    header['TEL_ALT1'] = gnc_info.get('TEL_ALT1', header['TEL_ALT0'])
+    header['POS_ANG1'] = gnc_info.get('POS_ANG1', header['POS_ANG0'])
+    header['POSI1_X'] = gnc_info.get('POSI1_X', header['POSI0_x'])
+    header['POSI1_Y'] = gnc_info.get('POSI1_Y', header['POSI0_y'])
+    header['POSI1_Z'] = gnc_info.get('POSI1_Z', header['POSI0_z'])
+    header['VELO1_X'] = gnc_info.get('VELO1_X', header['VELO0_x'])
+    header['VELO1_Y'] = gnc_info.get('VELO1_Y', header['VELO0_y'])
+    header['VELO1_Z'] = gnc_info.get('VELO1_Z', header['VELO0_z'])
+    header['EULER1_1'] = gnc_info.get('EULER1_1', header['EULER0_1'])
+    header['EULER1_2'] = gnc_info.get('EULER1_2', header['EULER0_2'])
+    header['EULER1_3'] = gnc_info.get('EULER1_3', header['EULER0_3'])
+    header['RA_PNT1'] = gnc_info.get('RA_PNT1', header['RA_PNT0'])
+    header['DEC_PNT1'] = gnc_info.get('DEC_PNT1', header['DEC_PNT0'])
+
+    header['EXPTIME'] = (exp_end - exp_start).total_seconds()
+    header['EPOCH'] = float(exp_start.year)
+
+    header['CHECKSUM'] = '0000000000000000'
+    header['DATASUM'] = '0000000000'
+
+    check_and_update_fits_header(header)
+
+    # other information
+    hdu = fits.PrimaryHDU(header=header)
+    hdu.add_checksum()
+
+    if filename_output:
+        return hdu, folder, filename
+    else:
+        return hdu
+
+
+def load_camera_and_optics_config(band):
+    """
+    Load camera and optics configuration from reference data.
+
+    Parameters
+    ----------
+    band : str
+        Band name.
+
+    Returns camera_config, optics_config : dict, dict
+    """
+    camera = which_focalplane(band)
+
+    if camera == 'vis':
+        config_file = 'emccd_config.yaml'
+    elif camera == 'nir':
+        raise ValueError('NIR camera is not supported yet')
+        config_file = 'nir_config.yaml'
+
+    with open(f"{cpism_refdata}/camera/{config_file}", 'r') as fid:
+        camera_config = yaml.load(fid, Loader=yaml.FullLoader)
+
+    with open(f"{cpism_refdata}/optics/optics_config.yaml", 'r') as fid:
+        optics_config = yaml.load(fid, Loader=yaml.FullLoader)[camera]
+
+    return camera_config, optics_config
+
+
+def frame_header(obs_info, index, bunch_start, primary_header={}):
+    """
+    Generate the header for a single frame.
+
+    Parameters
+    ----------
+    obs_info : dict
+        Dictionary of parameters. See `save_fits` function.
+    index : int
+        Frame index.
+    bunch_start : str
+        Start time of the bunch.
+    primary_header : dict (optional)
+        Primary header. default: {}
+
+    Returns
+    -------
+    astropy.io.fits.Header
+
+    """
+
+    header = fits.Header()
+    camera_config, optics_config = load_camera_and_optics_config(
+        obs_info['band'])
+
+    plszx = camera_config['plszx']
+    plszy = camera_config['plszy']
+    pscan1 = camera_config['pscan1']
+    pscan2 = camera_config['pscan2']
+    oscan1 = camera_config['oscan1']
+    oscan2 = camera_config['oscan2']
+    udark = camera_config['udark']
+    bdark = camera_config['bdark']
+    ldark = camera_config['ldark']
+    rdark = camera_config['rdark']
+    imgszx = plszx + pscan1 + oscan1 + ldark + rdark
+    imgszy = plszy + pscan2 + oscan2 + udark + bdark
+
+    header['XTENSION'] = 'IMAGE'
+    header['BITPIX'] = 16
+    header['NAXIS'] = 2
+    header['NAXIS1'] = 1080
+    header['NAXIS2'] = 1050
+    header['EXTNAME'] = 'IMAGE'
+    header['EXTVER'] = 1
+    header['BSCALE'] = 1.0
+    header['BZERO'] = 32768.0
+    header['BUNIT'] = 'ADU'
+
+    header['FILTER'] = obs_info['band']
+    header['DETSN'] = '00000000000'
+    header['DETNAME'] = camera_config['detector_name']
+    header['CHIPLAB'] = camera_config['ccd_label']
+    header['CHIPTEMP'] = float(camera_config['chip_temp'])
+    header['DEWTEMP'] = float(camera_config['dewar_temp'])
+    header['DETSIZE'] = f"{imgszx} * {imgszy}"
+    header['IMGINDEX'] = index
+
+    frame_time = obs_info['expt'] + camera_config['readout_time']
+    bunch_start = datetime.fromisoformat(bunch_start)
+    expstart = bunch_start + timedelta(seconds=frame_time * (index - 1))
+    bunch_start_mjd = datetime_obj_to_mjd(bunch_start)
+
+    ra0 = primary_header.get('RA_PNT0', -1.0)
+    dec0 = primary_header.get('DEC_PNT0', -1.0)
+    pa0 = primary_header.get('POS_ANG0', -1.0)
+    cab0 = primary_header.get('CABSTART', bunch_start_mjd)
+
+    delta_t = frame_time * index
+    bunch_end = bunch_start + timedelta(seconds=delta_t)
+    bunch_end_mjd = datetime_obj_to_mjd(bunch_end)
+
+    ra1 = primary_header.get('RA_PNT1', ra0)
+    dec1 = primary_header.get('DEC_PNT1', dec0)
+    pa1 = primary_header.get('POS_ANG1', pa0)
+    cab1 = primary_header.get('CABEND', bunch_end_mjd)
+
+    img_cab = datetime_obj_to_mjd(expstart)
+
+    ratio = (img_cab - cab0)/(cab1 - cab0)
+    ra = ra0 + (ra1 - ra0) * ratio
+    dec = dec0 + (dec1 - dec0) * ratio
+    pa = pa0 + (pa1 - pa0) * ratio
+
+    header['IMG_EXPT'] = expstart.isoformat(timespec='seconds')
+    header['IMG_CABT'] = header['IMG_EXPT']
+    header['IMG_DUR'] = float(obs_info['expt'])
+
+    header['IMG_PA'] = ra
+    header['IMG_RA'] = dec
+    header['IMG_DEC'] = pa
+
+    header['DATASECT'] = f"{plszx} * {plszy}"
+    header['PIXSCAL'] = optics_config['platescale']
+    header['PIXSIZE'] = camera_config['pitch_size']
+    header['NCHAN'] = 1
+    header['PSCAN1'] = pscan1
+    header['PSCAN2'] = pscan2
+    header['OSCAN1'] = oscan1
+    header['OSCAN2'] = oscan2
+    header['UDARK'] = udark
+    header['BDARK'] = bdark
+    header['LDARK'] = ldark
+    header['RDARK'] = rdark
+
+    # WCS
+    cstar = {'ra': '0d', 'dec': '0d'}
+    if obs_info['target'] != {}:
+        cstar = obs_info['target']['cstar']
+
+    radec = SkyCoord(cstar['ra'], cstar['dec'])
+    shift = obs_info['shift']
+    platescale = optics_config['platescale']
+    rotation = np.radians(obs_info['rotation'])
+
+    header['WCSAXES'] = 2
+    header['CRPIX1'] = (plszx + 1)/2 + pscan1 + ldark + shift[0] / platescale
+    header['CRPIX2'] = (plszy + 1)/2 + pscan2 + udark + shift[0] / platescale
+    header['CRVAL1'] = radec.ra.degree
+    header['CRVAL2'] = radec.dec.degree
+    header['CTYPE1'] = 'RA---TAN'
+    header['CTYPE2'] = 'DEC--TAN'
+    header['CD1_1'] = np.cos(rotation)
+    header['CD1_2'] = -np.sin(rotation)
+    header['CD2_1'] = np.sin(rotation)
+    header['CD2_2'] = np.cos(rotation)
+    header['others'] = 'other'
+
+    # Readout information
+    header['EMGAIN'] = float(obs_info['emgain'])
+    header['GAIN'] = float(camera_config['ph_per_adu'])
+    header['DET_BIAS'] = float(camera_config['bias_level'])
+    header['RON'] = float(camera_config['readout_noise'])
+    header['READTIME'] = float(camera_config['readout_time'])
+    header['ROSPEED'] = float(camera_config['readout_speed'])
+
+    # CPIC information
+    header['LS_STAT'] = 'OFF'
+    header['IWA'] = optics_config['mask_width'] / 2
+    header['WFSINFO1'] = -1.0
+    header['WFSINFO2'] = -1.0
+
+    header['CHECKSUM'] = '0000000000000000'
+    header['DATASUM'] = '0000000000'
+
+    header = check_and_update_fits_header(header)
+
+    return header
+
+
+def save_fits_simple(images, obs_info):
+    """
+    Save the image to a fits file with a simple header to TMP directory.
+
+    Parameters
+    ----------
+    images : numpy.ndarray
+        Image array to be written.
+    obs_info : dict
+        Dictionary of observation informations. See `save_fits` function.
+
+    Returns
+    ----------
+    Filename of the saved fits file.
+
+    """
+    target = obs_info['target']
+    target_info = 'NO_TARGET'
+    if 'cstar' in target.keys():
+        target_info = ''
+        target_info = f"S{target['cstar']['magnitude']:.1f}"
+        target_info += f"_P{len(target.get('planets', '[]'))}"
+
+    name = target_info
+    if 'name' in target.keys():
+        name = target['name']
+
+    name = name.replace('/', '_')
+    name = name.replace(',', '_')
+
+    now = datetime.now()
+    time = now.strftime("%Y%m%d%H%M%S")
+    filename = f"{name}_{time}.fits"
+
+    header = fits.Header()
+    header['skybg'] = obs_info['skybg']
+    header['name'] = name
+    header['exptime'] = obs_info['expt']
+    header['nframe'] = obs_info['nframe']
+    header['band'] = obs_info['band']
+    header['emgain'] = obs_info['emgain']
+    header['obsid'] = obs_info['obsid']
+    header['rotation'] = obs_info['rotation']
+    shift = obs_info['shift']
+    header['shift'] = f"x:{shift[0]},y:{shift[1]}"
+
+    fullname = f"{tmp_folder_path}/{filename}"
+    fits.writeto(fullname, images, overwrite=True, header=header)
+    return fullname
+
+
+def save_fits(images, obs_info, gnc_info, csst_format=True):
+    """
+    Save the image to a fits file.
+
+    Parameters
+    ----------
+    images : numpy.ndarray
+        Image array to be saved.
+    obs_info : dict
+        Dictionary of observation informations.
+        Must contain the following keys
+
+        - band: str
+            - Band of the image.
+        - expt: float
+            - Exposure time of the each image.
+        - nframe: int
+            - Number of frames in the image.
+        - emgain: int
+            - EM gain of the camera.
+        - obsid: str
+            - Observation ID. Obsid must be 11 digits and start with 5 for CPIC. See pharse_obsid() for details.
+        - rotation: float
+            - Rotation angle of the image.
+        - shift: list
+            - Shift of the image.
+
+    gnc_info : dict
+        Dictionary of GNC information.
+        Contains the keywords in the primary header. See primary_hdu() for details.
+
+    csst_format : bool, optional
+        Whether to save the fits file in CSST format, by default True.
+    """
+
+    if not csst_format:
+        save_fits_simple(images, obs_info)
+        return
+
+    hdu_header, folder, filename = primary_hdu(obs_info, gnc_info, True)
+    hdu_list = fits.HDUList([hdu_header])
+
+    if len(images.shape) == 2:
+        images = np.array([images])
+
+    for index in range(images.shape[0]):
+        header = frame_header(
+            obs_info,
+            index + 1,
+            hdu_header.header['DATE-OBS'],
+            primary_header=hdu_header.header
+        )
+        frame_hdu = fits.ImageHDU(images[index, :, :], header=header)
+        frame_hdu.add_checksum()
+        hdu_list.append(frame_hdu)
+
+    folder = f"{output_dir}/{folder}"
+    if not os.path.exists(folder):
+        os.makedirs(folder)
+
+    hdu_list.writeto(f"{folder}/{filename}", overwrite=True)

CpicImgSim/main.py

+import numpy as np
+import re
+
+from .target import spectrum_generator
+from .optics import make_focus_image, focal_mask, optics_config
+from .psf_simulation import simulate_psf
+from .camera import EMCCD, CosmicRayFrameMaker, sky_frame_maker
+from .io import save_fits, log
+from .config import which_focalplane
+
+
+def psf_function(band, cstar_spectrum, shape, error=0.1):
+    cstar = True
+    if shape < 300:
+        cstar = False
+    return simulate_psf(error, band, cstar_spectrum, nsample=1, cstar=cstar)
+
+
+def observation_simulation(
+        target: dict,
+        skybg: float,
+        expt: float,
+        nframe: int,
+        band: str,
+        emgain: float,
+        obsid: int = 51900000000,
+        rotation: float = 0,
+        shift: list = [0, 0],
+        gnc_info: dict = {},
+        csst_format: bool = True,
+        psf_function: callable = psf_function):
+    """
+    Simulate the observation. All-In-One function of the package.
+
+    Parameters
+    -----------
+    target: dict
+        The target information. See target.py for details.
+    skybg: float
+        magnitude of the skybackground at the input b and. (abmag system)
+    expt: float
+        exposure time in second.
+    nframe: int
+        number of frames to be simulated.
+    band: str
+        the band of the observation. (e.g. 'f661')
+    emgain: float
+        the EM gain of the camera.
+    obsid: int
+        the observation id. Default is 51900000000.
+    rotation: float
+        the rotation angle of the target in degree. Default is 0.
+    shift: list
+        the shift of the target in arcsec. Default is [0, 0].
+    gnc_info: dict
+        the gnc information. Default is {}. See io.py for details.
+    csst_format: bool
+        whether to save the fits file in CSST format. Default is True.
+    psf_function: callable
+        the function to generate the psf. See optics.py for details.
+
+    Returns
+    -----------
+    np.ndarray of the simulated images with shape (nframe, 1088, 1050).
+
+    """
+    target_list = []
+    if 'cstar' in target.keys():
+        target_list = spectrum_generator(target)
+
+    focal_name = which_focalplane(band)
+    this_focal_config = optics_config[focal_name]
+    telescope_config = optics_config['telescope']
+    area = telescope_config['aperature_area']
+
+    if focal_name == 'vis':
+        camera = EMCCD()
+    else:
+        raise ValueError('Only VIS focal plane is supported.')
+
+    platescale = this_focal_config['platescale']
+    iwa = this_focal_config['mask_width'] / 2
+
+    crmaker = CosmicRayFrameMaker()
+    images = []
+
+    params = {
+        'target': target,
+        'skybg': skybg,
+        'expt': expt,
+        'nframe': nframe,
+        'band': band,
+        'emgain': emgain,
+        'obsid': obsid,
+        'rotation': rotation,
+        'shift': shift,
+    }
+    paramstr = ', '.join([f'{k}={v}' for k, v in params.items()])
+    log.debug(f"parameters: {paramstr}")
+
+    for i in range(nframe):
+
+        log.info(f'Simulation Running: Frame {i+1}/{nframe}')
+
+        focal_frame = make_focus_image(
+            band,
+            target_list,
+            psf_function,
+            platesize=camera.flat_shape,
+            rotation=rotation,
+            init_shifts=shift,
+        )
+
+        if skybg is None or skybg > 100:
+            sky_bkg_frame = 0
+        else:
+            sky_bkg_frame = sky_frame_maker(
+                band,
+                skybg,
+                platescale,
+                camera.flat_shape
+            )
+
+        focal_frame = (focal_frame + sky_bkg_frame) * area
+        focal_frame = focal_mask(focal_frame, iwa, platescale)
+
+        cr_frame = crmaker.make_cr_frame(camera.dark_shape, expt)
+
+        img = camera.readout(
+            focal_frame,
+            emgain,
+            expt,
+            image_cosmic_ray=cr_frame
+        )
+
+        images.append(img)
+
+    images = np.array(images)
+
+    save_fits(images, params, gnc_info, csst_format=csst_format)
+
+    return images
+
+
+def quick_run(
+        target_str: str,
+        skymag: float,
+        band: str,
+        expt: float,
+        nframe: int,
+        emgain: float,
+        rotation: float = 0,
+        shift: list = [0, 0]) -> np.ndarray:
+    """
+    A quick run function to simulate the observation.
+
+    Parameters
+    -----------
+    target_str: str
+        The target information in string format.
+        In the format of "\*5.1/25.3(1.3,1.5)/22.1(2.3,-4.5)" which means a central star
+        with magnitude 5.1, and two substellar with magnitude 25.3 and 22.1, respectively.
+        The first number in the parenthesis is the x position in arcsec, and the second is the y position.
+    skybg: float
+        magnitude of the skybackground at the input band. (abmag system)
+    band: str
+        the band of the observation. (e.g. 'f661')
+    expt: float
+        exposure time in second.
+    nframe: int
+        number of frames to be simulated.
+    emgain: float
+        the EM gain of the camera.
+    rotation: float (optional)
+        the rotation angle of the target in degree. Default is 0.
+    shift: list (optional)
+        the shift of the target in arcsec. Default is [0, 0].
+
+    Returns
+    -----------
+    np.ndarray of the simulated images, with shape (nframe, 1088, 1050)
+
+    Notes
+    -----------
+    1. stars are assumed to be G0III with distance 10pc.
+    2. magnitude of the star and substellar are assumed to be in the same band.
+    3. Csst format is not supported.
+    4. The psf is assumed to be the default one.
+    5. fits file will be saved in the current directory.
+
+
+    """
+    log.info(f'Quick Run: {target_str}')
+    target_dict = {
+        'name': 'cal',
+    }
+
+    if (target_str != '') and (target_str[0] == '*'):
+        objects = target_str[1:].split('/')
+        cstar_mag = float(objects[0])
+
+        cstar = {
+            'magnitude': cstar_mag,
+            'ra': '0d',
+            'dec': '0d',
+            'sptype': 'G0III',
+            'distance': 10,
+            'mag_input_band': band
+        }
+
+        stars = []
+        for sub_stellar in objects[1:]:
+
+            float_regex = R"[+-]?\d+(?:\.\d+)?"
+            match = re.match(
+                rf"({float_regex})\(({float_regex}),({float_regex})\)", sub_stellar)
+            if not match:
+                raise ValueError('Wrong format for sub stellar.')
+            mag = float(match.group(1))
+            x = float(match.group(2))
+            y = float(match.group(3))
+            pangle = np.arctan2(x, y) * 180 / np.pi
+            separation = np.sqrt(x**2 + y**2)
+            stars.append({
+                'magnitude': mag,
+                'pangle': pangle,
+                'separation': separation,
+                'sptype': 'G0III',
+                'mag_input_band': band
+            })
+        target_dict = {
+            'name': target_str[1:],
+            'cstar': cstar,
+            'stars': stars,
+        }
+
+    return observation_simulation(
+        target=target_dict,
+        skybg=skymag,
+        expt=expt,
+        nframe=nframe,
+        band=band,
+        emgain=emgain,
+        csst_format=False,
+        shift=shift,
+        rotation=rotation,
+    )
+
+
+# observation_simulation(
+#     target={},
+#     skybg=15,
+#     expt=10,
+#     nframe=2,
+#     band='f661',
+#     emgain=30,
+#     obsid=50112345678,
+# )
+
+# quick_run('*5.1/25.3(0.8,0.8)', None, 'f661', 10, 1, 10)
+# quick_run('*5/20(0.8,0.8)', None, 'f883', 10, 1, 10)
+
+# # quick *5.1/25.3(1.3,1.5) expt nframe emgain band rotation shift
+# # quick target_name expt nframe emgain band rotation shift
+# # plan plan_file_or_folder
+if __name__ == '__main__':  # pragma: no cover
+    target_example = {
+        'cstar': {
+            'magnitude': 1,
+            'ra': '120d',
+            'dec': '40d',
+            'distance': 10,
+            'sptype': 'F0III',
+        },
+        'stars': [
+            {
+                'magnitude': 20,
+                'pangle': 60,
+                'separation': 1,
+                'sptype': 'F0III'
+            }
+        ]
+    }
+#     quick_run('', 10, 'f661', 1, 1, 30)
+#     quick_run('*2.4/10(3,5)/15(-4,2)', 21, 'f661', 1, 1, 30)
+
+#     # normal target
+    observation_simulation(
+        target=target_example,
+        skybg=21,
+        expt=1,
+        nframe=2,
+        band='f661',
+        emgain=30,
+        obsid=51012345678,
+    )
+
+#     # bias
+    # observation_simulation(
+    #     target=target_example,
+    #     skybg=999,
+    #     expt=1,
+    #     nframe=2,
+    #     band='f661',
+    #     emgain=1,
+    #     obsid=51012345678,
+    #     shift=[3, 3],
+    #     rotation=60
+    # )
+
+#     # bias-gain
+#     observation_simulation(
+#         target={},
+#         skybg=999,
+#         expt=0.01,
+#         nframe=2,
+#         band='f661',
+#         emgain=1000,
+#         obsid=50012345678,
+#     )
+
+#     # dark
+#     observation_simulation(
+#         target={},
+#         skybg=999,
+#         expt=100,
+#         nframe=2,
+#         band='f661',
+#         emgain=30,
+#         obsid=50112345678,
+#     )
+
+#     # flat
+#     observation_simulation(
+#         target={},
+#         skybg=15,
+#         expt=10,
+#         nframe=2,
+#         band='f661',
+#         emgain=30,
+#         obsid=50112345678,
+#     )

CpicImgSim/optics.py

+import os
+import yaml
+import time
+
+import scipy as sp
+
+import numpy as np
+from CpicImgSim.config import cpism_refdata, which_focalplane, S  # S is synphot
+from CpicImgSim.config import optics_config
+from CpicImgSim.utils import region_replace
+from CpicImgSim.io import log
+from astropy.convolution import convolve_fft
+from scipy.signal import fftconvolve
+
+FILTERS = {
+    'f565': S.FileBandpass(f'{cpism_refdata}/throughtput/f565_total.fits'),
+    'f661': S.FileBandpass(f'{cpism_refdata}/throughtput/f661_total.fits'),
+    'f743': S.FileBandpass(f'{cpism_refdata}/throughtput/f743_total.fits'),
+    'f883': S.FileBandpass(f'{cpism_refdata}/throughtput/f883_total.fits'),
+    'f940': S.FileBandpass(f'{cpism_refdata}/throughtput/f940_total.fits'),
+    'f1265': S.FileBandpass(f'{cpism_refdata}/throughtput/f1265_total.fits'),
+    'f1425': S.FileBandpass(f'{cpism_refdata}/throughtput/f1425_total.fits'),
+    'f1542': S.FileBandpass(f'{cpism_refdata}/throughtput/f1542_total.fits'),
+}
+
+def filter_throughput(filter_name):
+    """
+    Totally throughput of the each CPIC band.
+    Including the throughput of the filter, telescope, cpic, and camera QE.
+    If the filter_name is not supported, return the throughput of the default filter(f661).
+
+    Parameters
+    -----------
+    filter_name: str
+        The name of the filter.
+        One of ['f565', 'f661'(default), 'f743', 'f883', 'f940', 'f1265', 'f1425', 'f1542']
+
+    Returns
+    --------
+    synphot.Bandpass
+        The throughput of the filter.
+
+    """
+    filter_name = filter_name.lower()
+    filter_name = 'f661' if filter_name == 'default' else filter_name
+    if filter_name not in FILTERS.keys():
+        log.warning(f"滤光片名称错误({filter_name})，返回默认滤光片(f661)透过率")
+        filter_name = 'f661'
+
+    return FILTERS[filter_name]
+
+def example_psf_func(band, spectrum, frame_size, error=0.1):
+    """
+    Example psf generating function.
+
+    Parameters
+    -------------
+    band: str
+        The name of the band.
+    spectrum: synphot.Spectrum or synphot.SourceSpectrum
+        The spectrum of the target.
+    frame_size: int
+        The size of the frame.
+    error: float
+        Phase RMS error.
+
+    Returns
+    ---------------
+    2D array
+        psf image with shape of `frame_size`
+
+    """
+    pass
+
+
+def example_monochromatic_psf(wavelength, error=0.1):
+    pass
+
+
+def rotate_and_shift(shift, rotation, init_shifts):
+    rotation_rad = rotation / 180 * np.pi
+    return np.array([
+        shift[0] * np.cos(rotation_rad) + shift[1] * np.sin(rotation_rad),
+        -shift[0] * np.sin(rotation_rad) + shift[1] * np.cos(rotation_rad)
+    ]) + np.array(init_shifts)
+
+
+from scipy.ndimage import rotate
+
+
+def ideal_focus_image(
+        bandpass: S.spectrum.SpectralElement,
+        targets: list,
+        platescale,
+        platesize: list = [1024, 1024],
+        init_shifts: list = [0, 0],
+        rotation: float = 0,):
+    
+    focal_image = np.zeros(platesize)
+    focal_shape = np.array(platesize)[::-1] # x, y
+
+    if not targets:
+        return focal_image
+    
+    for target in targets:
+        sub_x, sub_y, sub_spectrum, sub_image = target
+        sub_shift = rotate_and_shift([sub_x, sub_y], rotation, init_shifts) / platescale
+        sed = (sub_spectrum * bandpass).integrate()
+
+        if sub_image is None:
+
+            x = (focal_shape[0] - 1)/2 + sub_shift[0]
+            y = (focal_shape[1] - 1)/2 + sub_shift[1]
+
+            int_x = int(x)
+            int_y = int(y)
+            if int_x < 0 or int_x >= focal_shape[0] - 1 or int_y < 0 or int_y >= focal_shape[1] - 1:
+                continue
+
+            dx1 = x - int_x
+            dx0 = 1 - dx1
+            dy1 = y - int_y
+            dy0 = 1 - dy1
+
+            sub = np.array([
+                [dx0*dy0, dx1*dy0],
+                [dx0*dy1, dx1*dy1]]) * sed
+            
+            focal_image[int_y: int_y+2, int_x: int_x+2] += sub
+        else:
+            # sub_image = sub_image
+            sub_image = np.abs(rotate(sub_image, rotation, reshape=False))
+            sub_image = sub_image / sub_image.sum()
+            sub_img_shape = np.array(sub_image.shape)[::-1]
+            sub_shift += (focal_shape-1)/2 - (sub_img_shape-1)/2
+            focal_image = region_replace(
+                focal_image,
+                sub_image * sed,
+                sub_shift,
+                subpix=True
+            )
+    return focal_image
+
+from scipy.signal import fftconvolve
+def sp_convole_fft(image, kernal):
+    kernal = kernal / kernal.sum()
+    # y0 = kernal.shape[0] // 2
+    # x0 = kernal.shape[1] // 2
+    outimg = fftconvolve(image, kernal, mode='same')
+    # return outimg[y0:y0+image.shape[0], x0:x0+image.shape[1]]
+    return outimg
+def convolve_psf(
+        band: str,
+        targets: list,
+        psf_function: callable,
+        init_shifts: list = [0, 0],
+        rotation: float = 0,
+        nsample: int = 5,
+        error: float = 1,
+        platesize: list = [1024, 1024]) -> np.ndarray :
+
+    config = optics_config[which_focalplane(band)]
+    platescale = config['platescale']
+
+    filter = filter_throughput(band)
+    wave = filter.wave
+    throughput = filter.throughput
+    min_wave = wave[0]
+    max_wave = wave[-1]
+    
+    all_fp_image = []
+    for i_wave in range(nsample):
+            d_wave = (max_wave - min_wave) / nsample
+            wave0 = min_wave + i_wave * d_wave
+            wave1 = min_wave + (i_wave + 1) * d_wave
+            center_wavelength = (wave0 + wave1) / 2 * 1e-10
+
+            i_throughput = throughput.copy()
+            i_throughput[(wave > wave1) | (wave < wave0)] = 0
+            i_band = S.ArrayBandpass(wave, i_throughput, waveunits='angstrom')
+
+            i_fp_image = ideal_focus_image(i_band, targets, platescale, platesize, init_shifts, rotation)
+            psf = psf_function(center_wavelength, error=error)
+            t0 = time.time()
+
+
+            # c_fp_image = convolve_fft(i_fp_image, psf, allow_huge=True)
+            c_fp_image = sp_convole_fft(i_fp_image, psf)
+
+            print(f"Convolution time: {time.time()-t0}")
+            all_fp_image.append(c_fp_image)
+
+    return np.array(all_fp_image).mean(axis=0)
+    
+
+def make_focus_image(
+        band: str,
+        targets: list,
+        psf_function: callable,
+        init_shifts: list = [0, 0],
+        rotation: float = 0,
+        platesize: list = [1024, 1024]) -> np.ndarray:
+    """
+    Make the focus image of the targets.
+
+    Parameters
+    -----------
+    band: str
+        The name of the band.
+    targets: list
+        The list of the targets.
+        Each element of the list is a tuple of (x, y, spectrum).
+
+        - x, y: float
+            - The position of the target in the focal plane.
+        - spectrum: synphot.Spectrum or synphot.SourceSpectrum
+            - The spectrum of the target.
+    psf_function: callable
+        The function to generate the PSF, with same parameters and return as `example_psf_func`.
+
+    init_shifts: list
+        The initial shifts of the center targets. Unit: arcsec.
+        The default is [0, 0].
+    rotation: float
+        The rotation of the focal plane. Unit: degree.
+        The default is 0 degree.
+    platesize: list
+        The size of the focal plane. Unit: pixel.
+        The default is [1024, 1024].
+
+    Returns
+    --------
+    np.ndarray
+        The focus image of the targets.
+        2D array with the shape of platesize.
+    """
+
+    config = optics_config[which_focalplane(band)]
+    platescale = config['platescale']
+
+    focal_image = np.zeros(platesize)
+    if not targets:
+        return focal_image
+
+    cstar_x, cstar_y, cstar_spectrum = targets[0]
+    cstar_shift = rotate_and_shift([cstar_x, cstar_y]) / platescale
+
+    error_value = 0  # nm
+
+    cstar_psf = psf_function(band, cstar_spectrum, config['cstar_frame_size'],
+                             error=error_value)
+
+    platesize = np.array(platesize)[::-1]
+    psf_shape = np.array(cstar_psf.shape)[::-1]
+    cstar_shift += (platesize-1)/2 - (psf_shape-1)/2
+
+    focal_image = region_replace(
+        focal_image,
+        cstar_psf,
+        cstar_shift,
+        padded_in=False,
+        padded_out=False,
+        subpix=True)
+
+    for i_target in range(1, len(targets)):
+        sub_x, sub_y, sub_spectrum = targets[i_target]
+        pdout = False if i_target == len(targets)-1 else True
+        pdin = False if i_target == 1 else True
+        log.debug(f"input target {sub_x=:}, {sub_y=:}")
+        sub_shift = rotate_and_shift([sub_x, sub_y], rotation, init_shifts) / platescale
+        log.debug(f"after rotate and shift {sub_shift=:}")
+        sub_psf = psf_function(
+            band,
+            sub_spectrum,
+            config['substellar_frame_size'],
+            error=error_value
+        )
+        psf_shape = np.array(sub_psf.shape)[::-1]
+        sub_shift += (platesize-1)/2 - (psf_shape-1)/2
+        log.debug(f"input shift of region_replace: {sub_shift=:}")
+        focal_image = region_replace(
+            focal_image,
+            sub_psf,
+            sub_shift,
+            padded_in=pdin,
+            padded_out=pdout,
+            subpix=True
+        )
+
+    return focal_image
+
+def focal_mask(image, iwa, platescale, throughtput=1e-6):
+    """
+    Mask the image outside the inner working angle.
+
+    Parameters
+    -----------
+    image: np.ndarray
+        The image to be masked.
+    iwa: float
+        The inner working angle. Unit: arcsec.
+    platescale: float
+        The platescale of the image. Unit: arcsec/pixel.
+    throughtput: float
+        The throughtput of the mask. The default is 1e-6.
+
+    Returns
+    --------
+    np.ndarray
+        The masked image.
+    """
+    xx, yy = np.mgrid[0:image.shape[0], 0:image.shape[1]]
+    center = np.array([(image.shape[0]-1)/2, (image.shape[1]-1)/2])
+    mask = (abs(xx - center[0]) < iwa /
+            platescale) | (abs(yy - center[1]) < iwa / platescale)
+    image_out = image.copy()
+    image_out[mask] *= throughtput
+    return image_out

CpicImgSim/psf_simulation.py

+import numpy as np
+from astropy.io import fits
+
+from hcipy import Field, Wavefront, DeformableMirror, FraunhoferPropagator
+from hcipy import SurfaceApodizer, SurfaceAberrationAtDistance
+from hcipy import make_pupil_grid, make_circular_aperture, make_focal_grid
+from hcipy import make_xinetics_influence_functions
+from hcipy import read_fits
+
+from .config import cpism_refdata, S
+from .optics import filter_throughput
+
+
+# initial psf simulation
+apm, apm_header = fits.getdata(
+    f'{cpism_refdata}/optics/apm.fits', header=True)
+actuator = read_fits(f'{cpism_refdata}/optics/actuator.fits')
+surface_aberration = read_fits(
+    f'{cpism_refdata}/optics/initial_phase_aberration.fits')
+
+wavelength = 625e-9  # m
+pupil_diameter = 2  # m
+focal_length = pupil_diameter * 83
+pupil_grid = make_pupil_grid(apm.shape[0], apm.shape[0] * apm_header['PHRATE'])
+aperture = make_circular_aperture(pupil_diameter)(pupil_grid)
+aperture = aperture * Field(apm.flatten(), pupil_grid)
+
+
+emccd_pixel_size = 13e-6  # m
+arcsec2rad = np.radians(1 / 3600)
+emccd_pixel_scale = emccd_pixel_size / \
+    (arcsec2rad * focal_length)  # arcsec / pixel
+spatial_resolution = focal_length * wavelength / \
+    pupil_diameter  # meter per airy disk
+
+q = spatial_resolution / emccd_pixel_size  # pixels per airy disk
+num_airy = 512 / q  # airy disk per frame (2 * 512 = 1024 pix)
+focal_full_frame = make_focal_grid(
+    q, num_airy, spatial_resolution=spatial_resolution)
+prop_full_frame = FraunhoferPropagator(
+    pupil_grid, focal_full_frame, focal_length)
+
+num_airy = 128 / q  # make a small frame for the planets
+focal_sub_frame = make_focal_grid(
+    q, num_airy, spatial_resolution=spatial_resolution)
+prop_sub_frame = FraunhoferPropagator(
+    pupil_grid, focal_sub_frame, focal_length)
+
+num_actuators_across = 32
+# dm spacing is little smaller than pupil
+actuator_spacing = 0.95 / 32 * pupil_diameter
+influence_functions = make_xinetics_influence_functions(
+    pupil_grid, num_actuators_across, actuator_spacing)
+deformable_mirror = DeformableMirror(influence_functions)
+
+aberration = SurfaceApodizer(
+    surface_sag=surface_aberration.flatten(), refractive_index=-1)
+# arbitrary distance for the aberration to propagate
+aberration_distance = 80 * focal_length
+aberration = SurfaceAberrationAtDistance(aberration, aberration_distance)
+
+
+def single_band_psf(wavelength, waveerror=0, aber_phase=None):
+    error = np.random.normal(0, waveerror*1e-9, actuator.shape)
+    deformable_mirror.actuators = actuator + error
+    wf = Wavefront(aperture, wavelength)
+    wf = aberration(wf)
+    if aber_phase is not None:
+        other_error = SurfaceApodizer(
+            surface_sag=aber_phase.flatten(), refractive_index=-1)
+        wf = other_error(wf)
+    img = prop_full_frame(deformable_mirror(wf)).intensity.shaped 
+    return img
+
+def simulate_psf(error_value, band, spectrum, nsample=5, cstar=True, aber_phase=None):
+    filter = filter_throughput(band)
+    wave = filter.wave
+    throughput = filter.throughput
+    min_wave = wave[0]
+    max_wave = wave[-1]
+
+    sed = []
+    sed_center_wavelength = []
+    for i_wave in range(nsample):
+        d_wave = (max_wave - min_wave) / nsample
+        wave0 = min_wave + i_wave * d_wave
+        wave1 = min_wave + (i_wave + 1) * d_wave
+        sed_center_wavelength.append((wave0 + wave1) / 2 * 1e-10)
+
+        i_throughput = throughput.copy()
+        i_throughput[(wave > wave1) | (wave < wave0)] = 0
+        i_band = S.ArrayBandpass(wave, i_throughput, waveunits='angstrom')
+        i_sed = (spectrum * i_band).integrate()
+        sed.append(i_sed)
+
+    error = np.random.normal(0, error_value*1e-9, actuator.shape)
+
+
+    imgs = []
+    deformable_mirror.actuators = actuator + error
+
+    prop = prop_full_frame
+    if not cstar:
+        prop = prop_sub_frame
+
+    for i_sample in range(nsample):
+        wf = Wavefront(aperture, sed_center_wavelength[i_sample])
+        wf = aberration(wf)
+        if aber_phase is not None:
+            other_error = SurfaceApodizer(
+                surface_sag=aber_phase.flatten(), refractive_index=-1)
+            wf = other_error(wf)
+
+        img = prop(deformable_mirror(wf)).intensity.shaped
+        imgs.append(img / img.sum() * sed[i_sample])
+
+    return np.array(imgs).sum(axis=0)

CpicImgSim/target.py

+import os
+import re
+import json
+import numpy as np
+from scipy import constants
+from astropy.io import fits
+from astropy.coordinates import SkyCoord
+
+from CpicImgSim.config import cpism_refdata, MAG_SYSTEM
+from CpicImgSim.config import S  # pysynphot
+from CpicImgSim.optics import filter_throughput
+from CpicImgSim.io import log
+
+PLANK_H = constants.h  # ~6.62607015e-34
+LIGHT_SPEED = constants.c  # ~3e8
+DEG2RAD = np.pi / 180
+R_JUPITER_METER = 6.99e4
+AU_METER = 1.49e8
+
+DEFAULT_FILTER_FILE = f'{cpism_refdata}\\throughtput\\f661_total.fits'
+HYBRID_MODEL_FILE = f'{cpism_refdata}\\target_model\\hybrid_model.fits'
+BCC_MODEL_FOLDER = f'{cpism_refdata}\\target_model\\bccmodels'
+
+CATALOG_CACHE = {}
+
+class AlbedoCat(S.spectrum.TabularSpectralElement, S.catalog.Icat):
+    """
+    This class is used to generate a spectrum from the planet reflection model.
+    
+    Parameters
+    ----------
+    phase : float
+        [degree], The phase angle of the planet, range from 0 to 180.
+    metallicity : float
+        The metallicity of the planet. log(Z/Zsun), range from 0 to 2.
+    f_sed : float
+        The sedimentation efficiency of the cloud. log(f_sed), range from -2 to 2. 2 is for cloud free case.
+
+    Notes
+    -----
+    f_sed is only reliable between -2 and log10(6) and the cloud free case (2).
+    Values between log10(6) and 2 are interpolated from the cloud free case (2) and log10(6).
+    
+    """
+    
+    def __init__(self, 
+                phase: float, 
+                metallicity: float,
+                f_sed: float,
+            ):
+        catname = 'BCCalbedo'
+        self.isAnalytic=False
+        self.name = f"{catname}_{phase}_{metallicity}_{f_sed}"
+        self.parameter_names = ['phase', 'metallicity', 'f_sed']
+        self.catalog_folder = BCC_MODEL_FOLDER
+        filename = os.path.join(self.catalog_folder, 'catalog.fits')
+
+        if filename in CATALOG_CACHE:
+            indices = CATALOG_CACHE[filename]
+        else:
+            with fits.open(filename) as table:
+                indexList = table[1].data.field('INDEX')
+                filenameList = table[1].data.field('FILENAME')
+            
+            indices = self._getArgs(indexList, filenameList)
+            CATALOG_CACHE[filename] = indices
+
+        list0,list1 = self._breakList(indices, 0, phase)
+
+        list2,list3 = self._breakList(list0, 1, metallicity)
+        list4,list5 = self._breakList(list1, 1, metallicity)
+
+        list6,list7   = self._breakList(list2, 2, f_sed)
+        list8,list9   = self._breakList(list3, 2, f_sed)
+        list10,list11 = self._breakList(list4, 2, f_sed)
+        list12,list13 = self._breakList(list5, 2, f_sed)
+
+        sp1, wave, waveunits = self._getSpectrum(list6[0],  catname, wave_output=True)
+        sp2 = self._getSpectrum(list7[0],  catname)
+        sp3 = self._getSpectrum(list8[0],  catname)
+        sp4 = self._getSpectrum(list9[0],  catname)
+        sp5 = self._getSpectrum(list10[0], catname)
+        sp6 = self._getSpectrum(list11[0], catname)
+        sp7 = self._getSpectrum(list12[0], catname)
+        sp8 = self._getSpectrum(list13[0], catname)
+
+        spa1 = self._interpolateSpectrum(sp1, sp2, f_sed)
+        spa2 = self._interpolateSpectrum(sp3, sp4, f_sed)
+        spa3 = self._interpolateSpectrum(sp5, sp6, f_sed)
+        spa4 = self._interpolateSpectrum(sp7, sp8, f_sed)
+
+        spa5 = self._interpolateSpectrum(spa1, spa2, metallicity)
+        spa6 = self._interpolateSpectrum(spa3, spa4, metallicity)
+
+        spa7 = self._interpolateSpectrum(spa5, spa6, phase)
+
+        sp = spa7[0]
+
+        self._wavetable = wave * 1e4
+        self._throughputtable = sp
+        self.waveunits = S.units.Units(waveunits.lower())
+        self.warnings = {}
+
+    def _getSpectrum(self, parlist, catdir, wave_output=False):
+        name = parlist[3]
+
+        filename = name.split('[')[0]
+        column = name.split('[')[1][:-1]
+
+        filename = f"{self.catalog_folder}/{filename}"
+        
+        # sp = S.spectrum.TabularSpectralElement(filename, thrucol=column)
+
+        with fits.open(filename) as td:
+            sp = td[1].data.field(column)
+            wave = td[1].data.field('wavelength')
+            # waveunits = td[1].header['tunit1']
+            waveunits = 'micron'
+
+        result = []
+        for member in parlist:
+            result.append(member)
+
+        result.pop()
+        result.append(sp)
+
+        if wave_output:
+            return result, wave, waveunits
+
+        return result
+
+
+
+def _sptype2num(spectral_type):
+    """
+    convert spectral type string to number, for interpretation
+
+    case0: normal case
+    - M1V: 6.1, 5
+    - O5IV: 0.5, 4
+    - F3V: 3.3, 5
+    - K4.5II: 5.45, 2
+
+    case 1: star type or subtype missing
+    zero or V will return
+    - M1: 6.1, 5
+    - M: 6.0, 5
+
+    case 2: spectral type + subtype
+    subtype will be ignored
+    - K3Vvar: 5.3, 5
+    - F6Vbwvar: 3.6, 5
+    - K0IV SB: 5.0, 4
+    - F5V+: 3.5, 5
+
+    case 3: multi spectral type
+    only the first sptype is used
+    - G5IV/V +K1IV: 4.5, 4
+    - F7IV-V: 3.7, 4
+    - O4/O5IV: 0.4, 0
+    - G3/G5V: 4.3, 0
+
+    case 4: illegal input
+    ValueError will be raised
+    """
+    obafgkm = 'OBAFGKML'
+    spectral_type = spectral_type.upper()
+    # match spectral type such as M1V, O5IV, F3V, K4.5II
+    matched = re.match(
+        RF'^([{obafgkm}])([0-9]\d*\.?\d*)*([IV]*)', spectral_type)
+
+    if not matched:
+        raise ValueError(f"illegal spectral type input: {spectral_type}")
+
+    shorttype = obafgkm.find(matched.group(1))
+
+    subtype = 0.0
+    if matched.group(2):
+        subtype = float(matched.group(2))
+
+    stlist = ['O', 'I', 'II', 'III', 'IV', 'V']
+    startype = dict(zip(stlist, range(len(stlist)))).get(matched.group(3), 5)
+
+    return shorttype + subtype / 10, startype
+
+
+def _spnum2teff(spn, stn):
+    """
+    interpret of spectral number (by __sptype2num) to get t_eff and log_g
+    look up table from the document of ck04model
+    """
+    with open(f'{cpism_refdata}\\target_model\\sptype2teff_lut.json', 'r') as fid:
+        sptype_teff_lut = json.load(fid)
+
+    def _interp(spn, stn):
+        lut = sptype_teff_lut[f'{stn}']
+        teff = np.interp(spn, lut[0], lut[1])
+        logg = np.interp(spn, lut[0], lut[2])
+        return teff, logg
+
+    stn = 5 if stn not in [1, 2, 3, 4, 5] else stn
+
+    if stn in [1, 3, 5]:
+        return _interp(spn, stn)
+    else:
+        teff_low, logg_low = _interp(spn, stn-1)
+        teff_high, logg_high = _interp(spn, stn+1)
+        return (teff_high + teff_low)/2, (logg_high + logg_low)/2
+
+
+def star_photlam(
+        magnitude: float,
+        sptype: str,
+        is_blackbody: bool = False,
+        mag_input_band: str = 'f661'):
+    """
+    genrate flux spectrum of a star by its observal magnitude and spectral type
+
+    Parameters
+    ----------
+    magnitude: float
+        magnitude of the star
+    sptype: str
+        spectral type of the star
+    is_blackbody: bool
+        if True, use blackbody spectrum instead of ck04model
+    mag_input_band: str
+        bandpass of the input magnitude, default is f661
+
+    Returns
+    -------
+    pysynphot.spectrum.ArraySpectrum
+        spectrum of the star in photlam unit
+    """
+    spn, stn = _sptype2num(sptype)
+    t_eff, log_g = _spnum2teff(spn, stn)
+
+    log.debug(f"{sptype} star => [{t_eff=:}, {log_g=:}]; {is_blackbody=:}")
+    filter = filter_throughput(mag_input_band)
+
+    if not is_blackbody:
+        METALLICITY = 0
+        spectrum = S.Icat('ck04models', t_eff, METALLICITY, log_g)
+    else:
+        spectrum = S.BlackBody(t_eff)
+
+    star_sp = spectrum.renorm(magnitude, MAG_SYSTEM, filter)
+    star_sp.convert('photlam')
+    return star_sp
+
+
+def bcc_spectrum(
+        coe_cloud: float,
+        coe_metalicity: float
+):
+    """Albedo spectrum of planet using BCC model (Batalha et al. 2018),
+
+    Parameters
+    ----------
+    coe_cloud: float
+        The sedimentation efficiency of the cloud. log(f_sed), range from -2 to 2. 2 is for cloud free case.
+    coe_metalicity: float
+        The metallicity of the planet. log(Z/Zsun), range from 0 to 2.
+
+    Returns
+    -------
+    pysynphot.spectrum.ArrayBandpass
+        albedo spectrum of the planet
+    """
+
+    spec = AlbedoCat(0, coe_metalicity, coe_cloud)
+    spec.convert('angstrom')
+    return spec
+
+
+def hybrid_albedo_spectrum(
+        coe_b: float,
+        coe_r: float):
+    """Albedo spectrum of planet using hybrid-jupiter-neptune model (Lacy et al. 2018)
+    jupiter and neptune spectrum is from Karkoschka’s 1994
+
+    Parameters
+    ----------
+    coe_b: float
+        coefficient of blue spectrum, 1 for jupiter, 0 for neptune
+    coe_r: float
+        coefficient of red spectrum, 1 for jupiter, 0 for neptune
+
+    Returns
+    -------
+    pysynphot.spectrum.ArrayBandpass
+        albedo spectrum of the planet
+    """
+    log.debug(f"planet hybrid spectrum with {coe_b=:}, {coe_r=:}")
+    model = fits.getdata(HYBRID_MODEL_FILE)
+    spec = model[1, :] * coe_r
+    spec += model[2, :] * coe_b
+    spec += model[3, :] * (1 - coe_r)
+    spec += model[4, :] * (1 - coe_b)
+
+    albedo = S.ArrayBandpass(
+        wave=model[0, :],
+        throughput=spec,
+        waveunits='nm'
+    )
+    albedo.convert('angstrom')
+    return albedo
+
+
+def extract_target_x_y(
+        target: dict,
+        ra0: str = None,
+        dec0: str = None):
+    """
+    extract x, y of target from target dict
+
+    Parameters
+    ----------
+    target: dict
+        target dict. must contain either (ra, dec) or (pangle, spearation)
+    ra0: str
+        ra of center star. must be provided if (ra, dec) of target is used
+    dec0: str
+        dec of center star. must be provided if (ra, dec) of target is used
+
+    Returns
+    -------
+    x, y: float
+        x, y of target in arcsec
+
+    Raises
+    ------
+    ValueError
+        if (ra, dec) of target is used but (ra, dec) of center star is not provided.
+
+    ValueError
+        one of (ra, dec) or (pangle, spearation) is not provided.
+    """
+
+    def _pa2xy(p_angle, separation):
+        p_angle_rad = p_angle * DEG2RAD
+        x = separation * np.sin(p_angle_rad)
+        y = separation * np.cos(p_angle_rad)
+        log.debug(f"({p_angle=:}, {separation=:}) => ({x=:}, {y=:})")
+        return x, y
+
+    if 'pangle' in target.keys() and 'separation' in target.keys():
+        return _pa2xy(target['pangle'], target['separation'])
+
+    if 'ra' not in target.keys() or 'dec' not in target.keys():
+        raise ValueError(
+            'either (ra, dec) or (pangle, separation) needed in target dict')
+
+    if ra0 is None or dec0 is None:
+        raise ValueError(
+            '(ra, dec) of center star must be provided if (ra, dec) of bkstar is used'
+        )
+
+    ra, dec = target['ra'], target['dec']
+    log.debug(f"target: {ra=:}, {dec=:}, center star: {ra0=:}, {dec0=:}")
+    cstar = SkyCoord(ra0, dec0)
+    bkstar = SkyCoord(ra, dec)
+    separation = cstar.separation(bkstar).arcsec
+    p_angle = cstar.position_angle(bkstar).degree
+    x, y = _pa2xy(p_angle, separation)
+
+    return x, y
+    
+
+class TargetOjbect(object):
+    """A helper class to generate target spectrum and albedo spectrum
+
+    Attributes
+    ----------
+    x: float
+        x of target in arcsec
+    y: fload
+        y of target in arcsec
+    ra: str
+        ra string for center star, such as '15d23m05s'
+    dec: str
+        dec string for center star
+    distance: float
+        distance of center star in pc
+    image: 2d np.array
+        image of the target
+    spectrum: pysynphot.spectrum.Spectrum
+        spectrum of the target
+    """
+
+    def __init__(self, info, cstar=None, sp_model=None):
+        """Initialize a target object
+
+        Parameters
+        ----------
+        info: dict
+            target info, see Example for more details
+        cstar: TargetOjbect or None
+            center star object bounded, if None, means the target is the center star itself
+            center star object is used to calculate the x, y of each target according to its ra and dec
+            also center star's distance is used to calculate seperation of planet
+        sp_model: str
+            keyword to compatible with V1.0 input. See Notes
+            spectral type of the target, can be 'star' or 'planet'
+
+        Notes
+        -----
+        keyargv 中的sp_model 完全用来用来兼容V1.0版的输入,执行逻辑如下:
+         1. 如果没有设置sp_model，对应了V2.0的输入，则sp_model从info中获取
+         2. 如果设置了sp_model，则区分为star和planet
+         3. 如果sp_model为star, 则判断是否设置了isblackbody，如果是的话，则sp_model改为blackbody
+         4. 如果sp_model为planet，则进一步区分为hybrid_planet 和 bcc_planet 和 template_planet
+         注意：2.0版本的输入尽量不要使用sp_model关键字,逻辑还是有点复杂的，容易乱，而应该直接在info中设置
+
+        Examples: 
+        --------
+
+        cstar = {
+            'magnitude': 0,
+            'ra': '120d',
+            'dec': '40d',
+            'distance': 10,
+            'sptype': 'G0III'
+        }
+             
+        stars = {
+            'magnitude': 15,
+            'ra': '120.001d',
+            'dec': '40.001d',
+            'sptype': 'F0III',
+            'sp_spectrum': 'blackbody',
+        }
+            
+        planets = {
+            'radius': 2,
+            'pangle': 60,
+            'coe_cloud': 0.3,
+            'coe_metal': 0.7,
+            'separation': 0.5,
+            'phase_angle': 90,                 
+            'sp_spectrum': 'hybrid_planet',
+            'image': 'extend_planet.fits'
+        }
+
+        # planet using input costum albedo spectrum!
+        # Note that the albedo spectrum is not normalized! 
+        # so the contrast of the planet sould be considered in the input file by user!
+        # The input file is in pysynphot.FileSpectralElement format.
+        # See the documentation of pysynphot for details. 
+
+        planets = {
+            'pangle': 60,
+            'separation': 0.5,              
+            'sp_spectrum': 'template_planet',
+            'template': 'planet_albedo.fits'
+        }
+        """
+
+        if sp_model is None:  
+            sp_model = info.get('sp_type', 'star') # 如果不设置，默认为star
+        elif sp_model == 'star' and info.get('isblackbody', False):
+            #原先的star 现在分为star和blackbody两个类型
+            sp_model = 'blackbody'
+        elif sp_model == 'planet':
+            # 原先的planet 现在分为hybrid_planet 和 bcc_planet两类
+            sp_model = info.get('model', 'hybrid_planet')
+            if sp_model not in ['bcc_planet', 'hybrid_planet', 'template_planet']:
+                Warning(f'Unknown planet model {sp_model}, using default[hybrid] planet model' )
+
+        self.sp_model = sp_model
+
+        if cstar is None:
+            self.x, self.y = 0, 0
+            self.ra = info['ra']
+            self.dec = info['dec']
+            self.distance = info.get('distance', None)
+        else:
+            self.x, self.y = extract_target_x_y(info, cstar.ra, cstar.dec)
+
+
+        self.image = None
+        if 'image' in info.keys():
+            file = info['image']
+            if os.path.exists(file):
+                self.image = fits.getdata(file)    
+            else:
+                self.image = None
+                Warning(f"extend Ojbect File {file} not found, return None!")
+        
+
+        if self.sp_model == 'blackbody':
+            self.spectrum = star_photlam(
+                info['magnitude'],
+                info['sptype'],
+                is_blackbody=True,
+                mag_input_band=info.get('mag_input_band', 'f661')
+            )
+
+        if self.sp_model == 'template_star':
+            self.spectrum = S.FileSpectrum(info['template'])
+
+        if self.sp_model == 'star':
+            self.spectrum = star_photlam(
+                info['magnitude'],
+                info['sptype'],
+                is_blackbody=False,
+                mag_input_band=info.get('mag_input_band', 'f661')
+            )
+
+        if self.sp_model in ['hybrid_planet', 'bcc_planet']:
+            planet = info
+            phase_angle = planet.get('phase_angle', 90)
+
+            radius = planet.get('radius', 1)
+
+            if cstar.distance is None:
+                raise ValueError('distance of center star must be provided if planet is added')
+
+
+            if planet.get('contrast', None) is not None:
+                contrast = planet['contrast']
+            else:
+                contrast = planet_contrast(
+                    self.x * cstar.distance,
+                    self.y * cstar.distance,
+                    phase_angle,
+                    radius,
+                )
+
+            if sp_model == 'hybrid_planet':
+                coe_blue, coe_red = planet.get('coe_b', 1), planet.get('coe_r', 1)
+                albedo_spect = hybrid_albedo_spectrum(coe_blue, coe_red) * contrast
+            elif sp_model == 'bcc_planet':
+                coe_c, coe_m = planet.get('coe_cloud', 2), planet.get('coe_metal', 0)
+                albedo_spect = bcc_spectrum(coe_c, coe_m) * contrast
+            if sp_model == 'template_planet':
+                albedo_spect = S.FileBandpass(planet['template'])
+
+            self.spectrum = cstar.spectrum * albedo_spect 
+
+
+
+def planet_contrast(
+        planet_x_au: float,
+        planet_y_au: float,
+        phase_angle: float,
+        radius: float):
+    """
+    calculate the contrast of a planet
+
+    Parameters
+    ----------
+    planet_x_au: float
+        x position of the planet in au
+    planet_y_au: float
+        y position of the planet in au
+    phase_angle: float
+        phase angle of the planet in degree
+    radius: float
+        radius of the planet in jupiter radius
+
+    Returns
+    -------
+    contrast: float
+        contrast of the planet
+    """
+    separation = np.sqrt(planet_x_au**2 + planet_y_au**2)
+    phase_angle = phase_angle * DEG2RAD
+
+    if np.sin(phase_angle) < 1e-9:
+        raise ValueError('sin(phase_angle) can not be 0')
+
+    sep_3d = separation / np.sin(phase_angle)
+
+    # Lambert Scattering phase function
+    # from Madhusudhan and Burrows 2012 equation 33.
+    phase = (np.sin(phase_angle) + (np.pi - phase_angle)
+             * np.cos(phase_angle)) / np.pi
+    log.debug(f'alpha: {phase_angle/np.pi*180} {phase=}')
+
+    contrast = (radius / sep_3d * R_JUPITER_METER / AU_METER)**2 * phase
+    return contrast
+
+def spectrum_generator(
+        targets: dict):
+    """
+    generate the spectrum due to the input target list
+
+    Parameters
+    ----------
+    targets: dict
+        target dictionary which contains keyword 'cstar' (necessary), 'stars'(optional), 'planets' (optional).
+        The values are: 
+        
+        - cstar: dict
+            - center star information. must contain keywords ra, dec, distance, magnitude, sptype
+        - stars: list of dict, optional, not preferred, compatible with old version
+            - list of background stars. each dict must contain ra, dec, magnitude, sptype
+        - planets: list of dict, optional, not preferred, compatible with old version
+            - list of planets. each dict must contain pangle, separation, magnitude, radius
+        - objects: list of dict, optional, recommended! V2.0 new!
+            - list of targets. each dict must contain ra, dec, magnitude, sptype
+
+    Returns
+    -------
+    obj_sp_list: list
+        list of [x, y, spectrum, image] of each target
+    """
+
+    cstar = targets['cstar']
+    stars = targets.get('stars', [])
+    planets = targets.get('planets', [])
+    objects = targets.get('objects', [])
+
+    obj_sp_list = []
+
+    cstar_obj = TargetOjbect(cstar)
+    obj_sp_list.append([cstar_obj.x, cstar_obj.y, cstar_obj.spectrum, cstar_obj.image])
+
+    for star in stars:
+        star_obj = TargetOjbect(star, sp_model='star', cstar=cstar_obj)
+        obj_sp_list.append([star_obj.x, star_obj.y, star_obj.spectrum, star_obj.image])
+
+    for planet in planets:
+        planet_obj = TargetOjbect(planet, sp_model='planet', cstar=cstar_obj)
+        obj_sp_list.append([planet_obj.x, planet_obj.y, planet_obj.spectrum, planet_obj.image])
+
+    for target in objects:
+        target_obj = TargetOjbect(target, cstar=cstar_obj)
+        obj_sp_list.append([target_obj.x, target_obj.y, target_obj.spectrum, target_obj.image])
+
+    return obj_sp_list
\ No newline at end of file

CpicImgSim/utils.py

+import numpy as np
+import scipy.ndimage as nd
+import logging
+import random
+# DO NOT IMPORT CPICIMGSIM MODULES HERE
+
+
+class Logger(object):
+    def __init__(self, filename, level='INFO'):
+        self.logger = logging.getLogger('cpism_log')
+        self.logger.setLevel(logging.DEBUG)
+
+        shinfo = logging.StreamHandler()
+        onlyinfo = logging.Filter()
+        onlyinfo.filter = lambda record: (record.levelno < logging.WARNING)
+        fmtstr = '%(message)s'
+        shinfo.setFormatter(logging.Formatter(fmtstr))  # 设置屏幕上显示的格式
+        shinfo.setLevel(logging.INFO)
+        shinfo.addFilter(onlyinfo)
+
+        sh = logging.StreamHandler()
+        fmtstr = '!%(levelname)s!: %(message)s [%(filename)s - %(funcName)s (line: %(lineno)d)]: '
+        sh.setFormatter(logging.Formatter(fmtstr))  # 设置屏幕上显示的格式
+        sh.setLevel(logging.WARNING)
+
+        th = logging.FileHandler(filename)  # 往文件里写入#指定间隔时间自动生成文件的处理器
+
+        fmtstr = '%(asctime)s %(filename)s [%(funcName)s] - %(levelname)s: %(message)s'
+        th.setFormatter(logging.Formatter(fmtstr))  # 设置文件里写入的格式
+        th.setLevel(logging.__dict__.get(level.upper()))
+        self.logger.addHandler(shinfo)
+        self.logger.addHandler(sh)
+        self.logger.addHandler(th)
+
+
+def random_seed_select(seed=-1):
+    """
+    Select a random seed for numpy.random and return it.
+    """
+    if seed == -1:
+        seed = random.randint(0, 2**32-1)
+    np.random.seed(seed)
+    return seed
+
+
+def region_replace(
+    background: np.ndarray,
+    front: np.ndarray,
+    shift: list,
+    bmask: float = 1.0,
+    fmask: float = 1.0,
+    padded_in: bool = False,
+    padded_out: bool = False,
+    subpix: bool = False
+):
+    """
+    replace a region of the background with the front image.
+
+    Parameters
+    ----------
+    background: np.ndarray
+        The background image.
+    front: np.ndarray
+        The front image.
+    shift: list
+        The [x, y] shift of the front image. Unit: pixel.
+        Relative to the lower-left corner of the background image.
+        [0, 0] means the lower-left corner of the front image is at the lower-left corner of the background image. 
+    bmask: float
+        The mask of the background image. Default: 1.0
+        0.0 means the background image is masked.
+        1.0 means the background image is fully added.
+    fmask: float
+        The mask of the front image. Default: 1.0
+        0.0 means the front image is masked (not added).
+        1.0 means the front image is fully added.
+    padded_in: bool
+        Whether the input background image is padded. Default: False
+        In the function, the background image is padded by the size of the front image.
+        If True, means the background image is padded.
+    padded_out: bool
+        Whether the output image is padded. Default: False
+        In the function, the background image is padded by the size of the front image.
+        If True, means the output image is padded.
+        padded_in and padded_out are designed for the case that replace_region fuction is called multiple times.
+    subpix: bool
+        Whether the shift is subpixel. Default: False
+        If True, the shift is subpixel, using scipy.ndimage.shift to shift the front image.
+        If False, the shift is integer, using numpy slicing to shift the front image.
+
+    Returns
+    -------
+    np.ndarray
+        The output image.
+        shape = background.shape if padded_out = False
+        shape = background.shape + 2 * front.shape if padded_out = True
+    """
+
+    int_shift = np.array(shift).astype(int)
+    b_sz = np.array(background.shape)
+    f_sz = np.array(front.shape)
+
+    if padded_in:
+        padded = background
+        b_sz = b_sz - f_sz * 2
+    else:
+        padded = np.pad(background, ((f_sz[0], f_sz[0]), (f_sz[1], f_sz[1])))
+
+    if np.any((int_shift < -b_sz) | (int_shift > b_sz)):
+        if padded_out:
+            return padded
+        return background
+
+    if subpix:
+        subs = np.array(shift) - int_shift
+        front = nd.shift(front, (subs[0], subs[1]))
+
+    int_shift += f_sz
+    roi_y = int_shift[1]
+    roi_x = int_shift[0]
+    padded[roi_y: roi_y+f_sz[0], roi_x:roi_x+f_sz[1]] *= bmask
+    padded[roi_y: roi_y+f_sz[0], roi_x:roi_x+f_sz[1]] += fmask * front
+
+    if padded_out:
+        return padded
+
+    return padded[f_sz[0]:b_sz[0]+f_sz[0], f_sz[1]:b_sz[1]+f_sz[1]]

docs/Makefile

+# Minimal makefile for Sphinx documentation
+#
+
+# You can set these variables from the command line, and also
+# from the environment for the first two.
+SPHINXOPTS    ?=
+SPHINXBUILD   ?= sphinx-build
+SOURCEDIR     = source
+BUILDDIR      = build
+
+# Put it first so that "make" without argument is like "make help".
+help:
+	@$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
+
+.PHONY: help Makefile
+
+# Catch-all target: route all unknown targets to Sphinx using the new
+# "make mode" option.  $(O) is meant as a shortcut for $(SPHINXOPTS).
+%: Makefile
+	@$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)

docs/conf.py

+# Configuration file for the Sphinx documentation builder.
+#
+# For the full list of built-in configuration values, see the documentation:
+# https://www.sphinx-doc.org/en/master/usage/configuration.html
+
+# -- Project information -----------------------------------------------------
+# https://www.sphinx-doc.org/en/master/usage/configuration.html#project-information
+
+import os
+import sys
+
+project = 'CpicImgSim'
+copyright = '2023, CSST/CPIC Group'
+author = 'CPIC Group'
+release = '1.0.0'
+
+# -- General configuration ---------------------------------------------------
+# https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration
+
+extensions = [
+    'sphinx.ext.autodoc',
+    'sphinx.ext.napoleon',
+    'sphinx.ext.doctest',
+    'sphinx.ext.intersphinx',
+    'sphinx.ext.todo',
+    'sphinx.ext.coverage',
+    'sphinx.ext.mathjax',
+    'nbsphinx'
+]
+
+templates_path = ['_templates']
+exclude_patterns = []
+
+
+# -- Options for HTML output -------------------------------------------------
+# https://www.sphinx-doc.org/en/master/usage/configuration.html#options-for-html-output
+
+html_theme = 'sphinx_rtd_theme'
+# html_static_path = ['_static']

docs/cpism.rst

+CpicImgSim package
+==================
+
+Submodules
+----------
+
+CpicImgSim.main module
+----------------------
+
+.. automodule:: CpicImgSim.main
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+CpicImgSim.target module
+------------------------
+
+.. automodule:: CpicImgSim.target
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+CpicImgSim.optics module
+------------------------
+
+.. automodule:: CpicImgSim.optics
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+CpicImgSim.camera module
+------------------------
+
+.. automodule:: CpicImgSim.camera
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+CpicImgSim.io module
+--------------------
+
+.. automodule:: CpicImgSim.io
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+CpicImgSim.config module
+------------------------
+
+.. automodule:: CpicImgSim.config
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+CpicImgSim.utils module
+-----------------------
+
+.. automodule:: CpicImgSim.utils
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+Module contents
+---------------
+
+.. automodule:: CpicImgSim
+   :members:
+   :undoc-members:
+   :show-inheritance:

docs/index.rst

+.. CpicImgSim documentation master file, created by
+   sphinx-quickstart on Wed May 10 02:19:25 2023.
+   You can adapt this file completely to your liking, but it should at least
+   contain the root `toctree` directive.
+
+
+
+CPISM: CSST/CPI-C 观测图像仿真软件包
+=====================================================
+
+.. image:: ../logo.png
+   :scale: 40
+   :align: center
+   :alt: CPISM logo
+
+
+.. toctree::
+   :maxdepth: 1
+   :hidden:
+
+   程序安装 <installation>
+   使用说明 <tutorials>
+   代码说明 <cpism>
+
+
+
+Indices and tables
+==================
+
+* :ref:`genindex`
+* :ref:`modindex`
+* :ref:`search`

docs/installation.rst

+CPISM程序安装说明
+=================
+
+Python版本
+--------------
+
+Python >= 3.9
+
+使用setup.py安装
+-----------------
+
+联系作者获取软件安装文件，解压缩并执行以下命令：
+
+.. code-block:: bash 
+
+	cd cpic-img-sim
+	python setup.py install 
+
+设置环境变量 
+--------------
+
+Linux 或 Mac 用户
+```````````````````
+执行以下命令，打开bash配置文件：
+
+.. code-block:: bash
+
+	vi ~/.bash_profile
+
+在文件中增加如下一行：
+
+.. code-block:: bash
+
+	export cpism_refdata="/path/to/cpism/reference/"
+
+编辑bash profile文件后，需要source一下，或者重新打开一个terminal。
+
+.. code-block:: bash
+
+	source ~/.bash_profile
+
+可以用如下命令检查是否设置成功：
+
+.. code-block:: bash
+
+	echo $cpism_refdata
+	cd $cpism_refdata
+	ls
+
+Windows 用户
+`````````````
+执行以下操作：
+
+- 打开系统控制面板。
+- 选择“系统” 。
+- 选择“高级系统设置”。
+- 转到“ 高级 ”选项卡。
+- 选择 “环境变量...”。
+- 点击用户变量下的“新建...”。
+- 变量名为“cpism_refdata”，变量值为cpism_refdata文件夹的完整地址。
+- 连续点击“确定”，直到所以对话框都关闭。
+
+
+下载和安装 Pysynphot 恒星数据
+----------------------------------------
+
+CPISM中使用了 `Pysynphot package <https://pysynphot.readthedocs.io/en/latest/appendixa.html>`_ 来处理恒星光谱。为了得到恒星光谱，需要下载Pysynphot的恒星光谱数据。
+``CPISM`` 中使用了 Castelli-Kurucz 模型: `ck04models <https://archive.stsci.edu/hlsps/reference-atlases/cdbs/grid/ck04models/>`_. 
+
+可以使用如下命令进行安装
+
+.. code-block:: bash
+
+	wget http://ssb.stsci.edu/trds/tarfiles/synphot3.tar.gz
+
+下载完成后并解压，文件夹结果如下所示, 需要将trd文件夹加入到环境变量中。
+
+ ``<path>/grp/redcat/trds/grid/ck04models``
+
+Linux 或者 Mac 用户 
+````````````````````````````````
+
+.. code-block:: bash
+
+	vi ~/.bash_profile
+
+增加如下一行
+
+.. code-block:: bash
+
+	export PYSYN_CDBS="<your_path>/grp/redcat/trds"
+
+记得要source一下，或者重新打开一个terminal。 
+
+.. code-block:: bash
+
+	source ~/.bash_profile
+
+现在可以检查一下是否设置成功
+
+.. code-block:: bash
+
+	cd $PYSYN_CDBS
+	ls
+
+Windows 用户
+````````````````
+执行以下操作：
+
+- 打开系统控制面板。
+- 选择“系统” 。
+- 选择“高级系统设置”。
+- 转到“ 高级 ”选项卡。
+- 选择 “环境变量...”。
+- 点击用户变量下的“新建...”。
+- 变量名为“PYSYN_CDBS”，变量值为trds文件夹的完整地址。
+- 连续点击“确定”，直到所以对话框都关闭。
+
+程序运行配置
+------------------
+打开cpism_refdata文件夹，编辑cpism_config.yaml文件。可以进行cpism程序的运行配置。
+
+配置输出文件夹
+````````````````
+使用output_dir参数和relative_path参数配置初始文件夹的位置。
+
+- 当relative_path设置为True时，可以通过更改output_dir来更改输出文件夹的相对地址。例如，如果你想将输出文件夹放在cpism_refdata的父文件夹下的output目录中，你可以将output_dir设置为output。
+
+- 当relative_path设置为False时，可以通过设置output_dir来更改输出文件夹的绝对地址。
+
+其他配置
+````````````````
+- log_level: 日志的输出级别。可选择DEBUG、INFO、WARNING、ERROR、CRITICAL。日志文件存储在output文件夹的LOG目录下。
+- tmp_dir: 快速仿真结果的输出文件夹。默认为TMP，存储在output目录的TMP文件夹下。可选择current, 将文件存储于当前文件夹下。
+- check_fits_header: 输出文件时检查fits文件头的关键字是否和data model规定的一致，不一致时日志中输出warning信息。默认为False。
+
+
+测试是否安装成功
+------------------
+执行example.py 文件
+
+.. code-block:: bash
+
+	python example.py
+
+如果安装成功，会在output目录下SCI、CAL、TMP目录下分别生成科学观测图像文件、定标文件与快速模拟结果。文件结构详见星冕仪0级数据说明。
\ No newline at end of file

docs/make.bat

+@ECHO OFF
+
+pushd %~dp0
+
+REM Command file for Sphinx documentation
+
+if "%SPHINXBUILD%" == "" (
+	set SPHINXBUILD=sphinx-build
+)
+set SOURCEDIR=.
+set BUILDDIR=build
+
+%SPHINXBUILD% >NUL 2>NUL
+if errorlevel 9009 (
+	echo.
+	echo.The 'sphinx-build' command was not found. Make sure you have Sphinx
+	echo.installed, then set the SPHINXBUILD environment variable to point
+	echo.to the full path of the 'sphinx-build' executable. Alternatively you
+	echo.may add the Sphinx directory to PATH.
+	echo.
+	echo.If you don't have Sphinx installed, grab it from
+	echo.https://www.sphinx-doc.org/
+	exit /b 1
+)
+
+if "%1" == "" goto help
+
+%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O%
+goto end
+
+:help
+%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O%
+
+:end
+popd