;-------------------------------------------------------------------------------
; IDL Procedure PHOT_NSFSUBARRAY.PRO
;
; This IDL procedure takes de-dithered, flatfielded observation made at the IRTF ; and places a box of width 2*APSIZE+1 around the star. A frame of pixels is 
; created around the object, with outer width of 4*APSIZE+1 and inner width of
; 2*APSIZE+1. The total counts within the box minus the (average counts/pixel in ; the frame) * (2*APSIZE+1)^2 are then ouput as the background corrected
; flux from the object, as well as the noise/pixel estimate given by the value
; of Sigma(Frame) * sqrt((2*APSIZE+1)^2).
; 
; CML 2/29/96
;-------------------------------------------------------------------------------

pro phot_nsfsubarray,array,apsize,framesize,xguess,yguess,xpos,ypos,fsvalue,bkvalue,bknoise

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

;Estimate Centroid If Not Already Given - Get Close Enough For Routine to Work

if((xguess eq 0) and (yguess eq 0)) then begin
	window,0,xsize=2*256,ysize=2*256,title='Object Image + Contour'
	tvscl,alog10(array>0+.01)
	center,array,xguess,yguess 			
	print,xguess,yguess
endif

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


;Make Best Estimate of Centroid

;Use for Phaethon
findmax,array,'max',xguess-5,xguess+5,yguess-5,yguess+5,xpos,ypos,maxval

;Use for HMP
;findmax,array,'max',xguess-30,xguess+30,yguess-30,yguess+30,xpos,ypos,maxval
;findmax,array,'med',xguess-30,xguess+30,yguess-30,yguess+30,xpos,ypos,maxval
;print,'Xcenter = ',xpos,'         Ycenter = ',ypos

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

;Estimate Signal = Sum Over Aperature of Radius Apsize

theta=findgen(720)/720.*2.*!pi
x = xpos+(apsize*cos(theta))
y = ypos+(apsize*sin(theta))
appix=(polyfillv(x,y,256,256))
fsvalue = total(array(appix))
Napp = n_elements(appix)
;print,'Apsize = ',apsize,' # of pixels = ',Napp

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

;Create Frame of Data Around Object For Background Analysis

d      = framesize	
;Largest Value Possible w/ HMP_LC.PRO Data - Otherwise Off Chip

goto,hmp
;Left Vertical Slat
array1  = array(xpos-4*l:xpos-2*l,ypos-4*l:ypos+4*l) & array1=array1(0:*)
;Right Vertical Slat
array2  = array(xpos+2*l:xpos+4*l,ypos-4*l:ypos+4*l) & array2=array2(0:*)
;Lower Slat
array3  = array(xpos-2*l:xpos+2*l,ypos-4*l:ypos-2*l) & array3=array3(0:*)
;Upper Slat
array4  = array(xpos-2*l:xpos+2*l,ypos+2*l:ypos+4*l) & array4=array4(0:*)
Frame   = [array1,array2,array3,array4]

hmp:
Frame1   = array(1:75,65:75) & Frame1 = Frame1(0:*)
Frame2   = array(1:75,5:15) & Frame2 = Frame2(0:*)
Frame    = [Frame1,Frame2]

Nframe  = n_elements(Frame)
;print,'Nframe = ',nframe

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

;Estimate Background in Signal = Background/pixel * (# of Pixels in Comet Box)
;Estimate Total Noise on Signal = (Noise/Pixel) * sqrt(# of pixels in Comet Box)

resistant_mean,Frame,3.0,bkvalue,bknoise,num_rej
bkvalue = bkvalue * Napp
bknoise = bknoise * sqrt(Nframe) * sqrt(Napp)	
					; N.B. - Sigma Returned By Res_Mean
					; is Sigma of the Mean, equal to
					; Sigma of the Dist/sqrt(Nframe)

;print,'Resistant Mean(Frame) = ',bkvalue/Napp
;print,'Median(Frame) = ',med(frame)

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

;Subtract Background From Observation
fsvalue = fsvalue - bkvalue

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

return
end