FUNCTION sxt_psd, index, xc, yc, $
	dim=dim, div=div, fourier=fourier, moffat=moffat, log10t=log10t

;+
;NAME:
;	sxt_psd
;
;PURPOSE:
;	Returns an image of the SXT point-spread-function.
;
;INPUT PARAMETERS:
;	index = SXT image index, used to get image center, filter.
;				OR
;       filter = number of the X-ray analysis filter, GT-FILTB convention.
;	xc, yc = center of FOV in CCD Full Resolution pixels.
;
;OUTPUT:
;       psf = image of the psf, full-resolution pixels.  Default is:
;               odd dimension (65x65), array centered,
;               normalized to unit integral.  This is
;               intended for use by IDL procedure CONVOL.
;
; OPTIONAL INPUT KEYWORDS:
;	dim = size of output PSF image in Full-resolution pixels.  Default = 65.
;
;       div = No of divisions per pixel	;JS
;
;	fourier = Return psf ready for use by IDL procedure FFT.
;		Psf image is dimensioned power of 2 (64 x 64),
;		corner centered, normalized to unit mean.
;
;	moffat = returns Moffat psf only.  Scattering wings not
;		included.
;
;	log10t = Log10 of temperature.  The Moffat parameters
;		are wavelength, and therefore, temperature
;		dependent.  Default uses values derived
;		assuming equal emission measures at all
;		temperatures  5.5 <= log10t <= 8.0.
;
;
; CALLING EXAMPLES:
;	psf = SXT_PSF(index)
;	psf = SXT_PSF(filter, xc, yc)
;
; METHOD:
;	Uses the Moffat function fits described by Piet Martins.
;	Adds scattering wings described by Acton.
;	Wavelength weighting done by SXT_WEIGHT procedure.
;
; RESTRICTION:
;	Scattering is only known for thin filters.  Scattering
;	for the thick filters is unknown but is expected to be larger.
;
; HISTORY:
;	Written April 4, 1994   Barry LaBonte
;	Added /DIM keyword  May 23, 1994  BJL
;	Added /FOURIER, /MOFFAT keywords, scattering wings.  June 27, 1994 BJL
;	Modified to use Moffat parameters weighted
;		by SXT passband - Mewe spectrum product,
;		added /LOG10T keyword.  March 24, 1995  BJL
;		added /DIV  keyword and made it double precision JS Jan. 98
;
;-


; MOFFAT Least-squares Fit parameters for equal emission measure
; Filter order is given by GT_FILTB():
;  1   Open       Op
;  2   Al.1       A1
;  3   AlMg       AM
;  4   Be119      Be
;  5   Al12       A2
;  6   Mg3        Mg
aeqem = double([0.77, 0.76, 0.76, 0.96, 0.71, 0.75])
beqem = double([1.53, 1.52, 1.52, 1.66, 1.48, 1.51])
xceqem = double([529.0, 528.8, 528.7, 533.0, 527.8, 528.7])
yceqem = double([585.7, 584.9, 583.7, 586.3, 579.6, 583.8])
eps0 = double([-5.54e-4, -3.18e-3, -5.05e-3, 8.00e-2, -2.47e-2, -6.82e-3])
eps1 = double([1.40, 1.41, 1.42, 1.12, 1.49, 1.42])
eps2 = double([-3.66, -3.69, -3.71, -3.32, -3.86, -3.72])
eps3 = double([2.69, 2.70, 2.73, 2.72, 2.81, 2.72])


; Temperature correction coefficients
atc = double([ [ 0.223,  0.230,  0.123, -0.208,  0.174,  0.072], $
        [-0.447, -0.475, -0.235,  0.118, -0.465, -0.162], $
        [ 0.180,  0.194,  0.096,  0.005,  0.208,  0.074] ])
btc = double([ [ 0.218,  0.224,  0.124, -0.155,  0.170,  0.074], $
        [-0.426, -0.453, -0.225,  0.084, -0.443, -0.155], $
        [ 0.168,  0.182,  0.088,  0.007,  0.196,  0.066] ])
xctc = double([ [  5.39,   5.53,  3.01, -4.41,   4.19,  1.79], $
         [-10.65, -11.32, -5.61,  2.45, -11.08, -3.86], $
         [  4.24,   4.59,  2.25,  0.15,   4.93,  1.71] ])
yctc = double([[ 31.93,  32.78,  19.59, -5.25,  24.88,  12.43], $
        [-59.11, -62.77, -31.45,  0.60, -60.86, -21.45], $
        [ 22.06,  23.93,  10.81,  1.50,  26.11,   7.68]])
e0tc = double([[0.108,    0.111,   0.060,  -0.088,   0.084,   0.036], $
        [-0.213,   -0.226,   -0.112,   0.049,   -0.222,  -0.077], $
        [0.085,   0.092,   0.045,  0.003,   0.099,   0.034]])
e1tc = double([[-0.371,   -0.381,   -0.207,    0.309,   -0.289,   -0.123], $
        [0.734,    0.781,    0.387,   -0.172,    0.765,    0.266], $
        [-0.293,   -0.317,   -0.156,  -0.010,   -0.341,   -0.118]])
e2tc = double([[0.934,    0.959,    0.548,   -0.452,    0.727,    0.338], $
        [-1.786,    -1.898,   -0.946,    0.216,    -1.849,   -0.648], $
        [0.689,    0.746,    0.352,   0.035,     0.808,    0.259]])
e3tc = double([[-0.667,   -0.685,   -0.413,   0.066,   -0.520,   -0.264], $
        [1.226,     1.302,    0.653,   0.017,     1.261,    0.445], $
        [-0.454,   -0.493,   -0.221,  -0.033,   -0.539,   -0.155]])
etc = [ [[e0tc]], [[e1tc]], [[e2tc]], [[e3tc]] ]

; Scattering wings parameters vs increasing energy (Not implemented.
;	 uses thin parameters for all filters.)
slevel = 0.00168d0
spower = -1.89d0


; Handle inputs
sz = SIZE(index)
type = sz(N_ELEMENTS(sz)-2)
IF( type eq 8 ) THEN BEGIN
	filter = GT_FILTB(index(0))
	xy = GT_CENTER(index(0))
	xc = xy(0)
	yc = xy(1)
   ENDIF ELSE BEGIN
	filter = index
ENDELSE
; Which fit to use?
idex = filter - 1

; Get parameter values
afit = aeqem(idex)
bfit = beqem(idex)
xcent = xceqem(idex)
ycent = yceqem(idex)
eps = [eps0(idex), eps1(idex), eps2(idex), eps3(idex)]

; Corrections to equal emission measure values as quadratic
; in Log10(temperature)
IF( KEYWORD_SET(log10t) ) THEN BEGIN
	dt = log10t - 6.
	dt2 = dt^2
	afit = afit + atc(idex,0) + atc(idex,1)*dt + atc(idex,2)*dt2
	bfit = bfit + btc(idex,0) + btc(idex,1)*dt + btc(idex,2)*dt2
	xcent = xcent + xctc(idex,0) + xctc(idex,1)*dt + xctc(idex,2)*dt2
	ycent = ycent + yctc(idex,0) + yctc(idex,1)*dt + yctc(idex,2)*dt2
	eps = eps + etc(idex,0,*) + etc(idex,1,*)*dt + etc(idex,2,*)*dt2
ENDIF

; How big to make output?
orsize = 65
IF( KEYWORD_SET(fourier) NE 0 ) THEN orsize = 64
IF( KEYWORD_SET(dim) NE 0 ) THEN  orsize = dim*div	;JS
if orsize/2*2 eq orsize then orsize=orsize+1

orhalf = orsize/2

; Coordinate images
x = FINDGEN(orsize)/div - FLOAT(orhalf)/div		;JS
x = REBIN(x, orsize, orsize)
y = TRANSPOSE(x)

; Polar angle
phi = (ATAN( yc - ycent, xc - xcent ))
; Distance
d = SQRT( ( xc - xcent )^2 + ( yc - ycent )^2 ) / 500.

; Ellipticity
epsilon = eps(0) + eps(1) * d + eps(2) * d^2	$
	 + eps(3) * d^3

; "Radius"
r2 = x^2 * (1. + epsilon *(1. + COS(2.*phi)))		$
	+ y^2 * (1. + epsilon*(1. - COS(2.*phi)))	$
	+ x * y * 2. * epsilon * SIN(2.*phi)

; Point-spread-function
psf = 1. / (1. + r2/afit^2 )^bfit

; Add scattering wings if needed
IF( KEYWORD_SET(moffat) EQ 0 ) THEN BEGIN
	r2 = x^2 + y^2
	r2(orhalf, orhalf) = 1.
	scat = slevel * 0.5d0^spower * r2 ^(spower/2)
	wing = WHERE(r2 GT 10.)
	if wing(0) ne -1 then psf(wing) = psf(wing) > scat(wing)
ENDIF

; What style of output?
IF( KEYWORD_SET(fourier) EQ 0 ) THEN BEGIN
;	Normalize to unit total
	psf = psf / TOTAL(psf)
  ENDIF ELSE BEGIN
;	Normalize to unit average
	pjunk = STDEV(psf, average)
	psf = psf / average
;	Shift to corner centered
	psf = SHIFT(psf, -orhalf, -orhalf)
ENDELSE

RETURN, psf
END