Pro NOISE1
    ;test how to add poisson&dark noise to an image, conv w/ IPC, add readnoise, then deconv w/ IPC
    ;all the noises to the resampled airy pattern not the ideal airy pattern

    restore,'~/Desktop/IDL/idl_IPC/In_sam.sav' ;load resampled airy pattern
    imsize=size(In_sam,/dim)
    In_sam_noi=fltarr(imsize(0),imsize(1))
    N_rd=10; 10  ;read noise per pixel IUO pixel value
    N_dark=3     ;dark current per pix
    ap1=(imsize(0)+4)/2-1 & ap2=(imsize(1)+4)/2+1;1 ;aperture size of some sampling around the central peak
    Npix=(ap2-ap1+1)^2 ;# of pixel within the aperture

    In_pad=zeropad(In_sam, imsize(0)+4,imsize(1)+4,/center) ;sampled img w/o any noise
    sig=total(In_pad(ap1:ap2,ap1:ap2))
    F0=1000; true flux level in AP
    In_sam=In_sam*(F0/sig) ;change the flux level to F0 within AP

    restore, '~/Desktop/IDL/idl_IPC/surf_ipc.sav' ; load the ipc data from dark_19&20
    h=surf(3:7,3:7);/total(surf(3:7,3:7))
    if (min(h) LT 0.) then begin
      h=h+abs(min(h)) ;make sure all the pixel value larger than zero
    endif

    ;ceate an IPC filter w/ large coupling, not the real case
    h_t=fltarr(5,5);
    h_t(1:3,1:3)=5
    h_t(2,1:3)=10 & h_t(1:3,2)=10
    h_t(2,2)=100
    h=h_t

    h=h/total(h) ;normalized to 1 over total area
    ;try to flip h because convol do not flip a filter
    h1=h ;
    for i=0,4 do begin
      for j=0,4 do begin
        h1(i,j)=h(4-i,4-j)
      endfor
    endfor

    i0=0
    num_run=20000;0; the number of iteration when adding noise
    sig1=fltarr(num_run) & sig1_conv=fltarr(num_run)
    sig2=fltarr(num_run) & sig3=fltarr(num_run)
    diff1=fltarr(num_run) & diff1_CONV=fltarr(num_run)
    diff2=fltarr(num_run) & diff3=fltarr(num_run)
    SNR1=fltarr(num_run) & SNR1_conv=fltarr(num_run)
    SNR2=fltarr(num_run) & SNR3=fltarr(num_run)
    ;diff1_conv=0 & diff2=0 & diff3=0 & diff0=0
    ;noi_rd=0

    WHILE (i0 LT num_run) DO BEGIN

       ;add poisson noise to each pixel
       for i=0,imsize(0)-1 do begin
         for j=0, imsize(1)-1 do begin
          if(In_sam(i,j) GT 0.) then begin
              ;dummy=i+j
              In_sam_noi(i,j)=randomn(s1,1,poisson=In_sam(i,j))
          endif
         endfor
       endfor

    ;Do convoluion and deconvolution with IPC h on noisy sampled airy pattern
    hp=zeropad(h1, imsize(0)+4,imsize(1)+4,/center);/total(h)
    In_pad0=zeropad(In_sam, imsize(0)+4,imsize(1)+4,/center) ;sampled img w/o any noise

    noi_dk_raw=randomn(s2,imsize(0)+5,imsize(1),poisson=N_dark) ; ceate dark current pattern
    noi_dk=noi_dk_raw(0:imsize(0)-1,*)
    dark_pad=zeropad(noi_dk, imsize(0)+4,imsize(1)+4,/center);+dark current
    dark=mean(noi_dk_raw(imsize(0):imsize(0)+4,*))
    In_pad1=zeropad(In_sam_noi+noi_dk, imsize(0)+4,imsize(1)+4,/center);+dark_pad ;sampled img w/ poisson noise&dark
    ;In_pad_poi=zeropad(In_sam_noi-In_sam, imsize(0)+4,imsize(1)+4,/center) ;only poisson noise

    ;In_conv0=convol(In_pad0, h,total(h),/center,/EDGE_TRUNCATE )
    In_conv1=convol(In_pad1, h,total(h),/center,/EDGE_TRUNCATE )
    noi_rd_raw=randomn(s1,imsize(0)+5,imsize(1))*N_rd ;+noi_rd  produce random readnoise pattern
    noi_rd=noi_rd_raw(5:imsize(0)+4,*)
    readnoise=stddev(noi_rd_raw(0:4,*)) ;sample to obtain real readnoise

    In_sam_noi2=In_conv1+zeropad(noi_rd, imsize(0)+4,imsize(1)+4,/center) ;adding readnoise to IPC convolved image

    ;In_conv_poi=convol(In_pad_poi, h,total(h),/center,/EDGE_TRUNCATE )
    ;edge_truncate: see smooth for explanation, choose neighbouring pixels for pixel out of array
    ;edge_wrap:  wrapping" the pixel subscripts out of Array, e.g. size(a)=5, a(5)<<a(0)
    ;center=0: first flip h around its 1st element h0, then multiply with data element and sum
    ;center=1, not flip h, just center h over each data element; no difference for symmetric h, flip h
    ;before apply convol for non-sysmmetry kernel(h).

    i=i0
    sig0=total(In_pad0(ap1:ap2,ap1:ap2)) ;obtain signal/flux within the preset aperture
    sig1(i)=total(In_pad1(ap1:ap2,ap1:ap2));+sig1 ;cumulative signal
    ;print, sig1(i), in_pad1(ap1:ap2,ap1:ap2)
    sig1_conv(i)=total(In_conv1(ap1:ap2,ap1:ap2));+sig1_conv
    ;diff1=sig1-sig0
    ;diff1_conv=sig1_conv-sig1+diff1_conv

    sig2(i)=total(In_sam_noi2(ap1:ap2,ap1:ap2));+sig2
    ;diff2=sig2-sig1+diff2

    ;SNR1=0 & SNR1_CONV=0 & SNR2=0 & SNR3=0
    ;SNR1=sig1/sqrt(sig1)+SNR1 ;just poisson noise added
    ;SNR1_conv=sig1_conv/sqrt(sig1_conv)+SNR1_CONV ; after convolve with IPC
    ;SNR2=sig2/sqrt(sig2+(ap2-ap1+1)^2*readnoise)+SNR2 ;after convl and readnoise added


    ;Do decovolution on the readnoise added img
    hf=fft(hp)
    Inf=fft(In_sam_noi2) ;In_sam_noi2 or In_conv1: set to check if deconv operation works
    imf_res=Inf/hf
    img_res=fft(imf_res, /inverse) ; inverse FT of deconvolved spectrum
    In_deconv=shift(img_res,(imsize(0)+4)/2,(imsize(1)+4)/2)

    In_deconv=in_deconv/((imsize(0)+4)*(imsize(1)+4)) ;correction for a normalization factor
    F_deconv=abs(total(in_deconv(ap1:ap2,ap1:ap2))) ;total flux after deconvolving the noisy image
                                               ;has to take ABS after total b/c negative pix value by RD
    sig3(i)=F_deconv;+sig3
    ;SNR3=F_deconv/sqrt(F_deconv+(ap2-ap1+1)^2*readnoise)+SNR3 ;after decconvl
    ;diff3=F_deconv-sig1+diff3

    ;calculate flux changes relative to the signal just w/ poisson noise
    diff1(i)=(sig1(i)-Npix*dark-sig0)/sig0*100  ;@%
    diff1_conv(i)=(sig1_conv(i)-sig1(i))/(sig1(i)-Npix*dark)*100
    diff2(i)=(sig2(i)-sig1(i))/(sig1(i)-Npix*dark)*100
    diff3(i)=(sig3(i)-sig1(i))/(sig1(i)-Npix*dark)*100

    ;calculate SNR after adding poisson noise, IPC convl, then read noise, lastly deconv with IPC
    SNR1(i)=(sig1(i)-Npix*dark)/sqrt(sig1(i));+Npix*dark) ;just poisson noise added
    ;print, (mean(sig1)-dark*Npix)/sqrt(mean(sig1))
    SNR1_conv(i)=(sig1_conv(i)-Npix*dark)/sqrt(sig1_conv(i));+Npix*dark) ;after convolve with IPC
    ;print, (mean(sig1_conv)-dark*Npix)/sqrt(mean(sig1_conv))
    SNR2(i)=(sig2(i)-Npix*dark)/sqrt(sig2(i)+Npix*readnoise^2) ; +Npix*dark;after convl and readnoise added
    ;print, (mean(sig2)-dark*Npix)/sqrt(mean(sig2)+Npix*readnoise^2)

    SNR3(i)=(sig3(i)-Npix*dark)/sqrt(sig3(i)+Npix*readnoise^2) ; +Npix*dark;after deconvolved with IPC
    ;print, (mean(sig3)-dark*Npix)/sqrt(mean(sig3)+Npix*readnoise^2)

    ;print, in_pad1(ap1:ap2,ap1:ap2), 's',abs(in_deconv(ap1:ap2,ap1:ap2)),'s'
    ;try to check the deconv operation if ''Inf=fft(In_conv1)''
    ;print, in_sam_noi2(ap1:ap2,ap1:ap2)
    ;print, sig0,total(In_pad1(ap1:ap2,ap1:ap2)),total(In_conv1(ap1:ap2,ap1:ap2)),$
    ;      total(In_sam_noi2(ap1:ap2,ap1:ap2)), F_deconv

    i0=i0+1
    endwhile

    ;taking average of signal
    ;SIG0: TRUE signal w/o noise
    ;sig1: signal with poisson noise; sig1_conv: sig1 convolved with IPC
    ;sig2: sig1_conv with read noise; sig3: sig2 deconvolved with IPC
    ;sig1=sig1/num_run & sig1_conv=sig1_conv/num_run
    ;sig2=sig2/num_run & sig3=sig3/num_run


    Print, 'The flux level over AP & # of run:', sig0, num_run
    print, 'The relative changes of total flux after poisson, conv, RD, deconv: ', $
            mean(diff1), mean(diff1_CONV), mean(diff2), mean(diff3)
    print,'SNR w/ poisson noise, ipc convl, readnoise, deconv', mean(SNR1), $
                            mean(SNR1_conv), mean(SNR2), mean(SNR3)
    print, (mean(sig1)-dark*Npix)/sqrt(mean(sig1)),(mean(sig1_conv)-dark*Npix)/sqrt(mean(sig1_conv)), $
       (mean(sig2)-dark*Npix)/sqrt(mean(sig2)+Npix*readnoise^2), (mean(sig3)-dark*Npix)/sqrt(mean(sig3)+Npix*readnoise^2)

    stop

    end