;-------------------------------------------------------------------------------
; IDL Procedure PHOT.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_msx,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
arraysize=size(array)
xsize=arraysize(1)
ysize=arraysize(2)
angstep = apsize * 4.0 * !pi

theta=findgen(angstep)/angstep*2.*!pi
x = xpos+(apsize*cos(theta))
y = ypos+(apsize*sin(theta))
appix=(polyfillv(x,y,xsize,ysize))
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	

;Lower Slat
Frame1 = array(5:xsize-5,5:15) & Frame1 = Frame1(0:*)
;Upper Slat
Frame2 = array(5:xsize-5,ysize-10:ysize-5) & Frame2 = Frame2(0:*)
Frame    = [Frame1,Frame2]
;Left Vertical Slat
Frame3  = array(5:10,5:ysize-5) &  Frame3 = Frame3(0:*)
;Right Vertical Slat
Frame4  = array(xsize-10:xsize-5,5:ysize-5) &  Frame4 = Frame4(0:*)

Frame   = [Frame1,Frame2,Frame3,Frame4]

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