pro phot, data, starrad, skyinrad, skyoutrad, filename
;
;  Photometry Program 
;  Measures stars in array data (assumed to be displayed in window 0
;  Parameters are the data array, the stellar aperture, the inner and 
;     outer sky radii, and the output filename
;     
;
np = n_params()
if np le 3 then goto, ERRMES
sz = size(data)
if sz(0) ne 2 then goto, ERRMES
;
;  Set up for circles (positive or negative background)
;
grey = 0
QUESTION:
answer=' '
print,format="('Do you have a color screen?',$)"
read,answer
;read,answer,prompt='Do you have a color screen?:  '
if answer eq 'y' or answer eq 'Y' then grey=-1
if answer eq 'n' or answer eq 'N' then grey=1
if grey eq 0 then goto, QUESTION
color = 255
; if grey * !neg gt 0  then color = 255
; if grey * !neg lt 0  then color = 123
;  
;  Make a window for image and surface plot
;
if answer eq 'y' then begin
    window,2,xpos=280,xsize=280,ysize=300,title='Surface Plot'
    window,4,xpos=00,xsize=280,ysize=300,title='Radial Profile'
endif else begin
    window,2,xpos=375,xsize=300,ysize=400,title='Surface Plot'
    window,4,xpos=60,xsize=300,ysize=400,title='Radial Profile'
endelse 
set_viewport,.2,.9,.2,.9
wset,0
centroid=5
;
if np eq 5 then begin
    close,/all
    i = strlen(findfile(filename))
    openw,luny,filename,/get_lun,/append
    if i(0) eq 0 then begin
        printf,luny, ' xpos    ypos    starflux   skyflux   starmag    comment'
        printf,luny,' '
    endif
endif
print,'Place cursor at star position and click...'
print,string(7B),' '
print,string(7B),' '
print,string(7B),' '
print,string(7B)
print, ' xpos    ypos    starflux   skyflux   starmag'
print,string(7B)

;begin main loop for reading cursor and plotting 
;
comment=' '
repeat begin
    cursor,xdat,ydat,/device
    getbox, data, xdat, ydat, centroid, i1, i2, j1, j2
    twodcen, data, i1, i2, j1, j2, xdat, ydat
    getbox, data, xdat, ydat, centroid, i1, i2, j1, j2
    twodcen, data, i1, i2, j1, j2, xc, yc
    getbox, data, xc, yc, skyoutrad, i1, i2, j1, j2
    loadpix, data, xc, yc, i1, i2, j1, j2, dist_arr, data_arr
    srad = where (dist_arr ge skyinrad and dist_arr le skyoutrad)
    skyflux = median(data_arr(srad))
    srad = where (dist_arr le starrad)
    starflux = total(data_arr(srad))
    netflux = starflux - (n_elements(data_arr(srad))*skyflux)
    starmag = 26. - (2.5 * alog10(netflux))
    tvcircle,starrad,xc,yc,color
    wset,2
    surface,data(i1:i2,j1:j2)
    wset,4
    !psym=3
    ymax=max(data_arr)
    plot,dist_arr,data_arr,xrange=[0,skyoutrad],yrange=[0,ymax],xstyle=1, $
          ystyle=1, xtitle='Distance',ytitle='Counts'
          
    wset,0
    print,xc,yc,starflux,skyflux,starmag,format='(f6.2,f8.2,f11.1,f10.2,f9.2)'
    print,format="('Comment (or type STOP)',$)"
    read,comment
;
;    read,comment,prompt='Comment (or type STOP):  '
    if np eq 5 then begin
        printf,luny,xc,yc,starflux,skyflux,starmag,comment,format= $
                 '(f6.2,f8.2,f11.1,f10.2,f9.2,5x,a)'
    endif
    if comment eq 'STOP' then goto, DONE
;
endrep until 0
ERRMES:
print,' ERROR! This routine is called by'
print,string(7B)
print,' phot, data_array, aperture, in_sky_ap, out_sky_ap, filename'
print,string(7B)
print,' '
return
DONE:
close,/all
return
end
;
pro getbox, data, xdat, ydat, boxrad, i1, i2, j1, j2
i1 = fix(xdat - boxrad)
i2 = fix(xdat + boxrad)
j1 = fix(ydat - boxrad)
j2 = fix(ydat + boxrad)
if i1 lt 0 then i1 = 0
if j1 lt 0 then j1 = 0
sz = size(data)
if i2 ge sz(1) then i2 = sz(1) - 1
if j2 ge sz(2) then j2 = sz(2) - 1
return
end
;
;
pro twodcen, data, i1, i2, j1, j2, xc, yc
ic=i1
jc=j1
for i = i1, i2 do begin
   for j = j1, j2 do begin
	if data(i,j) gt data(ic,jc) then begin
	    ic = i
	    jc = j
	endif
   endfor
endfor
ioff = min([i2-ic,ic-i1])
joff = min([j2-jc,jc-j1])
sumx = 0.
sumy = 0.
tsum = 0.
for i = ic-ioff, ic+ioff do begin
   for j = jc-joff, jc+joff do begin
	sumx = sumx + (data(i,j) * i)
	sumy = sumy + (data(i,j) * j)
	tsum = tsum + data(i,j)
   endfor
endfor
xc = sumx / tsum
yc = sumy / tsum
return
end
;
;
pro loadpix, data, xc, yc, i1, i2, j1, j2, dist_arr, data_arr
num = (i2-i1+1)*(j2-j1+1)
data_arr = fltarr(num)
dist_arr = fltarr(num)
k = 0
for i = i1, i2 do begin
    for j = j1, j2 do begin
	dist = ((i-xc)^2) + ((j-yc)^2)
	dist_arr(k) = sqrt(dist)
	data_arr(k) = data(i,j)
	k = k + 1
    endfor
endfor
return
end