;+
;
; Make a scattered light image from an SXT composite image.
;

PRO sxt_scatter, composite, scatter, level=level, power=power

;
; INPUT PARAMETERS:
;	composite = composite image, Full frame.
;
; OUTPUT PARAMETERS:
;	scatter = scattered light image.
;
; OPTIONAL INPUT KEYWORDS:
;	level = level of scattered light.  Default = 0.00168.
;	power = power-law index of scattered light.  Default = -1.89
;
; RESTRICTIONS:
;	The scattering function is only valid for the thin filters:
;	Al.1, AlMg, Mg3.
;	Uses Fourier transforms for convolution.
;	USES BIG ARRAYS, TAKES TIME, USES MEMORY.
;
; HISTORY:
;	Written June 17, 1994   Barry LaBonte
;	Cleaned up the code   June 21, 1994  BJL
;	Default values are Acton's result for 6-Sep-92 image  June 23,1994 BJL
;
;-

; How big is the input?
sz = SIZE(composite)
s1 = sz(1)
s2 = sz(2)
IF( s1 NE s2 ) THEN  $
	PRINT, 'WARNING: Image does not include whole Sun.'
; Scale relative to Half-resolution pixels
pix = 512./FLOAT(s1)

; Handle keywords
IF( KEYWORD_SET(level) NE 0 ) THEN slevel = level ELSE slevel = 0.00168
IF( KEYWORD_SET(power) NE 0 ) THEN spower = power ELSE spower = -1.89

; SCATTER FUNCTION  -----------------------------------------------
; Radius^2 = x^2 + y^2
psf = ( FINDGEN( 2*s1 ) - s1 )^2
psf = REBIN( psf, 2*s1, 2*s2 )
ypsf = ( FINDGEN( 2*s2 ) - s2 )^2
ypsf = REFORM( ypsf, 1, 2*s2 )
ypsf = REBIN( ypsf, 2*s1, 2*s2 )
psf = TEMPORARY(psf) + ypsf
ypsf = 0

; Shift center of function to (0,0)
psf = SHIFT( TEMPORARY(psf), -s1, -s2 )
psf(0,0) = 1.

; Get the correct function
psf = slevel * pix^(2+spower) * TEMPORARY(psf)^(spower/2)
psf(0,0) = 1.

; Find the correct normalization
pnorm = REBIN( psf, 1, 1 )

; "Subtract" the original image
psf(0,0) = 0.

; Take the Fourier Transform
psf = FFT( TEMPORARY(psf)/pnorm(0), -1 )

; NOW THE IMAGE  ------------------------------------------------
; Insert into double sized array to avoid Fourier wraparound
big = FLTARR( 2*s1, 2*s2 )
big(s1/2:s1-1+s1/2, s2/2:s2-1+s2/2) = composite

; Transform,  convolve, and back transform
big = FFT( ( FFT( TEMPORARY(big), -1 ) * psf ), 1 )

; Extract the valid region
scatter = FLOAT( TEMPORARY( big(s1/2:s1-1+s1/2, s2/2:s2-1+s2/2) ) )

RETURN
END