Effects.py

import galsim
from matplotlib.pyplot import flag
import numpy as np
from numpy.core.fromnumeric import mean, size
from numpy.random import Generator, PCG64
import math
from numba import jit


def AddOverscan(GSImage, overscan=1000, gain=1, widthl=27, widthr=27, widtht=8, widthb=8, read_noise=5):
    """
    Add overscan/gain; gain=e-/ADU
    widthl: left pre-scan width
    widthr: right pre-scan width
    widtht: top over-scan width (the top of nd-array with small row-index)
    widthb: bottom over-scan width (the bottom of nd-array with large row-index)
    """
    imgshape = GSImage.array.shape
    newimg = galsim.Image(imgshape[1]+widthl+widthr, imgshape[0]+widtht+widthb, init_value=0)
    rng = galsim.UniformDeviate()
    NoiseOS = galsim.GaussianNoise(rng, sigma=read_noise)
    newimg.addNoise(NoiseOS)
    newimg = (newimg+overscan)/gain
    newimg.array[widtht:(widtht+imgshape[0]),widthl:(widthl+imgshape[1])] = GSImage.array
    newimg.wcs = GSImage.wcs
    # if GSImage.wcs is not None:
    #     newwcs = GSImage.wcs.withOrigin(galsim.PositionD(widthl,widtht))
    #     newimg.wcs = newwcs
    # else:
    #     pass
    return newimg


def DefectivePixels(GSImage, IfHotPix=True, IfDeadPix=True, fraction=1E-4, seed=20210304, biaslevel=500):
    # Also called bad pixels, including hot pixels and dead pixels
    # Hot Pixel > 20e-/s
    # Dead Pixel < 70%*Mean
    rgf = Generator(PCG64(int(seed*1.1)))
    if IfHotPix==True and IfDeadPix==True:
        HotFraction = rgf.random()             # fraction in total bad pixels
    elif IfHotPix==False and IfDeadPix==False:
        return GSImage
    elif IfHotPix==True:
        HotFraction = 1
    else:
        HotFraction = 0

    NPix = GSImage.array.size
    NPixBad = int(NPix*fraction)
    NPixHot = int(NPix*fraction*HotFraction)
    NPixDead = NPixBad-NPixHot
    
    NPix_y,NPix_x = GSImage.array.shape
    mean = np.mean(GSImage.array)
    rgp = Generator(PCG64(int(seed)))
    yxposfrac = rgp.random((NPixBad,2))
    YPositHot = np.array(NPix_y*yxposfrac[0:NPixHot,0]).astype(np.int32)
    XPositHot = np.array(NPix_x*yxposfrac[0:NPixHot,1]).astype(np.int32)
    YPositDead = np.array(NPix_y*yxposfrac[NPixHot:,0]).astype(np.int32)
    XPositDead = np.array(NPix_x*yxposfrac[NPixHot:,1]).astype(np.int32)

    rgh = Generator(PCG64(int(seed*1.2)))
    rgd = Generator(PCG64(int(seed*1.3)))
    if IfHotPix==True:
        GSImage.array[YPositHot,XPositHot] += rgh.gamma(2,25*150,size=NPixHot)
    if IfDeadPix==True:
        GSImage.array[YPositDead,XPositDead] = rgd.random(NPixDead)*(mean-biaslevel)*0.7+biaslevel+rgp.standard_normal()*5
    return GSImage


def BadColumns(GSImage, seed=20210309, chipid=1):
    # Set bad column values
    ysize,xsize = GSImage.array.shape
    subarr = GSImage.array[int(ysize*0.1):int(ysize*0.12), int(xsize*0.1):int(xsize*0.12)]
    meanimg = np.median(subarr)
    stdimg = np.std(subarr)
    seed += chipid
    rgn = Generator(PCG64(int(seed)))
    rgcollen = Generator(PCG64(int(seed*1.1)))
    rgxpos = Generator(PCG64(int(seed*1.2)))
    rgdn = Generator(PCG64(int(seed*1.3)))

    nbadsecA,nbadsecD = rgn.integers(low=1, high=5, size=2) 
    collen = rgcollen.integers(low=int(ysize*0.1), high=int(ysize*0.7), size=(nbadsecA+nbadsecD)) 
    xposit = rgxpos.integers(low=int(xsize*0.05), high=int(xsize*0.95), size=(nbadsecA+nbadsecD))
    signs = 2*rgdn.integers(0,2,size=(nbadsecA+nbadsecD))-1
    dn = rgdn.integers(low=meanimg*0.4, high=meanimg*0.8, size=(nbadsecA+nbadsecD))*signs
    for badcoli in range(nbadsecA):
        GSImage.array[(ysize-collen[badcoli]):ysize,xposit[badcoli]:(xposit[badcoli]+1)] += np.random.random((collen[badcoli],1))*stdimg*3+dn[badcoli]
    for badcoli in range(nbadsecD):
        GSImage.array[0:collen[badcoli+nbadsecA],xposit[badcoli+nbadsecA]:(xposit[badcoli+nbadsecA]+1)] += np.random.random((collen[badcoli+nbadsecA],1))*stdimg*3+dn[badcoli+nbadsecA]
    return GSImage


def AddBiasNonUniform16(GSImage, bias_level = 500, nsecy = 2, nsecx=8, seed=202102):
    # Generate Bias and its non-uniformity, and add the 16 bias values to the GS-Image
    rg = Generator(PCG64(int(seed)))
    Random16 = (rg.random(nsecy*nsecx)-0.5)*20
    if int(bias_level)==0:
        BiasLevel = np.zeros((nsecy,nsecx))
    elif bias_level>0:
        BiasLevel = Random16.reshape((nsecy,nsecx)) + bias_level
    arrshape = GSImage.array.shape
    secsize_x = int(arrshape[1]/nsecx)
    secsize_y = int(arrshape[0]/nsecy)
    for rowi in range(nsecy):
        for coli in range(nsecx):
            GSImage.array[rowi*secsize_y:(rowi+1)*secsize_y,coli*secsize_x:(coli+1)*secsize_x] += BiasLevel[rowi,coli]
    return GSImage


def MakeBiasNcomb(npix_x, npix_y, bias_level=500, ncombine=1, read_noise=5, gain=1, seed=202102):
    # Start with 0 value bias GS-Image
    ncombine=int(ncombine)
    BiasSngImg0 = galsim.Image(npix_x, npix_y, init_value=0)
    BiasSngImg = AddBiasNonUniform16(BiasSngImg0, 
                bias_level=bias_level, 
                nsecy = 2, nsecx=8, 
                seed=int(seed))
    BiasCombImg = BiasSngImg*ncombine
    rng = galsim.UniformDeviate()
    NoiseBias = galsim.GaussianNoise(rng=rng, sigma=read_noise*ncombine**0.5)
    BiasCombImg.addNoise(NoiseBias)
    if ncombine == 1:
        BiasTag = 'Single'
        pass
    elif ncombine >1:
        BiasCombImg /= ncombine
        BiasTag = 'Combine'
    # BiasCombImg.replaceNegative(replace_value=0)
    # BiasCombImg.quantize()
    return BiasCombImg, BiasTag


def ApplyGainNonUniform16(GSImage, gain=1, nsecy = 2, nsecx=8, seed=202102):
    # Generate Gain non-uniformity, and multipy the different factors (mean~1 with sigma~1%) to the GS-Image
    rg = Generator(PCG64(int(seed)))
    Random16 = (rg.random(nsecy*nsecx)-0.5)*0.04+1   # sigma~1%
    Gain16 = Random16.reshape((nsecy,nsecx))/gain
    print("Gain of 16 channels: ",Gain16)
    arrshape = GSImage.array.shape
    secsize_x = int(arrshape[1]/nsecx)
    secsize_y = int(arrshape[0]/nsecy)
    for rowi in range(nsecy):
        for coli in range(nsecx):
            GSImage.array[rowi*secsize_y:(rowi+1)*secsize_y,coli*secsize_x:(coli+1)*secsize_x] *= Gain16[rowi,coli]
    return GSImage


def MakeFlatSmooth(GSBounds, seed):
    rg = Generator(PCG64(int(seed)))
    r1,r2,r3,r4 = rg.random(4)
    a1 = -0.5 + 0.2*r1
    a2 = -0.5 + 0.2*r2
    a3 = r3+5
    a4 = r4+5
    xmin,xmax,ymin,ymax = GSBounds.getXMin(), GSBounds.getXMax(), GSBounds.getYMin(), GSBounds.getYMax()
    Flty, Fltx = np.mgrid[ymin:(ymax+1), xmin:(xmax+1)]
    rg = Generator(PCG64(int(seed)))
    p1,p2,bg=rg.poisson(1000, 3)
    Fltz = 1e-6*(a1 * (Fltx-p1) ** 2 + a2 * (Flty-p2) ** 2 - a3*Fltx - a4*Flty) + bg*20
    FlatImg = galsim.ImageF(Fltz)
    return FlatImg


def MakeFlatNcomb(flat_single_image, ncombine=1, read_noise=5, gain=1, overscan=500, biaslevel=500, seed_bias=20210311):
    ncombine=int(ncombine)
    FlatCombImg = flat_single_image*ncombine
    rng = galsim.UniformDeviate()
    NoiseFlatPoi = galsim.PoissonNoise(rng=rng, sky_level=0)
    FlatCombImg.addNoise(NoiseFlatPoi)
    NoiseFlatReadN = galsim.GaussianNoise(rng=rng, sigma=read_noise*ncombine**0.5)
    FlatCombImg.addNoise(NoiseFlatReadN)
    # NoiseFlat = galsim.CCDNoise(rng, gain=gain, read_noise=read_noise*ncombine**0.5, sky_level=0)
    for i in range(ncombine):
        FlatCombImg = AddBiasNonUniform16(
            FlatCombImg, 
            bias_level=biaslevel, 
            nsecy=2, nsecx=8, 
            seed=seed_bias)
    if ncombine == 1:
        FlatTag = 'Single'
        pass
    elif ncombine >1:
        FlatCombImg /= ncombine
        FlatTag = 'Combine'
    # FlatCombImg.replaceNegative(replace_value=0)
    # FlatCombImg.quantize()
    return FlatCombImg, FlatTag


def MakeDarkNcomb(npix_x, npix_y, overscan=500, bias_level=500, seed_bias=202102, darkpsec=0.02, exptime=150, ncombine=10, read_noise=5, gain=1):
    ncombine=int(ncombine)
    darkpix = darkpsec*exptime
    DarkSngImg = galsim.Image(npix_x, npix_y, init_value=darkpix)
    rng = galsim.UniformDeviate()
    NoiseDarkPoi = galsim.PoissonNoise(rng=rng, sky_level=0)
    NoiseReadN = galsim.GaussianNoise(rng=rng, sigma=read_noise*ncombine**0.5)
    DarkCombImg = DarkSngImg*ncombine
    DarkCombImg.addNoise(NoiseDarkPoi)
    DarkCombImg.addNoise(NoiseReadN)
    for i in range(ncombine):
        DarkCombImg = AddBiasNonUniform16(
            DarkCombImg, 
            bias_level=bias_level, 
            nsecy = 2, nsecx=8, 
            seed=int(seed_bias))
    if ncombine == 1:
        DarkTag = 'Single'
        pass
    elif ncombine >1:
        DarkCombImg /= ncombine
        DarkTag = 'Combine'
    # DarkCombImg.replaceNegative(replace_value=0)
    # DarkCombImg.quantize()
    return DarkCombImg, DarkTag


def PRNU_Img(xsize, ysize, sigma=0.01, seed=202101):
    rg = Generator(PCG64(int(seed)))
    prnuarr = rg.normal(1, sigma, (ysize,xsize))
    prnuimg = galsim.ImageF(prnuarr)
    return prnuimg


def NonLinearity(GSImage, beta1=1E-7, beta2=1E-10):
    NonLinear_f = lambda x, beta_1, beta_2: x - beta_1*x*x + beta_2*x*x*x
    GSImage.applyNonlinearity(NonLinear_f, beta1, beta2)
    return GSImage


def chargeflow(ndarr, fullwell=10E4):
    size_y,size_x = ndarr.shape
    satpos_y = np.where(ndarr>=fullwell)[0]
    satpos_x = np.where(ndarr>=fullwell)[1]
    Nsatpix = len(satpos_y)

    if Nsatpix==0:
        # make no change for the image array
        return ndarr

    else:    
        for i in range(Nsatpix):
            ClumpFullwell=True
            satcharge = ndarr[satpos_y[i],satpos_x[i]]-fullwell
            ndarr[satpos_y[i],satpos_x[i]] = fullwell
            
            satpos_yi0 = satpos_y[i]

            # Define the x,y=0,0 element of image array is the lower-left corner. So y decreases being downward.

            # print('Charge Clump moves down')
            chargedn = ((np.random.random()-0.5)*0.05+0.5)*satcharge
            chargeup = satcharge - chargedn
            fwi = 1
            aa = np.log(chargedn/fullwell)**3*0.9   # blooming length begin to has e- less than fullwell
            if aa < 0.05:
                ClumpFullwell=False
            # Test
            try:
                while ClumpFullwell==True:
                    if satpos_y[i]<=0:
                        break
                    if ndarr[satpos_y[i]-1,satpos_x[i]]<fullwell:
                        ndarr[satpos_y[i]-1,satpos_x[i]] = ndarr[satpos_y[i]-1,satpos_x[i]] + chargedn
                        if ndarr[satpos_y[i]-1,satpos_x[i]]>=fullwell:
                            fx = 0.5*(math.exp(math.log(fwi)**3/aa)+np.exp(-1*math.log(fwi)**3/aa))
                            if fx>5:
                                fx=5
                            faa= 0.5*(math.exp(aa/aa)+np.exp(-1*aa/aa))
                            rand_frac = 1-0.1*(fx-1)/(faa-1)

                            chargedn = ndarr[satpos_y[i]-1,satpos_x[i]] - fullwell
                            ndarr[satpos_y[i]-1,satpos_x[i]] = fullwell*rand_frac
                            satpos_y[i] = satpos_y[i]-1
                            if satpos_y[i]<0:
                                ClumpFullwell=False
                                break
                        else:
                            ClumpFullwell=False
                        fwi += 1
                    else:
                        satpos_y[i] = satpos_y[i]-1
                        if satpos_y[i]<0:
                            ClumpFullwell=False
                            break
                  
                
                # print('Charge Clump moves up')
                ClumpFullwell=True
                satpos_y[i] = satpos_yi0
                fwi = 1
                aa = np.log(chargeup/fullwell)**3*0.9   # blooming length at which it begins to have e- less than fullwell
                if aa < 0.05:
                    ClumpFullwell=False
                while ClumpFullwell==True:
                    if satpos_y[i]>=size_y-1:
                        break
                    if ndarr[satpos_y[i]+1,satpos_x[i]]<fullwell:
                        ndarr[satpos_y[i]+1,satpos_x[i]] = ndarr[satpos_y[i]+1,satpos_x[i]] + chargeup
                        if ndarr[satpos_y[i]+1,satpos_x[i]]>=fullwell:
                            fx = 0.5*(math.exp(math.log(fwi)**3/aa)+np.exp(-1*math.log(fwi)**3/aa))
                            if fx>5:
                                fx=5
                            faa= 0.5*(math.exp(aa/aa)+np.exp(-1*aa/aa))
                            rand_frac = 1-0.1*(fx-1)/(faa-1)
                            chargeup = ndarr[satpos_y[i]+1,satpos_x[i]] - fullwell
                            ndarr[satpos_y[i]+1,satpos_x[i]] = fullwell*rand_frac
                            satpos_y[i] = satpos_y[i]+1
                            if satpos_y[i]>=size_y-1:
                                ClumpFullwell=False
                                break
                        else:
                            ClumpFullwell=False
                        fwi += 1
                    else:
                        satpos_y[i] = satpos_y[i]+1
                        if satpos_y[i]>=size_y-1:
                            ClumpFullwell=False
                            break
            except Exception as e:
                print(e)
                print(fwi, aa)
                pass
          
    return ndarr


def SaturBloom(GSImage, nsect_x=1, nsect_y=1, fullwell=10e4):
    """
    To simulate digital detector's saturation and blooming effect. The blooming is along the read-out direction, perpendicular to the charge transfer direction. Charge clumpy overflows the pixel well will flow to two oposite directions with nearly same charges.
    Work together with chargeflow() function.
    Parameters:
      GSImage: a GalSim Image of a chip;
      nsect_x: number of sections divided along x direction;
      nsect_y: number of sections divided along y direction.
    Return: a saturated image array with the same size of input.
    """
    imgarr = GSImage.array
    size_y, size_x = imgarr.shape
    subsize_y = int(size_y/nsect_y)
    subsize_x = int(size_x/nsect_x)

    for i in range(nsect_y):
        for j in range(nsect_x):
            subimg = imgarr[subsize_y*i:subsize_y*(i+1), subsize_x*j:subsize_x*(j+1)]

            subimg = chargeflow(subimg, fullwell=fullwell)

            imgarr[subsize_y*i:subsize_y*(i+1), subsize_x*j:subsize_x*(j+1)] = subimg

    return GSImage


def readout16(GSImage, rowi=0, coli=0, overscan_value=0):
    # readout image as 16 outputs of sub-images plus prescan & overscan.
    # assuming image width and height sizes are both even.
    # assuming image has 2 columns and 8 rows of output channels. 
    # 00  01
    # 10  11
    # 20  21
    # ...
    # return: GS Image Object
    npix_y,npix_x = GSImage.array.shape
    subheight = int(8+npix_y/2+8)
    subwidth = int(16+npix_x/8+27)
    OutputSubimg = galsim.ImageUS(subwidth, subheight, init_value=overscan_value)
    if rowi<4 and coli==0:
        subbounds = galsim.BoundsI(1, int(npix_x/2),  int(npix_y/8*rowi+1), int(npix_y/8*(rowi+1)))
        subbounds = subbounds.shift(galsim.PositionI(GSImage.bounds.getXMin()-1, GSImage.bounds.getYMin()-1))
        subimg = GSImage[subbounds]
        OutputSubimg.array[27:int(npix_y/8)+27,8:int(npix_x/2)+8] = subimg.array
    elif rowi<4 and coli==1:
        subbounds = galsim.BoundsI(npix_x/2+1, npix_x,  npix_y/8*rowi+1, npix_y/8*(rowi+1))
        subbounds = subbounds.shift(galsim.PositionI(GSImage.bounds.getXMin()-1, GSImage.bounds.getYMin()-1))
        subimg = GSImage[subbounds]
        OutputSubimg.array[27:int(npix_y/8)+27,8:int(npix_x/2)+8] = subimg.array
    elif rowi>=4 and rowi<8 and coli==0:
        subbounds = galsim.BoundsI(1, npix_x/2,  npix_y/8*rowi+1, npix_y/8*(rowi+1))
        subbounds = subbounds.shift(galsim.PositionI(GSImage.bounds.getXMin()-1, GSImage.bounds.getYMin()-1))
        subimg = GSImage[subbounds]
        OutputSubimg.array[16:int(npix_y/8)+16,8:int(npix_x/2)+8] = subimg.array
    elif rowi>=4 and rowi<8 and coli==1:
        subbounds = galsim.BoundsI(npix_x/2+1, npix_x,  npix_y/8*rowi+1, npix_y/8*(rowi+1))
        subbounds = subbounds.shift(galsim.PositionI(GSImage.bounds.getXMin()-1, GSImage.bounds.getYMin()-1))
        subimg = GSImage[subbounds]
        OutputSubimg.array[16 :int(npix_y/8)+16,8:int(npix_x/2)+8] = subimg.array
    else:
        print("\n\033[31mError: "+"Wrong rowi or coli assignment. Permitted: 0<=rowi<=7, 0<=coli<=1."+"\033[0m\n")
        return OutputSubimg
    return OutputSubimg


def CTE_Effect(GSImage, threshold=27, direction='column'):
    # Devide the image into 4 sections and apply CTE effect with different trail directions.
    # GSImage: a GalSim Image object.
    size_y,size_x = GSImage.array.shape
    size_sect_y = int(size_y/2)
    size_sect_x = int(size_x/2)
    imgarr = GSImage.array
    if direction == 'column':
        imgarr[0:size_sect_y,:]      = CTEModelColRow(imgarr[0:size_sect_y,:],      trail_direction='down', direction='column', threshold=threshold)
        imgarr[size_sect_y:size_y,:] = CTEModelColRow(imgarr[size_sect_y:size_y,:], trail_direction='up',   direction='column', threshold=threshold)
    elif direction == 'row':
        imgarr[:,0:size_sect_x]      = CTEModelColRow(imgarr[:,0:size_sect_x],      trail_direction='right', direction='row', threshold=threshold)
        imgarr[:,size_sect_x:size_x] = CTEModelColRow(imgarr[:,size_sect_x:size_x], trail_direction='left',  direction='row', threshold=threshold)
    return GSImage


@jit()
def CTEModelColRow(img, trail_direction = 'up', direction='column', threshold=27):

    #total trail flux vs (pixel flux)^1/2 is approximately linear
    #total trail flux = trail_a * (pixel flux)^1/2 + trail_b
    #trail pixel flux = pow(0.5,x)/0.5, normalize to 1
    trail_a =  5.651803799619966
    trail_b = -2.671933068990729

    sh1 = img.shape[0]
    sh2 = img.shape[1]
    n_img = img*0
    idx = np.where(img<threshold)
    if len(idx[0]) == 0:
        pass
    elif len(idx[0])>0:
        n_img[idx] = img[idx]

    yidx,xidx = np.where(img>=threshold)
    if len(yidx) == 0:
        pass
    elif len(yidx)>0:
        # print(index)
        for i, j in zip(yidx,xidx):
            f = img[i,j]
            trail_f = (np.sqrt(f)*trail_a + trail_b)*0.5
            # trail_f=5E-5*f**1.5
            xy_num = 10
            all_trail = np.zeros(xy_num)
            
            xy_upstr = np.arange(1,xy_num,1)

            # all_trail_pix = np.sum(pow(0.5,xy_upstr)/0.5)
            all_trail_pix = 0
            for m in xy_upstr:
                a1=12.97059491  
                b1=0.54286652
                c1=0.69093105
                a2=2.77298856
                b2=0.11231055
                c2=-0.01038675
                # t_pow = 0
                am=1
                bm=1
                t_pow = am*np.exp(-bm*m)
                # if m < 5:
                #     t_pow = a1*np.exp(-b1*m)+c1
                # else:
                #     t_pow = a2*np.exp(-b2*m)+c2
                if t_pow <0:
                    t_pow = 0

                all_trail_pix += t_pow
                all_trail[m] = t_pow


            trail_pix_eff = trail_f/all_trail_pix

            all_trail = trail_pix_eff*all_trail

            all_trail[0] = f - trail_f
            
            for m in np.arange(0,xy_num,1):
                if direction == 'column':
                    if trail_direction == 'down':
                        y_pos = i + m
                    elif trail_direction == 'up':
                        y_pos = i - m
                    if y_pos < 0 or y_pos >=sh1:
                        break
                    n_img[y_pos,j] = n_img[y_pos,j] + all_trail[m]
                elif direction == 'row':
                    if trail_direction == 'left':
                        x_pos = j - m
                    elif trail_direction == 'right':
                        x_pos = j + m
                    if x_pos < 0 or x_pos >=sh2:
                        break
                    n_img[i,x_pos] = n_img[i,x_pos] + all_trail[m]

    return n_img



#---------- For Cosmic-Ray Simulation ------------
#---------- Zhang Xin ----------------------------

def getYValue(collection, x):
    index = 0;
    if (collection.shape[1] == 2):
        while(x>collection[index,0] and index < collection.shape[0]):
            index= index + 1;
        if (index == collection.shape[0] or index == 0):
            return 0;

        deltX = collection[index, 0] - collection[index-1, 0];
        deltY = collection[index, 1] - collection[index-1, 1];

        if deltX == 0:
            return (collection[index, 1] + collection[index-1, 1])/2.0
        else:
            a = deltY/deltX;
            return a * (x - collection[index-1,0]) + collection[index-1, 1];
    return 0;


def selectCosmicRayCollection(attachedSizes, xLen, yLen,cr_pixelRatio,CR_max_size):

    normalRay = 0.90
    nnormalRay = 1-normalRay
    max_nrayLen = 100
    pixelNum = int(xLen *  yLen * cr_pixelRatio * normalRay );
    pixelNum_n = int(xLen *  yLen * cr_pixelRatio * nnormalRay )
    CRPixelNum = 0;
    
    maxValue = max(attachedSizes[:,1])
    maxValue += 0.1;

    cr_event_num = 0;
    CRs = np.zeros(pixelNum);
    while (CRPixelNum < pixelNum):
        x = CR_max_size * np.random.random();
        y = maxValue * np.random.random();
        if (y <= getYValue(attachedSizes, x)):
            CRs[cr_event_num] = np.ceil(x);
            cr_event_num = cr_event_num + 1;
            CRPixelNum = CRPixelNum + round(x);

    while (CRPixelNum < pixelNum + pixelNum_n):
        nx = np.random.random()*(max_nrayLen-CR_max_size)+CR_max_size
        CRs[cr_event_num] = np.ceil(nx);
        cr_event_num = cr_event_num + 1;
        CRPixelNum = CRPixelNum + np.ceil(nx);

    return   CRs[0:cr_event_num];


def defineEnergyForCR(cr_event_size):
    import random
    sigma = 0.6 / 2.355;
    mean = 3.3;

    energys = np.zeros(cr_event_size);
    for i in np.arange(cr_event_size):
        energy_index = random.normalvariate(mean,sigma);
        energys[i] = pow(10, energy_index);

    return energys;

def convCR(CRmap=None, addPSF=None, sp_n = 4):
    sh = CRmap.shape

    # sp_n = 4
    subCRmap = np.zeros(np.array(sh)*sp_n)
    pix_v0 = 1/(sp_n*sp_n)
    for i in np.arange(sh[0]):
        i_st = sp_n*i
        for j in np.arange(sh[1]):
            if CRmap[i,j] ==0:
                continue
            j_st = sp_n*j
            pix_v1 = CRmap[i,j]*pix_v0
            for m in np.arange(sp_n):
                for n in np.arange(sp_n):
                    subCRmap[i_st+m, j_st + n] = pix_v1

    m_size = addPSF.shape[0]

    subCRmap_n = np.zeros(np.array(subCRmap.shape) + m_size -1)

    for i in np.arange(subCRmap.shape[0]):
        for j in np.arange(subCRmap.shape[1]):
            if subCRmap[i,j]>0:
                convPix = addPSF*subCRmap[i,j]
                subCRmap_n[i:i+m_size,j:j+m_size]+=convPix

    CRmap_n = np.zeros((np.array(subCRmap_n.shape)/sp_n).astype(np.int32))
    sh_n = CRmap_n.shape

    for i in np.arange(sh_n[0]):
        i_st = sp_n*i
        for j in np.arange(sh_n[1]):
            p_v = 0
            j_st=sp_n*j
            for m in np.arange(sp_n):
                for n in np.arange(sp_n):
                    p_v += subCRmap_n[i_st+m, j_st + n]

            CRmap_n[i,j] = p_v

    return CRmap_n


def produceCR_Map(xLen, yLen, exTime, cr_pixelRatio, gain, attachedSizes, seed=20210317):
    # Return: an 2-D numpy array
    # attachedSizes = np.loadtxt('./wfc-cr-attachpixel.dat');
    np.random.seed(seed)
    CR_max_size = 20.0;
    cr_size = selectCosmicRayCollection(attachedSizes, xLen, yLen, cr_pixelRatio, CR_max_size);

    cr_event_size = cr_size.shape[0];
    cr_energys = defineEnergyForCR(cr_event_size);

    CRmap = np.zeros([yLen, xLen]);

    ## produce conv kernel
    from astropy.modeling.models import Gaussian2D
    o_size = 4
    sp_n = 8

    m_size = o_size*sp_n+1
    m_cen = o_size*sp_n/2
    sigma_psf = 0.2*sp_n
    addPSF_ = Gaussian2D(1, m_cen, m_cen, sigma_psf, sigma_psf)
    yp, xp = np.mgrid[0:m_size, 0:m_size]
    addPSF = addPSF_(xp, yp)
    convKernel = addPSF/addPSF.sum()

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


    for i in np.arange(cr_event_size):
        xPos = round((xLen - 1)* np.random.random());
        yPos = round((yLen - 1)* np.random.random());
        cr_lens = int(cr_size[i]);
        if cr_lens ==0:
            continue;
        pix_energy = cr_energys[i]/gain/cr_lens;
        pos_angle = 1/2*math.pi*np.random.random();

        crMatrix = np.zeros([cr_lens+1, cr_lens + 1])

        for j in np.arange(cr_lens):
            x_n = int(np.cos(pos_angle)*j - np.sin(pos_angle)*0);
            if x_n < 0:
                x_n = x_n + cr_lens+1
            y_n = int(np.sin(pos_angle)*j + np.cos(pos_angle)*0);
            if x_n<0 or x_n >cr_lens or y_n < 0 or y_n > cr_lens:
                continue;
            crMatrix[y_n,x_n] = pix_energy;

        crMatrix_n = convCR(crMatrix, convKernel, sp_n)
        # crMatrix_n = crMatrix

        xpix = np.arange(crMatrix_n.shape[0]) + int(xPos)
        ypix = np.arange(crMatrix_n.shape[1]) + int(yPos)

        sh = CRmap.shape
        okx = (xpix >= 0) & (xpix < sh[1])
        oky = (ypix >= 0) & (ypix < sh[0])

        sly = slice(ypix[oky].min(), ypix[oky].max()+1)
        slx = slice(xpix[okx].min(), xpix[okx].max()+1)
        CRmap[sly, slx] += crMatrix_n[oky,:][:,okx]


    return CRmap.astype(np.int32);


def ShutterEffectArr(GSImage, t_shutter=1.3, dist_bearing=735, dt=1E-3):
    # Generate Shutter-Effect normalized image
    # t_shutter: time of shutter movement
    # dist_bearing: distance between two bearings of shutter leaves
    # dt: delta_t of sampling

    from scipy import interpolate

    SampleNumb = int(t_shutter/dt+1)
    DistHalf = dist_bearing/2

    t = np.arange(SampleNumb)*dt
    a_arr = 5.84*np.sin(2*math.pi/t_shutter*t)
    v = np.zeros(SampleNumb)
    theta = np.zeros(SampleNumb)
    x = np.arange(SampleNumb)/(SampleNumb-1)*dist_bearing
    s = np.zeros(SampleNumb)
    s1 = np.zeros(SampleNumb)
    s2 = np.zeros(SampleNumb)
    brt  = np.zeros(SampleNumb)
    idx = np.arange(SampleNumb)
    sidx = np.zeros(SampleNumb)
    s1idx = np.zeros(SampleNumb)
    s2idx = np.zeros(SampleNumb)

    v[0] = 0
    theta[0] = 0

    for i in range(SampleNumb-1):
        v[i+1] = v[i]+a_arr[i]*dt
        theta[i+1] = theta[i]+v[i]*dt
        s1[i] = DistHalf*np.cos(theta[i])
        s2[i] = dist_bearing-DistHalf*np.cos(theta[i])
        s1idx[i] = int(s1[i]/dist_bearing*(SampleNumb))
        s2idx[i] = int(s2[i]/dist_bearing*(SampleNumb))
        brt[(idx>s1idx[i]) & (idx<s2idx[i])] += dt

    brt = brt*2+(150-t_shutter*2)

    x = (x-dist_bearing/2)*100

    intp = interpolate.splrep(x, brt, s=0)

    xmin = GSImage.bounds.getXMin()
    xmax = GSImage.bounds.getXMax()
    ymin = GSImage.bounds.getYMin()
    ymax = GSImage.bounds.getYMax()
    if xmin<np.min(x) or xmax>np.max(x):
        raise LookupError("Out of focal-plane bounds in X-direction.")
    if ymin<-25331 or ymax>25331:
        raise LookupError("Out of focal-plane bounds in Y-direction.")
    sizex = xmax-xmin+1
    sizey = ymax-ymin+1
    xnewgrid = np.mgrid[xmin:(xmin+sizex)]
    expeffect = interpolate.splev(xnewgrid, intp, der=0)
    expeffect /= np.max(expeffect)
    exparrnormal = np.tile(expeffect, (sizey,1))
    # GSImage *= exparrnormal

    return exparrnormal