; Straylight masking procedure mask_straylight_v5p2.pro
; ---------------------------------------------------
;
; mask_straylight_v5p2,bcd_list,dmask_list,tmass_list,CHDMSK=chdmsk,UNAGRR=unaggr
;
; Uses procedures to mask each source of straylight:
;
; straymask[1-4]fpa: scattering from the hole in the FPA cover (especially 
; strong in channels 1 and 2)
; 
; straymask1c: unknown source of scattering in channel 1 off lower right 
; corner of the BCD ("zone 1C").
;
; straymask[1-4]det: scattering off edge of detector in channels 3 and 4
;
; straymask3s: unknown source of scattering producing "conic section"
; pattern in channel 3 off top right-hand edge of BCD.
;
; Thresholds for applying the masks are hard-coded in, but easy to alter
; if necessary by editing this file. Reddening/stellar color information
; from the J-K color is incorperated into the thresholding assuming:
;
;  ch1 = K - 0.298(J-K)
;  ch2-4 = K - 0.387(J-K)
;
;  in line with the extinction law of Indebetouw et al.
;
; Mask options: 
; 
; 1) Default - makes files named "smask_<filename>" which contain the 
; straylight masks. Useful for making "advisory" mask mosaics indicating
; where straylight is likely to be present.
;
; 2) CHDMSK keyword set - adds straylight masks to dmasks (using bit 3)
; the code overwrites the dmasks with this option set so make a copy of them 
; first! To make mosaics with the masking applied, add 8 to the 
; DCE_Status_Mask_Fatal_BitPattern in mopex. Bit 3 is actually reserved for 
; electronic crosstalk, but is unlikely to be used in the near future.
;
; 3) UNAGGR keyword set - uses the small masks even for bright objects where
; the larger area mask would usually be employed. Use for for low coverage 
; datasets where masking of large areas may lead to holes in coverage.
;
; Inputs: list of BCDs, matching list of dmasks, list of 2MASS stars
; from the post-BCD pipeline (SPITZER*_irsa.tbl)
;
; Output: modified (overwritten) dmask file with the straylight
; regions masked with bit 3 (8). This bit is reserved for electronic
; crosstalk, so is unlikely to be used in the near future.
; If scattering is found, the type of scattering and position of the 
; scattering star are written to the BCD header in the history.
;
; ML 6/16/04.
;
; Revisions:
; Bug fixed in straymask3 6/18/04
; Multiple potential scatterers now dealt with correctly,I/O bug fixed 6/28/04
; RGA 7/15/05 mask_straylight_v3:
; swapped multiple IF statments for CASE statements
; eliminated some never-used or one-time-used variables
; altered the position screening so that only sources off the FPA are checked
; altered the logic for contruction of the final dmask to not presume 
; initially unset bit.
; ML 7/28/04 mask_straylight_v4: 
; added RGA's new masks for channel 3 and 4 FPA scattering, 
; channel 3 "conic section" scattering. Tidied up and made more user-friendly
; for public release
; ML 10/29/04 mask_straylight_v5:
; tweaked zone 1C
; included ghostbuster code from RGA for channels 1 and 2
; ML 12/04: put in colour dependent thresholds
; ML 12/23/04: fixed bug so dmask files now get correct bit values
; ML 04/08/05: added UNAGGR keyword so that only the small masks are used 
; if the keyword is set, make both the CHDMASK and UNAGGR keywords simple
; switches -> version 5p1.
; ML 05/05/05: put in Rick's newer ghostbuster code which includes channel
; 3 and 4 ghosts, fixed bug so ch1 threshold is now consistent with comments 
; -> v5p2. 

pro straymask1fpa, x0, y0, mask, mask2, hdr

; INPUT
; (x0,y0) = detector pixel location of star
; OUTPUT
; mask = radial + azimuthal zone of most intense stray light
; mask2 = azimuthal zone of all stray light
; hdr = fits image header

x = lindgen(256,256) mod 256
y = lindgen(256,256) / 256


; If the source is outsid the "1a" and "1b" zones
; then return blank masks. 

if (y0 lt 250) or (y0 gt 313) or (x0 le -15) or (x0 gt 285) or $
  (sqrt((x0-(128+6))^2+(y0-(128-30))^2) lt 208) or $
  (sqrt((x0-(128+26))^2+(y0-(128-30))^2) gt 250) then begin 

    mask = bytarr(256,256)
    mask2 = bytarr(256,256)

endif else begin

; Calulate radial coords wrt the source 

  r = sqrt((x-x0)^2+(y-y0)^2)
  theta = atan(y-y0,x-x0)

; The mean stray light seems to point toward a fixed
; location slightly offset from the center of the array (~64,128).

  r0 = sqrt(float(x0-134)^2+float(y0-128)^2)
  theta0 = atan(y0-64,x0-128) 

; rstray, thetastray are functions of (x0,y0) (or r0,theta0)
; delr, deltheta may be functions of (x0,y0), but might also be fixed values.
; Also needs a cut on a specific range of (x0,y0) and (r0,theta0) 
; There may also be asymmetrical offsets involved.

  rstray= 55 + 1.4*(r0-175)
  delr = 20

; thetastray is theta0 - pi if theta0 gt 0, theta0+pi if theta0 lt 0

  thetastray = theta0-!pi*(theta0 gt 0)+!pi*(theta0 lt 0)

; dtor is the system variable to convert degrees to radians

  deltheta = 17*!dtor

; Define masks as regions specified in polar coords wrt source.
; mask = mask for the brightest stray light components
; mask2 = has no radial cut and flags all potential rays. 

  mask2 = abs(theta-thetastray) lt deltheta
  mask = abs(r-rstray) lt delr and mask2

	xs = string(x0)
	ys = string(y0)

	position = xs+','+ys
	radiuss = string(rstray) 

	print,'Channel 1 FPA cover scatterer found at ',x0,y0

	sxaddpar,hdr,'HISTORY', $
'mask_straylight: potential FPA cover scattering; star at'
	sxaddpar,hdr,'HISTORY',position
	sxaddpar,hdr,'HISTORY','mask_straylight: rstray='
	sxaddpar,hdr,'HISTORY',radiuss

endelse

end


pro straymask2fpa, x0, y0, mask, mask2, hdr

; INPUT
; (x0,y0) = detector pixel location of star
; OUTPUT
; mask = radial + azimuthal zone of most intense stray light
; mask2 = azimuthal zone of all stray light

x = lindgen(256,256) mod 256
y = lindgen(256,256) / 256

; If the source is outside the "2a" and "2b" zones
; then return blank masks.

if (y0 lt 240) or (y0 gt 321) or (x0 le -15) or (x0 gt 286) or $
  (sqrt((x0-(128-4))^2+(y0-(128-34))^2) lt 214) or $
  (sqrt((x0-(128+0))^2+(y0-(128-25))^2) gt 260) then begin
    mask = bytarr(256,256)
    mask2 = bytarr(256,256)
endif else begin

; Calulate radial coords wrt the source
  r = sqrt((x-x0)^2+(y-y0)^2)
  theta = atan(y-y0,x-x0)

; The mean stray light seems to point to a fixed
; location slightly the center of the array (~64,128).

  r0 = sqrt(float(x0-122)^2+float(y0-128)^2)
  theta0 = atan(y0-64,x0-128) 

; rstray, thetastray are functions of (x0,y0) (or r0,theta0)
; delr, deltheta may be functions of (x0,y0), but might also be fixed values.
; The coefficients here are just placeholders.
; Also needs a cut on a specific range of (x0,y0) and (r0,theta0) 
; There may also be asymmetrical offsets involved.

  rstray= 50 + 1.4*(r0-175)
  delr = 28
  thetastray = theta0-!pi*(theta0 gt 0)+!pi*(theta0 lt 0)
  deltheta = 17*!dtor

; Define masks as regions specified in polar coords wrt source.
; mask = mask for the brightest stray light components
; mask2 = has no radial cut and flags all potential rays. 

  mask2 = abs(theta-thetastray) lt deltheta
  mask = abs(r-rstray) lt delr and mask2

	xs = string(x0)
	ys = string(y0)
	position = xs+','+ys
	radiuss = string(rstray) 

	print,'Channel 2 FPA cover scatterer found at ',x0,y0

	sxaddpar,hdr,'HISTORY', $
'mask_straylight: potential FPA cover scattering; star at'
	sxaddpar,hdr,'HISTORY',position

endelse

end

pro straymask3det, x0, y0, mask, mask2, hdr
; INPUT
; (x0,y0) = detector pixel location of star
; OUTPUT
; mask = radial + azimuthal zone of most intense stray light
; mask2 = azimuthal zone of all stray light
; hdr = fits image header
;
; If the source is outside the scattered light zones
; then return blank masks. 
; scattered light zone for channels 3 and 4 is a border about 
; 4 pixels wide around the arrays

if ((x0 lt -17) or (x0 gt 274) or (y0 lt -17) or (y0 gt 274)) or $
((x0 gt -13) and (x0 lt 270) and (y0 gt -13) and (y0 lt 270)) then begin

    mask = bytarr(256,256)
    mask2 = bytarr(256,256)

endif else begin

; pattern in 3 and 4 is just a strip along the row/column affected. 
; pattern is concentrated in a well-defined region about 5 pixels wide
; Bad scattered light can extend across the whole array (=mask2)
; mask = mask for the brightest stray light components
; mask2 = has no radial cut and flags all potential rays. 
; start with scatterers on left-hand edge:


print,'star inside boundary at ',x0,y0

    mask = bytarr(256,256)
    mask2 = bytarr(256,256)

if (x0 le -13) and (y0 ge 0) and (y0 le 255) then begin

	ymin = y0 - 3
	ymax = y0 + 3

	if ymin lt 0 then ymin=0
	if ymax gt 255 then ymax=255

	for j=ymin,ymax do begin

		mask[0:5,j]=1
		mask[28:39,j]=1

	endfor

	ymin2 = y0 - 10
	ymax2 = y0 + 10

	if ymin2 lt 0 then ymin2=0
	if ymax2 gt 255 then ymax2=255

	for j=ymin2,ymax2 do begin

		mask2[*,j]=1

	endfor
	

endif

if (x0 ge 269) and (y0 ge 0) and (y0 le 255) then begin

	ymin = y0 - 3
	ymax = y0 + 3

	if ymin lt 0 then ymin=0
	if ymax gt 255 then ymax=255

	for j=ymin,ymax do begin

		mask[250:255,j]=1
		mask[231:237,j]=1
		mask[205:225,j]=1

	endfor


	ymin2 = y0 - 10
	ymax2 = y0 + 10

	if ymin2 lt 0 then ymin2=0
	if ymax2 gt 255 then ymax2=255

	for j=ymin2,ymax2 do begin

		mask2[*,j]=1

	endfor
	
endif


if (y0 le -13) and (x0 ge 0) and (x0 le 255) then begin

	xmin = x0 - 3
	xmax = x0 + 3

	if xmin lt 0 then xmin=0
	if xmax gt 255 then xmax=255

	for j=xmin,xmax do begin

		mask[j,0:5]=1
		mask[j,28:40]=1

	endfor

	xmin2 = x0 - 10
	xmax2 = x0 + 10

	if xmin2 lt 0 then xmin2=0
	if xmax2 gt 255 then xmax2=255

	for j=xmin2,xmax2 do begin

		mask2[j,*]=1

	endfor
	
endif

if (y0 ge 270) and (x0 ge 0) and (x0 le 255) then begin

;	print,'high-y scatterer found at ',x0,y0

	xmin = x0 - 3
	xmax = x0 + 3

	if xmin lt 0 then xmin=0
	if xmax gt 255 then xmax=255

	for j=xmin,xmax do begin

		mask[j,250:255]=1
		mask[j,231:241]=1
		mask[j,213:230]=1

	endfor


	xmin2 = x0 - 10
	xmax2 = x0 + 10

	if xmin2 lt 0 then xmin2=0
	if xmax2 gt 255 then xmax2=255

	for j=xmin2,xmax2 do begin

		mask2[j,*]=1

	endfor
endif

	xs = string(x0)
	ys = string(y0)
	position = xs+','+ys

	print,'Channel 3 detector edge scatterer at',x0,y0 

	sxaddpar,hdr,'HISTORY','mask_straylight: detector scattering; star at'
	sxaddpar,hdr,'HISTORY',position
endelse
end

pro straymask4det, x0, y0, mask, mask2, hdr
; INPUT
; (x0,y0) = detector pixel location of star
; OUTPUT
; mask = radial + azimuthal zone of most intense stray light
; mask2 = azimuthal zone of all stray light
; hdr = fits image header
;
; If the source is outside the scattered light zones
; then return blank masks. 
; scattered light zone for channels 3 and 4 is a border about 
; 4 pixels wide around the arrays

if ((x0 lt -14.5) or (x0 gt 275) or (y0 lt -14.5) or (y0 gt 276)) or $
((x0 gt -12) and (x0 lt 269) and (y0 gt -12) and (y0 lt 267)) then begin
    mask = bytarr(256,256)
    mask2 = bytarr(256,256)
endif else begin

;print,'scatterer at',x0,y0

; pattern in 3 and 4 is just a strip along the row/column affected. 
; pattern is concentrated in a well-defined region about 5 pixels wide
; Bad scattered light can extend across the whole array (=mask2)

; Define masks as regions specified in polar coords wrt source.
; mask = mask for the brightest stray light components
; mask2 = has no radial cut and flags all potential rays. 
; start with scatterers on left-hand edge:

    mask = bytarr(256,256)
    mask2 = bytarr(256,256)

if (x0 le -12) and (y0 ge 0) and (y0 le 255) then begin

	ymin = y0 - 3
	ymax = y0 + 3

	if ymin lt 0 then ymin=0
	if ymax gt 255 then ymax=255

	for j=ymin,ymax do begin

		mask[0:5,j]=1
		mask[28:39,j]=1

	endfor

	ymin2 = y0 - 10
	ymax2 = y0 + 10

	if ymin2 lt 0 then ymin2=0
	if ymax2 gt 255 then ymax2=255

	for j=ymin2,ymax2 do begin

		mask2[*,j]=1

	endfor
	
endif

if (x0 ge 269) and (y0 ge 0) and (y0 le 255) then begin

	ymin = y0 - 3
	ymax = y0 + 3

	if ymin lt 0 then ymin=0
	if ymax gt 255 then ymax=255

	for j=ymin,ymax do begin

		mask[250:255,j]=1
		mask[212:226,j]=1

	endfor


	ymin2 = y0 - 10
	ymax2 = y0 + 10

	if ymin2 lt 0 then ymin2=0
	if ymax2 gt 255 then ymax2=255

	for j=ymin2,ymax2 do begin

		mask2[*,j]=1

	endfor
	
endif


if (y0 le -12) and (x0 ge 0) and (x0 le 255) then begin

;    print,'low-y scatterer at ',x0,y0

	xmin = x0 - 3
	xmax = x0 + 3

	if xmin lt 0 then xmin=0
	if xmax gt 255 then xmax=255

	for j=xmin,xmax do begin

		mask[j,0:17]=1
		mask[j,33:46]=1

	endfor

	xmin2 = x0 - 10
	xmax2 = x0 + 10

        if xmin2 lt 0 then xmin2=0
	if xmax2 gt 255 then xmax2=255

	for j=xmin2,xmax2 do begin

		mask2[j,*]=1

	endfor
	
endif

if (y0 ge 267) and (x0 ge 0) and (x0 le 255) then begin

;    print,'high-y scatterer at ',x0,y0

	xmin = x0 - 3
	xmax = x0 + 3

	if xmin lt 0 then xmin=0
	if xmax gt 255 then xmax=255

	for j=xmin,xmax do begin

		mask[j,206:255]=1
;		mask[j,212:226]=1

	endfor


	xmin2 = x0 - 10
	xmax2 = x0 + 10

	if xmin2 lt 0 then xmin2=0
	if xmax2 gt 255 then xmax2=255

	for j=xmin2,xmax2 do begin

		mask2[j,*]=1

	endfor


endif

	xs = string(x0)
	ys = string(y0)
	position = xs+','+ys

	print,'Channel 4 detector edge scatterer at',x0,y0 
	sxaddpar,hdr,'HISTORY','mask_straylight: detector scattering; star at'
	sxaddpar,hdr,'HISTORY',position

endelse

end


pro straymask1c, x0, y0, mask, mask2,hdr
; 
; Stray light masking for Zone 1C.
; In contrast to the masks for the 1A+1B zones, this mask is always set at 
; the same location. MASK2 is set whenever the source is within the zone, but
; MASK is only set when the source is in the lower portion of the zone where
; it tends to produce stronger stray light. 
; 
; INPUT
; (x0,y0) = detector pixel location of star
; OUTPUT
; mask = mask for stronger portion stray light zone 1c
; mask2 = larger mask for all stray light zone 1c

x = lindgen(256,256) mod 256
y = lindgen(256,256) / 256

; If the source is outside the "1c" zone
; then return blank masks. 

;if (y0 lt -64) or (x0 le 210) or (x0 gt 279+0.3*(y0+40)) or $
;  (y0 gt 0.4*x0-120) then begin 

if (y0 lt -64) or (x0 le 210) or (x0 gt 271+0.3*(y0+40)) or $
  (y0 gt 0.4*x0-120) then begin 
    mask = bytarr(256,256)
    mask2 = bytarr(256,256)
endif else begin

  mask2 = y+x*0.68 gt 43 and y+x*0.68 lt 150
  mask = mask2*(y0 lt 0.4*x0-135)

	xs = string(x0)
	ys = string(y0)

	position = xs+','+ys
	sxaddpar,hdr,'HISTORY','mask_straylight: potential zone 1C scattering; star at'
	sxaddpar,hdr,'HISTORY',position

endelse

end



pro straymask3fpa, x0, y0, mask, mask2,hdr
;
; Defines the masks for reflections from the FPA covers for Ch. 3.
; These are relatively weak compared to the chip-edge form of straylight.
;
; INPUT
; (x0,y0) = detector pixel location of star
; OUTPUT
; mask = radial + azimuthal zone of most intense stray light
; mask2 = azimuthal zone of all stray light

x = lindgen(256,256) mod 256
y = lindgen(256,256) / 256

; If the source is outside the "3B" zone AND outside the 
; "3C" zone AND outside the "3D" zone 
; then return blank masks. 

if ((y0 lt (256-x0)+310) or (x0 gt 288) or (y0 gt 320)) $
  and ((y0 gt (x0-256)*0.8-30) or (x0 gt 285) or (y0 lt -64)) $
  and ((y0 gt -x0-55) or (x0 lt -10) or (y0 lt -64)) $

  then begin

    mask = bytarr(256,256)
    mask2 = bytarr(256,256)

endif else begin



; Calulate radial coords wrt the source 

  r = sqrt((x-x0)^2+(y-y0)^2)
  theta = atan(y-y0,x-x0)

; The mean stray light seems to point toward a fixed
; location slightly offset from the center of the array.

  inzone3c = x0 gt 128 and y0 lt 128
  r0 = sqrt(float(x0-128)^2+float(y0-145)^2)
  theta0 = atan(y0-145,x0-128) 

  if (inzone3c) then begin

    r0 = sqrt(float(x0-128)^2+float(y0-120)^2)
    theta0 = atan(y0-120,x0-128) 

  endif



; rstray, thetastray are functions of (x0,y0) (or r0,theta0)
; delr, deltheta may be functions of (x0,y0), but might also be fixed values.
; The coefficients here are just placeholders.
; Also needs a cut on a specific range of (x0,y0) and (r0,theta0) 
; There may also be asymmetrical offsets involved.

  rstray= 25 + 1.4*(r0-175) + 35*inzone3c
  delr = 30 + 10*inzone3c
  thetastray = theta0-!pi*(theta0 gt 0)+!pi*(theta0 lt 0)
  deltheta = !dtor*10

; Define masks as regions specified in polar coords wrt source.
; mask = mask for the brightest stray light components
; mask2 = has no radial cut and flags all potential rays. 

  mask2 = abs(theta-thetastray) lt deltheta
  mask = abs(r-rstray) lt delr and mask2

	xs = string(x0)
	ys = string(y0)

	position = xs+','+ys

;	print,'Channel 3 FPA scatterer found at ',x0,y0

	sxaddpar,hdr,'HISTORY',$
'mask_straylight: potential FPA cover scattering; star at'
	sxaddpar,hdr,'HISTORY',position

endelse

end



pro straymask3s, x0, y0, mask, mask2,hdr
;
; Masking for the conic section splotches that are rarely seen in Ch. 3.
; The campaign VW test sees this stray light >50% of the time when 
; 283<x0<285 and -17<y0<273, but it's not clear how to predict in better 
; detail whether or not the stray light will be present. 
; The GFLS (4958976) has an example on the other side of the array, x=-28 to
; -35 at y=194-195. 
; Other columns in the test at x0~278.5 and x0~287 did not show stray light,
; but there were no measurments in the ranges 279-283 and 285-296. 
; The pattern and brightness of the straylight within the parabolic mask is
; highly variable.
;
; INPUT
; (x0,y0) = detector pixel location of star
; OUTPUT
; mask = parabolic zone of possible stray light.
; mask2 = (Not calculated)
;
; If the source is outside the possible scattering zones return blank
; masks:
; 

if ((abs(x0-284) gt 1 or abs(y0-128) gt 145) and (abs(x0+32) ge 4 or $
((y0 le 190) or (y0 ge 200)))) then begin

    mask = bytarr(256,256)
    mask2 = bytarr(256,256)

endif else begin



; make mask according to which side of the array the scatter is on

  if x0 gt 0 then begin

   x = lindgen(256,256) mod 256 - 84.
   y = lindgen(256,256) / 256 
   R = -70.

   mask =  (y-y0)^2 lt R*(x-x0)

  endif else begin

   x = lindgen(256,256) mod 256 - 84
   y = lindgen(256,256) / 256
   R = 70.

   mask =  (y-y0)^2 lt R*(x-x0)

   endelse

	xs = string(x0)
	ys = string(y0)

	position = xs+','+ys

;	print,'Channel 3 conic section scatterer found at ',x0,y0
	sxaddpar,hdr,'HISTORY',$
'mask_straylight: potential conic section scattering; star at'
	sxaddpar,hdr,'HISTORY',position

endelse

end


pro straymask4fpa, x0, y0, mask, mask2,hdr
;
; This routine masks the FPA cover stray light in Ch. 4.
; The stray light is only strong from zone 4D
;  
;
; INPUT
; (x0,y0) = detector pixel location of star
; OUTPUT
; mask = radial + azimuthal zone of most intense stray light
; mask2 = azimuthal zone of all stray light

x = lindgen(256,256) mod 256
y = lindgen(256,256) / 256

; If the source is outside the "4d" zone AND outside the 
; "4c" zone AND outside the "4a" zone then return blank masks. 

if ((y0 lt -64) or (x0 le -24+0.3*(y0+48)) or (x0 gt 10) or (y0 gt -30-x0*0.6)) $
  and (abs(x0+42) gt 5 or abs(y0+19) gt 8) $  ; extra part of "4D"
  and (sqrt((x0-274)^2+(y0+40)^2) gt 6.5)  $  ; "4C"
  and (abs(x0-215) gt 45 or abs(y0+32) gt 12)  $  ; weak part of "4C"
  and (sqrt((x0+17)^2+(y0-308)^2) gt 5) then begin  ; "4A"

    mask = bytarr(256,256)
    mask2 = bytarr(256,256)

endif else begin

; Calulate radial coords wrt the source 

  r = sqrt((x-x0)^2+(y-y0)^2)
  theta = atan(y-y0,x-x0)

; The mean stray light seems to point toward a fixed
; location slightly offset from the center of the array.

  r0 = sqrt(float(x0-128)^2+float(y0-145)^2)
  theta0 = atan(y0-145,x0-128) 

; rstray, thetastray are functions of (x0,y0) (or r0,theta0)
; delr, deltheta may be functions of (x0,y0), but might also be fixed values.
; The coefficients here are just placeholders.
; Also needs a cut on a specific range of (x0,y0) and (r0,theta0) 
; There may also be asymmetrical offsets involved.

  rstray= 45 + 1.4*(r0-175) + 40*(x0 gt 128) -30*(x0 gt 260) - 35*(y0 gt 260)
  delr = 20
  thetastray = theta0-!pi*(theta0 gt 0)+!pi*(theta0 lt 0)
  deltheta = 10*!dtor

; Define masks as regions specified in polar coords wrt source.
; mask = mask for the brightest stray light components
; mask2 = has no radial cut and flags all potential rays. 

  mask2 = abs(theta-thetastray) lt deltheta
  mask = abs(r-rstray) lt delr and mask2

	xs = string(x0)
	ys = string(y0)
	position = xs+','+ys

;	print,'Channel 4 FPA scatterer found at ',x0,y0
	sxaddpar,hdr,'HISTORY',$
'mask_straylight: potential FPA cover scattering; star at'
	sxaddpar,hdr,'HISTORY',position

endelse

end


pro optics10_ghost, nsize, lambda, D, f_ratio, obstfrac, spiderfrac, $
  um_pix=um_pix, a_zcoeff=a_zcoeff, beam0=beam0

; INPUT
; nsize       = size of grid for mapping aperture and beams (pixels)
; lambda      = wavelength of light (meters)
; D           = diameter of aperture (meters) only used to translated scales 
;               from pixels to m, or pixels to arcsec (see f_ratio)
; f_ratio     = f-ratio of the input beam (focal length / D), which sets the
;               basic geometry for the calculations. 
; obstfrac    = fractional diameter of circular secondary obstruction
; spiderfrac  = defines fractional width (and optionally positions) of 
;               secondary support structure
;               spiderfrac[0] = >0 = fractional linear width of spider vanes 
;                               <0 = angular width (rad) of spider vanes
;               spiderfrac[1:*] = if present these define the number and 
;               position angles (radians) of the spider arms. If omitted, 
;               the default is four arms at [0,1,2,3]*!pi/2.
; um_pix      = um/pix scale for beam and shutter images
; a_zcoeff    = Zernike coefficients for the aperture
;               Each of these sets of coefficients is applied as a phase 
;               modification at the relevant location. The arrays can be from 
;               1 - 9 values in length (missing trailing values are assumed 0)
;               specifying: 
;                   [0]: piston = 1.
;                   [1]: x-tilt = r cos(theta)
;                   [2]: y-tilt = r sin(theta)
;                   [3]: defocus = 2r^2-1
;                   [4]: 0-90-astigmatism = r^2 cos(2 theta)
;                   [5]: 45-135-astigmatism = r^2 sin(2 theta)
;                   [6]: x-coma = (3r^2-2)r cos(theta)
;                   [7]: y-coma = (3r^2-1)r sin(theta)
;                   [8]: spherical aberration = 6r^4-6r^2+1
;                   [9]: x-trefoil: r^3 cos(3 theta)
;                  [10]: y-trefoil: r^3 sin(3 theta)
;                In the absence of the higher order wavefront aberrations:
;              The tilt coefficients are releated to the tilt (alpha) of the 
;              mirror by Z[1],Z[2] = D*alpha/lambda, where D = aperture 
;              diameter. The image displacement will be 2*alpha.
;              The defocus coefficient is related to the the distance from 
;              the true focus (z) by Z[3] = z/(16*(f#)^2*lambda), where
;              f# = focal length/aperture diameter.
;                I'm calculating the phases as exp(2*!pi*i*Z). If I really 
;              should be using exp(i*Z) then the coeffecients as used here are 
;              a factor of 2*!pi smaller than they should be. 

; OUTPUT
; beam0       = this contains the ideal beam produced by the aperture


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; This is a very stripped down version of optics10.pro, which can do the same 
; calculation, but will be slower.
; Only the essential features for computing the IRAC ghost images are retained, 
; and assumptions are made that the input will be appropriate.
; Roughly half the CPU time is spent doing the FFT. 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;-- convert physical units to pixels for mapping

asec_um = 3600.*!radeg/(D*f_ratio)*1e-6
asec_pix = um_pix*asec_um
aprad = asec_pix*nsize*D/lambda/2./3600*!dtor
;spider = aprad*spiderfrac[0]

;-- calculate aperture masks

x = lindgen(nsize,nsize) mod nsize
y = lindgen(nsize,nsize)  /  nsize
x0 = x-nsize/2
y0 = y-nsize/2
th = shift(atan(y0,x0),nsize/2,nsize/2)
r = dist(nsize)/aprad

aperture = r le 1. and r ge obstfrac

nspid  = n_elements(spiderfrac)
for i=1,nspid-1 do begin 
  th0 = spiderfrac[i] mod (2*!pi)
  th0 = th0 -2*!pi*(th0 gt !pi)
  aperture=aperture and (abs(th-th0) ge -spiderfrac[0])
  if (th0-spiderfrac[0] gt !pi) then $
    aperture=aperture and (abs(th-th0+2*!pi) ge -spiderfrac[0])
  if (th0+spiderfrac[0] le -!pi) then $
    aperture=aperture and (abs(th-th0-2*!pi) ge -spiderfrac[0])
endfor

; -- apply shifts, defocus, and trefoil aberration

zc = fltarr(11)
zc[0:n_elements(a_zcoeff)-1] = a_zcoeff
r2 = r^2
r3 = r^3
sinth = sin(th)
costh = cos(th)
aperture=aperture*exp((2*!pi*complex(0,1))* $
  (zc[0]+zc[1]*r*costh+zc[2]*r*sinth+zc[3]*(2*r2-1)+ $
   zc[4]*r2*cos(2*th)+zc[5]*r2*sin(2*th)+ $
;   zc[6]*r*(3*r2-2)*costh+zc[7]*r*(3*r2-2)*sinth+ $
;   zc[8]*(6*r^4-6*r2+1)+zc[9]*r3*cos(3*th)+$
   zc[10]*r3*sin(3*th)))

;-- calculate FFTs
beam0 = fft(aperture,-1)*nsize

end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; ghostmasker1234 adapted from Rick's final ghostbuster code.

function ghostmasker1234, channel, x_source_pix, y_source_pix, $
  defocus_ghost=defocus_ghost, astigma4_ghost=astigma4_ghost, $
  astigma5_ghost=astigma5_ghost

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; INPUT:
; channel      = (INT) 1 or 2, for IRAC 3.6 or 4.5 um channels
; x_source_pix = (FLOAT) x-position of point source on array (units: pixels)
; y_source_pix = (FLOAT) y-position of point source on array (units: pixels)
;
; OPTIONAL INPUT KEYWORDS:
; defocus_ghost  = Zernike coefficient for defocussing Ch. 3 or 4 ghost image.
; astigma4_ghost = Zernike coeff. for 0-90 deg astigmatism of Ch. 3 or 4 ghost
; astigma5_ghost = Zernike coeff. for 45-135 deg astigmatism of Ch. 3 or 4 ghost

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; PROCEDURES CALLED:
; OPTICS10_GHOST - A stripped-down version of my "OPTICS10" FFT beam modeling 
; program tailored specifically to use here.

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; This function calculates the ghost image for an IRAC Channel 1 (3.6um), 
; Channel 2 (4.5 um), Channel 3 (5.8 um), or Channel 4 (8um)
; point source at the specified location. 

; The ghost image is computed by calculating the FFT of the Spitzer 
; aperture (which alone would yield the beam), with aberrations to shift and
; defocus the image. A slight blurring is also applied as reported
; in the Hoffman et al. 2004 SPIE paper.

; The returned image is normalized such that multiplying it by the integrated 
; flux of the source should yield an appropriately scaled ghost image. 
; The main difficulty here is that obviously ghosted sources are often 
; saturated, and determining their true integrated flux may be difficult.

; Notes:
; 1) Ghosts should be removed from single frames before mosiacking.
; 2) Ghost positions and shapes are not perfectly accurate. The aberrations 
; were fit by eye, and were not quantitatively optimized. 
; 3) The burden of determining source position and intensity is placed on the 
; user. It would be nice to have the code calculate these automatically given
; a single frame, but this really the business of a point source extractor 
; and I choose not to embed one here.
; 4) The code is slow. For each ghost, it takes ~2 seconds on my 
; dual-1.8GHz G5 Mac. Roughly half of that is spent doing the FFT.
; The fastest shortcut would probably be to read a pre-computed ghost image
; and simply interpolate that image to the desired location on the array.
; 5) By editing the value of nsamples, the code can average the ghost over
; a specified number of wavelengths within the band. However, the code will
; run proportionally slower, and present results do not suggest that such 
; additional effective smoothing of the ghost image is necessary.
; 6) The default shapes for the Ch. 3 and 4 ghosts don't seem quite right
; for some cases, which is why the code currently has optional defocus and 
; astigmatism keywords. Adding a small defocus changes a "+" shaped astigmatic
; ghost into a more linear "-" shape or "|" shape depending upon the sign
; of the defocus. 

; Version 1.0 : R. Arendt 2004/09/23
; version 1.1 : R. Arendt 2004/11/08 - added astigmatic ghosts for Ch. 3 & 4

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; set default values (not really very useful to call without parameters)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

if (n_elements(channel) ne 1) then channel = 1
if (n_elements(x_source_pix) ne 1) then x_source_pix = 128.
if (n_elements(y_source_pix) ne 1) then y_source_pix = 128.

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; define array sizes
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

n = 1024 ; 10'x10' array to avoid aliasing
um_pix = 15. ; calculate on 1/2 size pixels

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Telescope parameters
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

d = 0.85
f = 6. ; telescope is f/12, but IRAC reimages to f/6.
obscur = 0.3765
spider = [-0.035,2*!pi/12*[3,7,11]+[+0.002,-0.015,-0.005]]

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Channel-dependent parameters
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

case (channel) of
  1: begin
      wave0 = 3.58e-6
      bandwidth = 0.750e-6
      scale_ghost = 0.005
      defocus_ghost = 0.8
      astigma4_ghost = 0
      astigma5_ghost = 0
      trefoil_ghost =0.10
      blur = [0.43,0.43] * (30/um_pix)
      ;; Fit to C. Marx's models (quoted in SOM)
      ;shift_ghost = ([x_source_pix-128,y_source_pix-128] + $
      ;         [-0.04351*(255-x_source_pix)/256.-0.00288,$
      ;           0.04761*y_source_pix/256.+0.00211]*256 )*30e-6
      ; Fit to 155 hand measurements on data from GlIMPSE AOR 9234432
      ; This seems to work better. Differs by 1-3 pix in x, 0-0.2 pix in y
      shift_ghost = ([x_source_pix-128,y_source_pix-128] + $
               [0.049686*x_source_pix-29.222/2,$
                0.048301*y_source_pix+1.0921/2])*30e-6
     end
  2: begin
      wave0 = 4.52e-6
      bandwidth = 1.015e-6
      scale_ghost = 0.008
      defocus_ghost = 0.56
      astigma4_ghost = 0
      astigma5_ghost = 0
      trefoil_ghost =0.10
      blur = [0.37,0.37] * (30/um_pix)
      ;; Fit to C. Marx's models (quoted in SOM, but with typo) 
      ;shift_ghost = ([x_source_pix-128,y_source_pix-128] + $
      ;         [ 0.04956*(255-x_source_pix)/256.+0.00105,$
      ;           0.04964*y_source_pix/256.+0.000387]*256 )*30e-6
      ; Fit to 137 hand measurements on data from GlIMPSE AOR 9234432
      ; Measured these as for Ch. 1.  Differs by 1-1.5 pix in x, +/-0.5 pix in y
      shift_ghost = ([x_source_pix-128,y_source_pix-128] + $
               [0.051957*x_source_pix+2.4135/2,$
                0.053381*y_source_pix-0.90516/2])*30e-6
     end
  3: begin
      wave0 = 5.693e-6
      bandwidth = 1.425e-6
      scale_ghost = 0.001
      if (n_elements(defocus_ghost) ne 1) then defocus_ghost = 0.00
      if (n_elements(astigma4_ghost) ne 1) then astigma4_ghost = -0.4
      if (n_elements(astigma5_ghost) ne 1) then astigma5_ghost = 0.03
      trefoil_ghost =0.0
      blur = [0.65,0.65] * (30/um_pix)
      ; Fit to [42,42] hand measurements on data from GlIMPSE AOR 9234432
;      shift_ghost = [-35.90+x_source_pix-128,$
;        0.036478521*y_source_pix-4.8522816/2+y_source_pix-128]*30e-6
      shift_ghost = [-34.90+x_source_pix-128,$
        0.036478521*y_source_pix-4.8522816/2+y_source_pix-128]*30e-6
     end
  4: begin
      wave0 = 7.982e-6
      bandwidth = 2.905e-6
      scale_ghost = 0.001
      if (n_elements(defocus_ghost) ne 1) then defocus_ghost = -0.10
      if (n_elements(astigma4_ghost) ne 1) then astigma4_ghost = -0.35
      if (n_elements(astigma5_ghost) ne 1) then astigma5_ghost = -0.05
      trefoil_ghost =0.0
      blur = [0.36,0.36] * (30/um_pix)
      ; Fit to [46,48] hand measurements on data from GlIMPSE AOR 9234432
      shift_ghost = [35.78+x_source_pix-128,$
        0.037653032*y_source_pix-3.6676230/2+y_source_pix-128]*30e-6
     end
endcase

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; calculation of ghost image
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

nsamples = 1
gpsf = fltarr(n,n)
for i = -(nsamples-1.)/2.,(nsamples-1.)/2. do begin
  wave = wave0+i*bandwidth/nsamples
  optics10_ghost,n,wave,d,f,obscur,spider,um=um_pix,beam0=ghost,$
    a_zcoeff=[0,shift_ghost[0]/(2*f*wave),shift_ghost[1]/(2*f*wave),$
    defocus_ghost,astigma4_ghost,astigma5_ghost,0,0,0.,0,trefoil_ghost]
  gpsf = gpsf + abs(temporary(ghost))^2
endfor

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;  smooth, extract, and rescale detector-scale region  
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

kshift = fix(blur*5)
kwidth = kshift*2+1
xx = (indgen(kwidth[0],kwidth[1]) mod kwidth[0]) - kshift[0]
yy = (indgen(kwidth[0],kwidth[1]) / kwidth[0]) - kshift[1]
kern = exp(-0.5*((xx/blur[0])^2+(yy/blur[1])^2))
kern = kern/total(kern)

gpsf = shift(gpsf,256,256)
gpsf = convol(gpsf,kern,/edge_wrap)
gpsf = congrid(gpsf[0:511,0:511],256,256)
gpsf = gpsf*(scale_ghost/nsamples/total(gpsf,/double))

; this part makes a mask from the ghostimage, scaled to typical dynamic
; range required for masking the ghosts (100000). Mask bit for optical 
; ghosts is 2.

gmptr = where(gpsf gt 10.0e-6)

gmask = intarr(256,256)

if gmptr[0] ge 0 then gmask[gmptr] = 4

return, gmask

end




pro mask_straylight_v5p2,bcd_list,dmask_list,tmass_list,CHDMSK=chdmsk,UNAGGR=unaggr

; Driver for straylight masking program. Uses masking subroutines
; to make a straylight mask given a list of BCDs and a 2mass list
; covering the FIF. See top of file for details of I/O.

; explain what's being written out

print,'Program sets bit 3 (8) in areas likely to be affected by stray light'
print,'Program sets bit 2 (4) in areas likely to have a filter ghost'

; check for dmask overwriting keyword

if (keyword_set(chdmsk)) then print,'Will write masks into dmasks/imasks'

; smallonly determines whether or not the large masks are used

smallonly = 0

if (keyword_set(unaggr)) then smallonly = 1
if (smallonly eq 1) then print,'Using unaggressive option. (Small masks only)'

; File I/O:

bcd_list = strcompress(bcd_list,/remove_all)
tmass_list = strcompress(tmass_list,/remove_all)
dmask_list = strcompress(dmask_list,/remove_all)

readcol,bcd_list,filename,format='a'
readcol,tmass_list,ra,dec,emaj,emin,eang,jmag,hmag,kmag,$
format='d,d,f,f,i,f,f,f',skipline=9
readcol,dmask_list,dmaskname,format='a'

nframes = n_elements(filename)
nstars  = n_elements(ra)

; Loop over files:

for i=0,nframes-1 do begin

; read in FITS file, extract astrometry from header and 
; determine channel number and exposure time

  print,'reading file ',filename[i]

  image = readfits(filename[i],hdr,/silent)

  dmask = readfits(dmaskname[i],dhdr,/silent)

  chan = fix(sxpar(hdr,'CHNLNUM'))

  exptime = sxpar(hdr,'EXPTIME')

  framtime = sxpar(hdr,'FRAMTIME')

; according to Bill G best exposure time estimate is mean of frametime
; and exptime:

  otime = (exptime+framtime)/2.

; define the blank mask for each frame

  tmask=bytarr(256,256)

; Depending on which channel the data is from, set the thresholding 
; levels and call the appropriate procedures. Procedures ending fpa
; mask FPA hole scattering, those ending det detector scattering. 
; Plus 1c and 3s scattering zones which are probably due to other mechanisms.
; Two thresholding levels (for mask and mask2) are set for each different 
; mechanism, with a similar naming convention.

  case (chan) of

    1: begin 

;      Set the thresholding levels and adjust for exposure time:

	ch1ghst = 8.7 + 2.5*alog10(otime/11.2)
	ch1fpa  = 9.3 + 2.5*alog10(otime/11.2)
	ch1fpa2 = 7.2 + 2.5*alog10(otime/11.2)
	ch1c    = 7.2 + 2.5*alog10(otime/11.2)
	ch1c2   = 4.9 + 2.5*alog10(otime/11.2)

;      Run through 2mass list for stars brighter than the faintest scattering 
;      threshold

	  for j=0,nstars - 1 do begin

	    jmk = jmag[j] - kmag[j]

	    ch1 = kmag[j] - 0.298*jmk

;      Check the faintest scattering threshold

	    if (ch1 le ch1fpa) then begin

;      adxy converts RA, Dec to x,y coords in the BCD's frame

            adxy,hdr,ra[j],dec[j],x,y

;      create blank masks to feed the subroutines

            mask = bytarr(256,256)
            mask2= bytarr(256,256)

; first check if the stars have ghosts (only for stars within 10 pixels
; of the FPA)

       if ((x gt -10) and (y gt -10) and (x le 265) and (y le 265) and ch1 le ch1ghst) then begin



	    gmask = ghostmasker1234(1,x[0,0],y[0,0])


   	if total(gmask) gt 0 then print,'Adding channel 1 ghost mask..'

            tmask = tmask or gmask

	    endif

;      only process stars with (x,y) within 100 pixels of edge of FPAs

            if (abs(x-128) lt 228) and (abs(y-128) lt 228) and $
            ((abs(x-128) gt 128) or (abs(y-128) gt 128)) then begin

;       Call the appropriate straymask programs:
;
;       (Apparently older versions of IDL and/or ADXY (or AD2XY) will return X
;       and Y as single-element arrays. Redefining them as 2-d arrays here 
;       prevents the straymaskN routines from calculating 
;       array[256,256]-array[1] = array[1] )
;

	    straymask1fpa,x[0,0],y[0,0],mask,mask2,hdr

;       Add the new mask:

          tmask = tmask or 8*mask

; if the star is (very bright) above the threshold for the 2nd mask add it in

          if (ch1 le ch1fpa2) and (smallonly eq 0) then tmask = tmask or 8*mask2

;       reset the mask arrays, call the next straymask routine:

          mask = mask*0
	  mask2= mask2*0

          straymask1c,x[0,0],y[0,0],mask,mask2,hdr

; if the star is bright enough, add this mask

	  if ((ch1 le ch1c) and (total(mask) gt 0))  then begin

	  tmask = tmask or 8*mask

          print,'Channel 1 zone 1C scatterer added to mask'

	 endif

; if the star is incredibly bright add mask 2:

	  if (ch1 le ch1c2) and (smallonly eq 0) then tmask = tmask or 8*mask2

        endif  ; end of work for stars close to the FPA

      endif  ; end of work for bright stars

     endfor ; end loop over stars

   end ; end case of chan 1

    2: begin 

;      Set the thresholding levels and adjust for exposure time:

	ch2ghst = 7.4 + 2.5*alog10(otime/11.2)
	ch2fpa  = 8.8 + 2.5*alog10(otime/11.2)
	ch2fpa2 = 6.8 + 2.5*alog10(otime/11.2)

;      Run through 2mass list for stars brighter than the faintest scattering 
;      threshold

	  for j=0,nstars - 1 do begin

	  jmk = jmag[j] - kmag[j]
	  ch2 = kmag[j] - 0.387*jmk

;      Check the faintest scattering threshold

	    if (ch2 le ch2fpa) then begin

;      adxy converts RA, Dec to x,y coords in the BCD's frame

            adxy,hdr,ra[j],dec[j],x,y

;      create blank masks to feed the subroutines

            mask = bytarr(256,256)
            mask2= bytarr(256,256)


; first check if the stars have ghosts (only for stars within 10 pixels
; of the FPA)

      if ((x gt -10) and (y gt -10) and (x le 265) and (y le 265) and ch2 le ch2ghst) then begin

	    gmask = ghostmasker1234(2,x[0,0],y[0,0])

   	if total(gmask) gt 0 then print,'Adding channel 2 ghost mask..'

            tmask = tmask or gmask

       endif

;      only process stars with (x,y) within 100 pixels of edge of FPAs

            if (abs(x-128) lt 228) and (abs(y-128) lt 228) and $
            ((abs(x-128) gt 128) or (abs(y-128) gt 128)) then begin

;       Call the appropriate straymask programs:
;
;       (Apparently older versions of IDL and/or ADXY (or AD2XY) will return X
;       and Y as single-element arrays. Redefining them as 2-d arrays here 
;       prevents the straymaskN routines from calculating 
;       array[256,256]-array[1] = array[1] )
;

	    straymask2fpa,x[0,0],y[0,0],mask,mask2,hdr

;         Add the new mask:

            tmask = tmask or 8*mask

; if the star is (very bright) above the threshold for the 2nd mask add it in

            if (ch2 le ch2fpa2) and (smallonly eq 0) then tmask = tmask or 8*mask2

        endif  ; end of work for stars close to the FPA

      endif  ; end of work for bright stars

     endfor ; end loop over stars

   end ; end case of chan 2


    3: begin 

;      Set the thresholding levels and adjust for exposure time:

	ch3ghst = 7.0 + 2.5*alog10(otime/11.2)
	ch3det  = 7.8 + 2.5*alog10(otime/11.2)
	ch3det2 = 6.5 + 2.5*alog10(otime/11.2)
	ch3fpa  = 4.8 + 2.5*alog10(otime/11.2)
	ch3fpa2 = 2.8 + 2.5*alog10(otime/11.2)
	ch3s    = 4.8 + 2.5*alog10(otime/11.2)
	ch3s2   = 2.8 + 2.5*alog10(otime/11.2)

;      Run through 2mass list for stars brighter than the faintest scattering 
;      threshold

        for j=0,nstars - 1 do begin

	  jmk = jmag[j] - kmag[j]
	  ch3 = kmag[j] - 0.387*jmk

;      Check the faintest scattering threshold

	    if (ch3 le ch3det) then begin

;      adxy converts RA, Dec to x,y coords in the BCD's frame

            adxy,hdr,ra[j],dec[j],x,y

;      create blank masks to feed the subroutines

            mask = bytarr(256,256)
            mask2= bytarr(256,256)

; first check if the stars have ghosts (only for stars within 10-40 pixels
; of the FPA)

      if ((x gt -40) and (y gt -10) and (x le 295) and (y le 265) and ch3 le ch3ghst) then begin

	    gmask = ghostmasker1234(3,x[0,0],y[0,0])

   	if total(gmask) gt 0 then print,'Adding channel 3 ghost mask..'

            tmask = tmask or gmask

       endif

;      only process stars with (x,y) within 100 pixels of edge of FPAs

            if (abs(x-128) lt 228) and (abs(y-128) lt 228) and $
            ((abs(x-128) gt 128) or (abs(y-128) gt 128)) then begin

;       Call the appropriate straymask programs:
;
;       (Apparently older versions of IDL and/or ADXY (or AD2XY) will return X
;       and Y as single-element arrays. Redefining them as 2-d arrays here 
;       prevents the straymaskN routines from calculating 
;       array[256,256]-array[1] = array[1] )
;

	    straymask3det,x[0,0],y[0,0],mask,mask2,hdr

;       Add the new mask:

          tmask = tmask or 8*mask

; if the star is (very bright) above the threshold for the 2nd mask add it in

          if (ch3 le ch3det2) and (smallonly eq 0) then tmask = tmask or 8*mask2

;       reset the mask arrays, call the next straymask routine:

          mask = mask*0
	  mask2= mask2*0

          straymask3fpa,x[0,0],y[0,0],mask,mask2,hdr

; if the star is bright enough, add this mask

	  if ((ch3 le ch3fpa) and (total(mask) gt 0)) then begin

		tmask = tmask or 8*mask
		print,'Added channel 3 FPA scatterer mask'

	  endif

; if the star is incredibly bright add mask 2:

	  if (ch3 le ch3fpa2) and (smallonly eq 0) then tmask = tmask or 8*mask2

;       reset the mask arrays, call the next straymask routine:

          mask = mask*0
	  mask2= mask2*0

          straymask3s,x[0,0],y[0,0],mask,mask2,hdr

; if the star is bright enough, add this mask

	  if ((ch3 le ch3s) and (total(mask) gt 0)) then begin

		tmask = tmask or 8*mask
		print,'Adding channel 3 conic section scatterer to mask'

	  endif

; if the star is incredibly bright add mask 2:

	  if (ch3 le ch3s2) and (smallonly eq 0) then tmask = tmask or 8*mask2

        endif  ; end of work for stars close to the FPA

      endif  ; end of work for bright stars

     endfor ; end loop over stars

   end ; end case of chan 3

    4: begin

;      Set the thresholding levels and adjust for exposure time:

	ch4ghst = 6.4 + 2.5*alog10(otime/11.2)
	ch4det  = 8.0 + 2.5*alog10(otime/11.2)
	ch4det2 = 5.7 + 2.5*alog10(otime/11.2)
	ch4fpa  = 5.5 + 2.5*alog10(otime/11.2)
	ch4fpa2 = 4.5 + 2.5*alog10(otime/11.2)

;      Run through 2mass list for stars brighter than the faintest scattering 
;      threshold

	  for j=0,nstars - 1 do begin

		jmk = jmag[j] - kmag[j]

		ch4 = kmag[j] - 0.387*jmk

;      Check the faintest scattering threshold

		    if (ch4 le ch4det) then begin

;      adxy converts RA, Dec to x,y coords in the BCD's frame

        	    adxy,hdr,ra[j],dec[j],x,y

;      create blank masks to feed the subroutines

	            mask = bytarr(256,256)
	            mask2= bytarr(256,256)

; first check if the stars have ghosts (only for stars within 10-40 pixels
; of the FPA)

      if ((x gt -30) and (y gt -10) and (x le 295) and (y le 265) and ch4 le ch4ghst) then begin

	    gmask = ghostmasker1234(4,x[0,0],y[0,0])

   	if total(gmask) gt 0 then print,'Adding channel 4 ghost mask..'

            tmask = tmask or gmask

       endif

;      only process stars with (x,y) within 100 pixels of edge of FPAs

            if (abs(x-128) lt 228) and (abs(y-128) lt 228) and $
            ((abs(x-128) gt 128) or (abs(y-128) gt 128)) then begin

;       Call the appropriate straymask programs:
;
;       (Apparently older versions of IDL and/or ADXY (or AD2XY) will return X
;       and Y as single-element arrays. Redefining them as 2-d arrays here 
;       prevents the straymaskN routines from calculating 
;       array[256,256]-array[1] = array[1] )
;

	    straymask4det,x[0,0],y[0,0],mask,mask2,hdr

;       Add the new mask:

          tmask = tmask or 8*mask

; if the star is (very bright) above the threshold for the 2nd mask add it in

          if (ch4 le ch4det2) and (smallonly eq 0) then tmask = tmask or 8*mask2

;       reset the mask arrays, call the next straymask routine:

          mask = mask*0
	  mask2= mask2*0

          straymask4fpa,x[0,0],y[0,0],mask,mask2,hdr

; if the star is bright enough, add this mask

	  if ((ch4 le ch4fpa) and (total(mask) gt 0)) then begin

	     tmask = tmask or 8*mask

	     print,'Adding channel 4 FPA scatterer to mask'

	  endif

; if the star is incredibly bright add mask 2:

	  if (ch4 le ch4fpa2) and (smallonly eq 0) then tmask = tmask or 8*mask2

        endif  ; end of work for stars close to the FPA

      endif  ; end of work for bright stars

     endfor ; end loop over stars

   end ; end case of chan 4

  endcase

; Now decide whether to add the straymask to the dmask, or write out
; a separate smask file.

if keyword_set(chdmsk) then begin

; Alter the dmasks is chdmsk is set.
; Use bit 3 (8) as the mask bit as this is electrical cross-talk, 
; so unlikely ever to be implemented in the dmask.
; Use of OR eliminates problems if the bit is already set anyway.

     if (total(tmask) gt 0) then print,'Straylight mask added to dmask file'

        dmask = dmask or tmask

        writefits,dmaskname[i],dmask,dhdr

        writefits,filename[i],image,hdr

     endif else begin

	smaskname = 'smask_'+filename[i]

	print,'writing scattered light mask to smask file'

	writefits,smaskname,tmask,hdr

	writefits,filename[i],image,hdr

     endelse

endfor ; end loop over BCD files

end