;-----------------------------------------------------------------------------
FUNCTION fg_noise, im, DIAMETER=diameter, LAMBDA=lambda, PLATE_SCALE=plate_scale, $
                     FILTIM=filtim, $
                     QUIET=quiet
;+
; NAME:
;	FG_noise
;
; PURPOSE:
;	Estimate the photometric noise of a BFI or NFI filtergram. 
;
; CATEGORY:
;	Image processing
;
; SAMPLE CALLING SEQUENCE:
;       noise = fg_noise(image)
;       noise = fg_noise(image, FILTIM=fim)    where fim = image*0, e.g.
;
; INPUTS:
;	IM         - A 2-D image.
; 
; OUTPUTS: 
;       NOISEPOWER - The spectral noise power measured in the Fourier domain outside of the 
;                    Nyquist frequency of the theoretically perfect circular aperture of 
;                    diameter DIAM at wavelength LAMBDA.
;
; OPTIONAL INPUT KEYWORD PARAMETERS:
;	DIAMETER   - The diameter of the telescope aperture. Assumed circular; no allowance 
;                    for spider arms. Units of meters.
;       LAMDA      - The wavelength of the filtergram. Units of meters.
;
;       PLATE_SCALE- The plate scale of the filtergram in arcseconds per pixel.
;
;       Default values for all of the above are supplied as SOT/BFI G-band values.
;
; OPTIONAL OUTPUT KEYWORD PARAMETERS:
;	FILTIM     - Set this keyword to a defined image variable to have a Fourier filtered
;                    image returned. The image is filtered by a sharp-edged circular mask at 
;                    the Nyquist frequency of the theoretically perfect circular aperture of 
;                    diameter DIAM at wavelength LAMBDA.
;
; COMMON BLOCKS: none.
;
; RESTRICTIONS:
;
;
; PROCEDURE:
;        The routine returns the spectral noise power, the square root of the average 
;        power spectral density outside the theoretical telescope Nyquist limit. 
;
; MODIFICATION HISTORY:
;v0.1    Written by T. Berger, 2-Feb-07, LMSAL
;-
;-----------------------------------------------------------------------------
;

;Process keywords:
if KEYWORD_SET(diameter) then diam = diameter else diam = 0.50  ;default SOT 50-cm diameter telescope
if KEYWORD_SET(lambda) then lam = lambda else lam = 430.5e-09   ;default G-band
if KEYWORD_SET(plate_scale) then ps = plate_scale else ps = 0.054  ;default SOT BFI, arcseconds/pixel
if KEYWORD_SET(quiet) then loud = 0 else loud = 1
if KEYWORD_SET(filtim) then fil = 1 else fil = 0

;Image size
sz = SIZE(im)
xn = sz[1]
yn = sz[2]

;Apodize image. Insert in square array if needed for symmetric
;frequency sampling. Use Hamming function since it doesn't go to 0 at the edges.
if xn ne yn then begin
    if loud then PRINT,'Non-square image: embed image in square array to get equal frequency sampling...'
    if xn gt yn then begin
        nn = xn
        imn = FLTARR(xn,xn)
        imn[0,(nn-yn)/2]= HANNING(xn,yn,alpha=0.54)*im 
    end else begin
        nn = yn
        imn = FLTARR(yn,yn)
        imn[(nn-xn)/2,0]= HANNING(xn,yn,alpha=0.54)*im
    end
end else begin
    nn = xn
    imn = HANNING(nn,nn,alpha=0.54)*im
end
if loud then PRINT,'Image size: ',STRTRIM(nn,2),'x',STRTRIM(nn,2),' pixels'

;Spatial and frequency sampling:
deltax = ps*725          ;km
deltak = 1./(nn*deltax)  ;km-1
if loud then PRINT,'Spatial sampling (pixel scale) = ',STRTRIM(deltax,2),' km'
if loud then PRINT,'Frequency sampling in Fourier space = ',STRTRIM(deltak,2),' km-1'

;Nyquist frequency of telescope/optics
nyqt = diam/lam			        ;Nyquist frequency of Telescope, radians^{-1}
nyqk = nyqt*!DTOR/3600/725              ;1/deltax km^{-1} on Sun
nyqk = nyqk/2                           ;Nyquist = 1/(2*deltax)
nyqk = nyqk/deltak                      ;pixels in Fourier transform image
if loud then PRINT,'Nyquist frequency of the imaging system = ',STRTRIM(nyqk*deltak,2),' km-1'
if loud then PRINT,'Nyquist frequency of the imaging system = ',STRTRIM(nyqk,2),' pixels in FFT image'

;Transform and compute power spectrum:
f = FFT(imn,-1)
p = ABS(f)^2.
p = p[0:nn/2-1,0:nn/2-1]     ;extract unique quadrant
if loud then begin
    WINDOW,xs=nn/2,ys=nn/2
    TVSCL,ALOG(p)
    TVCIRCLE,nyqk,0,0
    XYOUTS,nyqk*COS(45*!DTOR)+2,nyqk*SIN(45*!DTOR)+2,'Nyquist limit of optical system',/DEVICE,chars=2
end

;Mask out the area inside the optical Nyquist circle and measure the
;noise power:
mask = 1 - aperture(2*nyqk,2*nyqk,nn/2,nn/2)
xcoord = INDGEN(nn/2)#REPLICATE(1,nn/2)
ycoord = TRANSPOSE(xcoord)
mask = mask AND (xcoord GT nn/12) AND (ycoord GT nn/12)
noisepower = TOTAL(mask*p)/TOTAL(mask)   ;average level in masked area.
noisepower = SQRT(noisepower)            ;IDL divides by NN in the FFT so no need here.
if loud then PRINT,'Noise power = ',STRTRIM(noisepower,2)

;Return filtered image if requested:
if fil then begin   ;return noise filtered image
    nq2 = 2*nyqk
    mask = aperture(nq2,nq2,nn,nn)
    mask = mask OR aperture(nq2,nq2,nn,nn,xc=0,yc=nn-1)
    mask = mask OR aperture(nq2,nq2,nn,nn,xc=nn-1,yc=nn-1)
    mask = mask OR aperture(nq2,nq2,nn,nn,xc=nn-1,yc=0) 
    filtim = ABS(FFT(mask*f,1)) 
    if xn gt yn then filtim = filtim[*,(nn-yn)/2:(nn+yn)/2-1]
    if yn gt xn then filtim = filtim[(nn-xn)/2:(nn+xn)/2-1,*]
    filtim /= HANNING(xn,yn,alpha=0.54)
end


RETURN,noisepower
END