astropy / ccd-reduction-and-photometry-guide

Read the CCD guide here:
http://www.astropy.org/ccd-reduction-and-photometry-guide/
BSD 3-Clause "New" or "Revised" License
101 stars 46 forks source link

the simulation of dark current in section 1.2 seems problematic #338

Open wa1kmoon opened 1 year ago

wa1kmoon commented 1 year ago

Hi, I'm a beginner on ccd-reduction and I really appreciate this friendly guide!

However, when I read the dark_current function defined in section 1.2, I find it difficult to understand following code:

    # dark current for every pixel; we'll modify the current for some pixels if 
    # the user wants hot pixels.
    base_current = current * exposure_time / gain

    # This random number generation should change on each call.
    dark_im = noise_rng.poisson(base_current, size=image.shape)

I think it is the current * exposure_time that obeys the poisson distribution, so in my opinion, gain should be divided after the random image(in electrons) is generated(like how the function sky_background does). So this code should be like below:

    # dark current for every pixel; we'll modify the current for some pixels if 
    # the user wants hot pixels.
    base_current = current * exposure_time

    # This random number generation should change on each call.
    dark_im = noise_rng.poisson(base_current, size=image.shape) / gain

I'm not very sure which one is right, so I made an 'experiment' using the code snippet below:

dark_cur = 1
exposure = 100
gain = 2
imgsize = (100,100)

seed = os.getenv('GUIDE_RANDOM_SEED', None)
rng = np.random.default_rng(seed)

darkimg_exp = rng.poisson(dark_cur * exposure / gain, imgsize) # the one used in `dark_current`
darkimg_mul = rng.poisson(dark_cur * exposure, imgsize) / gain # the one I supposed
darkimg_cum = rng.zeros_like(darkimg_exp) # simulate a cumulative process
for i in range(exposure):
    darkimg_cum = darkimg_cum + rng.poisson(dark_cur, imgsize) / gain

print(np.std(darkimg_exp),np.std(darkimg_mul),np.std(darkimg_cum))

then I got (7.118792059752835, 5.0494683244872425, 4.946861679640941), which means that probably I'm right?