function  sxt_vignette, index, energy, qtest=qtest
;+
; NAME:
;    sxt_vignette
; PURPOSE:
;    Return the SXT X-ray vignette function determined from WSMR1 data.
;
;    Called by sxt_off_axis to correct image shape (portion) given in index.
; CALLING SEQUENCE:
; 	vignette = sxt_vignette(index)		; Index must be a scalar structure.
; 	vignette = sxt_vignette(array)
;	vignette = sxt_vignette(index,energy)
;
; INPUT:
;  INDEX  - Standard SXT index structure.  Specifies portion of SXT FOV to return.
;	    Determine low or high energy from filter (if ENERGY argument not present).
;  ARRAY  - Vector of length 5 where:
;		  array(0) = x CCD coordinate & array(1) = y CCD coordinate
;		  array(2) = # of x pixels    & array(3) = # of y pixels
;		  array(4) = Resolution as: 0=full, 1=half, 2=quarter
; OPTIONAL INPUT:
;  ENERGY - If = 0 or not present, determine from index (filtb=4 or 5, energy=2)
;	    If = 1, use low-energy result  (overrides index)
;	    If = 2, use high-energy result (overrides index)
; OUTPUT:
;   2-d matrix describing the SXT vignette function normalized to 1.0 at the
;   position of highest effective area (gt_sxt_axis).
; PROCEDURE:
;   The SXT vignette function consists of two non-concentric cones whose slopes
;   and origins were determined from WSMR1 data and refined by analyis of flight
;   data.  See Calibration Notes 36 and 37.  The origins are defined as:
;                                   x     y
;   Inner-region (r < 21 arc-min)  515   633  (determined fr WSMR data)
;   Outer-region (r > 21 arc-min)  533   599  (determined fr WSMR data)
;
;   Note that the out-region centroid location and the cutoff radius (=21) may
;   be changed in the future based on more extensive analysis of flight data.
;
;   The slopes were derived from off-axis measurements at 
;       Low energy		Al-K (8.3 A) and C-K (44 A)
;       High energy		Ag L (4.16 A)
;
;   Be patient.  On a DECstations 5000/200 this routine takes about 
;       	 8 sec for a  512x 512 image
;		53 sec for a 1024x1024 image
; RESTRICTIONS:
;   If index is passed in, must be a scalar
; MODIFICATION HISTORY:
;   12-Jan-94, J. R. Lemen LPARL, Written
;    3-Feb-94, JRL, Check that index is a scalar structure
;    4-Feb-94, JRL, Changed the outer-region cone centroid.
;    4-Oct-95, JRL, Fixed a bug that effected half-res and full-res PFIs
;   16-Jun-08, LWA, Used round in cmask(0,round(j-y0)) = peak < wing. 
;   16-Oct-08, LWA, Revised 16-Jun-98 fix for general case.
;-

if N_elements(index) eq 0 then begin	; No input -- show the doc header
   doc_library,'sxt_vignette'
   return,-1
endif

; Determine the type of the input variable

sz = size(index)
typ = sz(1+sz(0))

; NOTE:  For half & quarter resolution, axis coordinates are divided by 2 
;        and 4, respectively.  This is slightly incorrect.

wingaxis = [533.,599.]	; Axis of outer-region (derived from WSMR data)
			; Axis of inner-region derived from gt_sxt_axis

if (typ EQ 8) then begin		; Yes - then index = SXT index structure 
     if n_elements(index) gt 1 then begin
	tbeep
	message,'Warning: index is NOT a scalar',/cont
	help,index
	message,'Will return vignette function for index(0)',/cont
     endif
     mode = gt_res(index(0))		; Get resolution from "index"
     maxaxis = gt_sxt_axis(index(0))	; Get the maximum effective axis
     maxaxis = maxaxis / (2^mode)	;*** Kluge for gt_sxt_axis
     wingaxis= wingaxis/ (2^mode)
     filtb = gt_filtb(index(0))		; Get the filter
     if (filtb eq 4) or (filtb eq 5) then ienergy = 2 else ienergy = 1
     x0 = gt_corner(index(0),/x)
     y0 = gt_corner(index(0),/y)
     dx = gt_shape(index(0),/x)		; Number of X pixels
     dy = gt_shape(index(0),/y)		; Number of Y pixels
endif else begin
     if n_elements(index) eq 5 then mode = index(4) else begin
	message,'INDEX must be an SXT index or a vector of length 5',/cont
	return, -1
     endelse
     maxaxis = gt_sxt_axis(0)		; Assume full resolution
     maxaxis = maxaxis / (2^mode)	;*** Kluge for gt_sxt_axis
     wingaxis= wingaxis/ (2^mode)
     ienergy = 1			; Default to low-energy case
     x0 = index(0) & y0 = index(1)
     dx = index(2) & dy = index(3)
; The next statement allows us to produce plots for display showing the
; the knee at 21 arc-min
     if keyword_set(qtest) then wingaxis = maxaxis
endelse

x0 = x0 / (2^mode)			;*** Another Kluge/approx (JRL: 4-Oct-95)
y0 = y0 / (2^mode)			;*** Another Kluge/approx (JRL: 4-Oct-95)

; Which energy should we use?

if n_elements(energy) ne 0 then eenergy = energy else eenergy = 0
if (eenergy eq 0) then eenergy = ienergy

; Establish the slope and intercept parameters

Rcut = 21.0		; 21 arc-min (where the two cones intersect)
if eenergy eq 1 then begin		; Low energy case
  slopein  = -0.01244	; relative intensity / arc-min
  slopeout = -0.04274	; relative intensity / arc-min

endif else begin			; High-energy case
  slopein  = -0.02335	; relative intensity / arc-min
  slopeout = -0.03328	; relative intensity / arc-min
endelse

; Set up the intercepts.  Solve for outer-region cone at x = Rcut
interceptin = 1.0	; Y intercept for inner-region
interceptout = interceptin + slopein * Rcut - slopeout * Rcut

; Convert slopes to per pixel [Full, half or quarter]:
slopein  = slopein / 60. * gt_pix_size(mode)	; rel int/pixel
slopeout = slopeout / 60. * gt_pix_size(mode)	; rel int/pixel
Rcut = Rcut * 60. / gt_pix_size(mode)		; = 513.7 pixels for full res

; - - - - - This table gives the values for Full-Resolution - - - - - -
;		 per pixel	per arc-min
;  Interceptin=	 1.0				Both Energies
;
;  slopein    =	-0.0005085 	0.01244 	Low Energy
;  slopeout   =	-0.0017472  	0.04274 
;  Interceptout= 1.63630			Rcut = 21 arc-min
;
;  slopein    =	-0.0009545	0.02335		High Energy
;  slopeout   =	-0.0013605	0.03328
;  Interceptout= 1.20853			Rcut = 21 arc-min
; - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

; Setup output variable

cmask = fltarr(dx,dy) 
xp = maxaxis(0)  & yp = maxaxis(1)	; Axis of inner-region
xw = wingaxis(0) & yw = wingaxis(1)	; Axis of outer-region

;--------------------------------------------------------

; Do the Hard work - Use a 1-d loop 
; (This is faster than vectorizing, because uses less memory)

xx = x0 + findgen(dx)			; X vector
yo=fix(y0)
for j=yo,yo+dy-1 do begin		; Loop on Y
  rp = sqrt((xx-xp)^2 + (float(j)-yp)^2)
  peak = interceptin + slopein * rp
  rw = sqrt((xx-xw)^2 + (float(j)-yw)^2)
  wing = interceptout + slopeout * rw
  cmask(0,j-yo) = peak < wing		; Take the smaller of the two
endfor
;--------------------------------------------------------

return, cmask
end