pro defocus_star_phot, infile_list, output_summary_file
;
;   performs aperture photometry on a single star in a group of
;   images. This version works on defocussed stars, as generated
;   when doing polarimetric calibration on standards
;
;   DPC     20051026    original code
;           20051223    changed bad pixel replacer to call to fix_bad_pixels.pro
;           20051228    reformed entire program
;
conversion_gain = 10.0   ; e/ADU
read_noise = sqrt(2.0) * 21.   ; appropriate for single CDS images
;
np = N_PARAMS()
if(np eq 0) then begin
    files = dialog_pickfile(Title="Select Image Files or Input List File",/MULTIPLE_FILES)
    nf = SIZE(files)
    nfiles = nf[1]
    ilist = 0
    trial_name = ''
    if(nfiles eq 1) then begin
        ;
        ;  if only one file indicated, parse name to see if image or list
        ;
        test = STRMATCH(files[0],'FITS',/FOLD_CASE)
        ;
        if(test eq 0) then begin       ; was not fits file
           ;
           ; is a list file name, set pointer
           ;
           ilist = 1
           infile_list = files[0]
           trial_name = FILE_BASENAME(files[0],".dat") + '_phot.dat'
        endif
    endif
    file_out = dialog_pickfile(Title="Select Output Summary File",file=trial_name)
endif else begin
    ilist = 1
    ;
    file_out = output_summary_file
endelse
;
;   read list if flag set
;
if(ilist eq 1) then begin
    ;
    ;  read file list
    ;
    openr, lu_in, infile_list, /GET_LUN
    files = strarr(300)
    nfiles = 0
    while (~EOF(lu_in)) do begin
       line = ''
       readf, lu_in, line
       files[nfiles] = STRTRIM(line,2)
       nfiles++
    endwhile
    FREE_LUN, lu_in
    files = files[0:nfiles-1]
    ;
endif
;
;   look at the first image in the stack and select the star of interest
;
image0 = readfits(files[0], header)
image0 = reset_high_low_pixels(image0)
image0 = fix_bad_pixels(image0)
filename = FILE_BASENAME(files[0])
if(ilist eq 1) then filename = filename + ' ' + infile_list
image_size = SIZE(image0)
ncol = image_size[1]
nrow = image_size[2]
;
defocus_image_display, filename, image0
defocus_select_star, image0, xcenter, ycenter
wdelete, 10
wdelete, 11
;
;   do a first pass to determine the largest aperture radius indicated by
;   the stars themselves in their images
;
;   at the same time, collect their individual centers, sky values, and HWP angles
;
centers = intarr(nfiles, 2)
sky = fltarr(nfiles,2)
HWP_angle = fltarr(nfiles)
;
rad_max = -1
rad_min = 1000
rad_mean = 0.
;
defocus_find_sky, image0, xcenter, ycenter, sky_mean, sky_sigma
;
;   prepare for cross correlation
;
;     - clear top and bottom few rows
;     - trim image to all zeros under some sigma cut
;
sub_image0 = image0 - sky_mean
;
;sub_image0[*,0:9] = sky_mean
;sub_image0[*,1016:1025] = sky_mean
ind_0 = WHERE(image0 lt 3.0*sky_sigma+sky_mean)
sub_image0[ind_0] = 0.0
;
for ifile = 0, nfiles-1 do begin
    image = readfits(files[ifile],header)
    image = reset_high_low_pixels(image)
    image = fix_bad_pixels(image)
    ;
    HWP_angle[ifile] = get_HWP_angle(header)
    ;
    x_start = xcenter
    y_start = ycenter
    ;
    print,"Centering Star in Frame = ",FILE_BASENAME(files[ifile])
    ;-----------------------------------------------------------------
    ;   correlate against first image to get offsets
    ;
    defocus_find_sky, image, xcenter, ycenter, sky_mean, sky_sigma
    ;
    sub_image = image - sky_mean
    ;sub_image[*,0:9] = sky_mean     ; trim bottom
    ;sub_image[*,1016:1025] = sky_mean   ; and top rows
    ;
    ind = WHERE(image lt 3.0*sky_sigma+sky_mean)
    sub_image[ind] = 0.0
    ;
    ;   correlate this frame with the reference one
    ;
    dpc_correlate, sub_image0, sub_image, x_offset, y_offset
    ;
    print, "Offsets X, Y = ",x_offset, y_offset
    ;---------------------------------------------------------------
    x_start = x_start - x_offset
    y_start = y_start - y_offset
    ;
    ;   then refine star coordinates with centroiding
    ;
    nsigma = 8.
    extent = 18
    defocus_center_locate, image, x_start, y_start, nsigma, extent, x_end, y_end
    ;
    ;   update sky values
    ;
    defocus_find_sky, image, x_end, y_end, sky_mean, sky_sigma
    ;
    ;   and measure the radius needed to collect some large fraction of the star's light
    ;   the actual test looks at the fractional flux increase with a step of 1 pixel in
    ;   radius. If the fractional increase is too small, the radius is then fixed.
    ;
    percentage = 2.5       ; % - roughly equivalent to 5-7 sigma cutoff
    defocus_star_radius, image, x_end, y_end, percentage, radius
    ;
    rad_max = rad_max > radius
    rad_min = rad_min < radius
    rad_mean = rad_mean + radius
    centers[ifile,0] = x_end
    centers[ifile,1] = y_end
    sky[ifile,0] = sky_mean
    sky[ifile,1] = sky_sigma
    ;
endfor
;
;stop
;
rad_mean = float(rad_mean) / float(nfiles)
;
;   report largest aperture radius
;
print, "Smallest, Mean, Largest Radii = ",rad_min, rad_mean, rad_max
;
;   tests - if range of radii more than 65%, have problem
;
if(nfiles gt 2) then begin
    if(abs(rad_max - rad_min)/rad_mean gt 0.65) then begin
        print,"ERROR - Min, Max Radius range too large for these data ....stopping"
        stop
    endif
endif
;
;   compute range of aperture radii to use (one pixel wide steps)
;
rad = (INDGEN(6)-1)+rad_max     ; from 1 pixels below max to 4 pixels beyond max
;
;   open output file and prepare for writing
;
openw, lu_out, file_out, /GET_LUN
;
;   for each image
;
for ifile = 0, nfiles-1 do begin
    ;
    ;  read the image again
    ;
    image = readfits(files[ifile],header)
    image = reset_high_low_pixels(image)
    image = fix_bad_pixels(image)
    ;
    ix0 = 0 > (centers[ifile,0] - rad[5])
    ix1 = (ncol-1) < (centers[ifile,0] + rad[5])
    iy0 = 0 > (centers[ifile,1] - rad[5])
    iy1 = (nrow-1) < (centers[ifile,1] + rad[5])
    nbad=0
    ;
    ;   determination of data range
    ;
    ;   for largest aperture, get min, max data values
    ;
    data_max = MAX(image[ix0:ix1,iy0:iy1])
    data_min = MIN(image[ix0:ix1,iy0:iy1])
    ;
    ;  for each aperture size
    ;
    flux = dblarr(6)
    npix_ap = lonarr(6)
    s_flux = dblarr(6,5)
    mag = dblarr(6)
    s_mag = dblarr(6)
    ;
    for iap = 0, 5 do begin
       xx_center = centers[ifile,0]
       yy_center = centers[ifile,1]
       ;
       ;    for a 7x7 grid of x and y center offsets, jog the aperture
       ;    to look for the peak integrated flux
       ;
       flux_max = -1.0d0
       pixel_max = 0L
       ;
       xbest = 0.0
       ybest = 0.0
       ;
       for xjog = 0, 6 do begin
        for yjog = 0, 6 do begin
            xx = xx_center + xjog - 3
            yy = yy_center + yjog - 3
            ;
            ; range to examine
            ;
            ix0 = 0 > (xx - rad[iap])
            ix1 = (ncol-1) < (xx + rad[iap])
            iy0 = 0 > (yy - rad[iap])
            iy1 = (nrow-1) < (yy + rad[iap])
            ;
            ; clear the accumulators
            ;
            total_flux = 0.0d0
            pixel_count = 0L
            ;
            ; go look at each pixel in the range
            ;
            for ix = ix0, ix1 do begin
               for iy = iy0, iy1 do begin
                   ;
                   ;  compute radius of center of this pixel
                   ;
                   pix_radius = sqrt((ix-centers[ifile,0])^2+(iy-centers[ifile,1])^2)
                   if(pix_radius le rad[iap]) then begin
                       if(image[ix,iy] gt -1.0e5 and image[ix,iy] lt 1.0e5) then begin
                           total_flux = total_flux + image[ix,iy]
                           pixel_count++
                       endif
                   endif
               endfor
            endfor
            ;
            if(total_flux gt flux_max) then begin
               flux_max = total_flux
               pixel_max = pixel_count
               xbest = xx
               ybest = yy
            endif
        endfor
       endfor
       ; ------------------------------------------------------------------------------
       ; maximum found, save results
       ;
       flux[iap] = flux_max
       npix_ap[iap] = pixel_max
       ;
       ;   next, compute the noise sources [in electrons]
       ;
       ;   1. poisson from counts in aperture
       ;
       s_flux[iap,1] = sqrt(flux[iap] * conversion_gain)
       ;
       ;   2. read noise
       ;
       s_flux[iap,2] = sqrt(npix_ap[iap]) * read_noise
       ;
       ;   4. bad pixel contribution
       ;
       s_flux[iap,4] = sqrt(nbad*conversion_gain) * 0.5 * (data_max - data_min)
       ;
       ;   3. sky contribution - from the mean sky value, not variance
       ;
       s_flux[iap,3] = sqrt(float(npix_ap[iap])*conversion_gain*sky[ifile,0])
       ;
       ;   total uncertainty (in ADUs)
       ;
       s_flux[iap,0] = sqrt(s_flux[iap,1]^2 + s_flux[iap,2]^2 + s_flux[iap,3]^2 + s_flux[iap,4]^2) / conversion_gain
       ;
       ;   now, subtract sky contribution to flux in aperture
       ;
       flux[iap] = flux[iap] - sky[ifile,0] * npix_ap[iap]
       ;
       ;   convert to magnitudes and mag uncertainty
       ;
       mag[iap] = -2.5 * alog10(flux[iap])
       s_mag[iap] = 1.086 * s_flux[iap,0] / flux[iap]
       ;
       ;    return constituent noises to ADUs
       ;
       s_flux[iap,1:4] = s_flux[iap,1:4]/conversion_gain
       ;
       ;   print
       ;
       ;             filename, hwp angle, rad, mag, s_mag, flux, total uncert, nbad, poisson, read, sky, bad
       ;
       ;printf, lu_out, format='(a,1h;,f6.1,1h;,f9.5,1h;,f8.5,1h;,f13.2,1h;,f9.2,1h;,i3,1h;,f9.2,1h;,f9.2,1h;,f9.2,1h;,f9,2)',files[ifile],$
       ;           HWP_angle[ifile],mag[iap],s_mag[iap],flux[iap],s_flux[iap,0],nbad,s_flux[iap,1],s_flux[iap,2],$
       ;              s_flux[iap,3],s_flux[iap,4]
       printf, lu_out, format='(a,1h;,f6.2,1h;,2(f6.1,1h;),f6.1,1h;,f9.5,1h;,f8.5,1h;,f13.2,1h;,f9.2,1h;,i3,4(1h;,f9.1))',files[ifile],$
                  HWP_angle[ifile],xbest,ybest,rad[iap],mag[iap],s_mag[iap],flux[iap],s_flux[iap,0],nbad,s_flux[iap,1],s_flux[iap,2],$
                     s_flux[iap,3],s_flux[iap,4]
       ;
       ;   next aperture
       ;
       ;       if the final x, y coordinates were very different than the starting ones, scream and moan
       ;
       if(abs(xbest - centers[ifile,0]) gt 10 or abs(ybest - centers[ifile,1]) gt 10) then begin
            print, "ERROR!!! - Photometry center very different from Original Center for file = ",FILE_BASENAME(files[ifile])
            print, "XBest, YBest = ",xbest,ybest," X,Y Centers = ",centers[ifile,0],centers[ifile,1]
            printf, lu_out, "Previous line has large delta-x or delta-y!!!"
       endif
    endfor
    ;
    ;   next image
    ;
endfor
;
;   close output summary file
;
FREE_LUN, lu_out
;
print, "Defocus Star Photometry is Done!"
end