pro prepare_polar_science_for_coadd, input_file_list, output_directory, NOGUI=nogui
;
;   reads the FITS headers of a group of polarimetry science images and
;   writes output list files with the same HWP angle and the same object name
;
;   DPC     20051102  copied and modified from "prepare_polar_flats_for_coadd.pro"
;
np = N_PARAMS()
if(N_ELEMENTS(NOGUI) eq 0) then NOGUI = 0
;
if(np eq 0) then begin
    if(NOGUI eq 0) then begin
       files_in = dialog_pickfile(Title="Select Polarimetry FITS Files or File List",/MULTIPLE_FILES)
       nfiles = size(files_in)
       nfiles = nfiles[1]
       if(nfiles eq 1) then begin
            input_file = files_in
       endif else input_file = ' '
       dir_out = dialog_pickfile(Title="Name a Directory for the Output List files",/DIRECTORY)
    endif else begin
       files_in = ' '
       input_file = ''
       print, "Enter pathname to list of Polarimetry FITS files:"
       read, input_file
       dir_out = ''
       print, "Enter pathname for a Directory for the Output List files:"
       read, dir_out
    endelse
endif else begin
    input_file = input_file_list
    files_in = ' '
    dir_out = output_directory
endelse
;
;   if an input file, read it
;
if(input_file ne ' ') then begin
    openr, lu_in, input_file, /GET_LUN
    nfiles = 0
    files_in = strarr(40000)
    onefile = ''
    while (~EOF(lu_in)) do begin
        readf, lu_in, format='(a)',onefile
        files_in[nfiles] = STRTRIM(onefile,2)
        nfiles++
    endwhile
    FREE_LUN, lu_in
    files_in = files_in[0:nfiles-1]
endif
;
;   ----- prepare the data structure(s) -----------
;
info = {filename:'', object_name:'', hwp_angle:0L}
summary = REPLICATE(info,nfiles)
;
;   read all the FITS headers, populating the summary data structure
;
;
for ifile = 0, nfiles-1 do begin
    summary[ifile].filename=files_in[ifile]
    ;
    header = headfits(files_in[ifile])
    ;
    ;   look for keywords
    ;
    ;   start by creating object names appropriate for filenames
    ;
    object = FXPAR(header,"OBJECT",COUNT=ocount)
    if(ocount ne 0) then begin
        ;
        ;  rules for object names:
        ;
        ;  1. trim, convert to uppercase
        ;
        object = STRTRIM(STRUPCASE(object),2)
        ;
        ;  2. remove any "=" signs
        ;
        parts = STRSPLIT(object,'=',/EXTRACT,COUNT=ecount)
        if(ecount gt 1) then begin
           object = parts[0]
           for ii = 1, ecount-1 do begin
             object = object + parts[ii]
           endfor
           object = STRTRIM(object,2)
        endif
        ;
        ;  3. convert blanks, semicolons, and commas to underscores
        ;
        parts = STRSPLIT(object,' ',/EXTRACT,COUNT=ecount)
        if(ecount gt 1) then begin
           object = parts[0]
           for ii = 1, ecount-1 do begin
             object = object + '_' + parts[ii]
           endfor
           object = STRTRIM(object,2)
        endif
       ;
       parts = STRSPLIT(object,';',/EXTRACT,COUNT=ecount)
        if(ecount gt 1) then begin
           object = parts[0]
           for ii = 1, ecount-1 do begin
             object = object + '_' + parts[ii]
           endfor
           object = STRTRIM(object,2)
        endif
        ;
        parts = STRSPLIT(object,',',/EXTRACT,COUNT=ecount)
        if(ecount gt 1) then begin
           object = parts[0]
           for ii = 1, ecount-1 do begin
             object = object + '_' + parts[ii]
           endfor
           object = STRTRIM(object,2)
        endif
        ;
        ;  store the resulting object name in the summary structure
        ;
        summary[ifile].object_name = object
    endif else print, "ERROR - No object name for file ",FILE_BASENAME(files_in[ifile])
    ;
    ; ------ next, determine hwp angle
    ;
    comment = FXPAR(header,"COMMENT",COUNT=count)
    ;
    ;   search the comments for "BD1"
    ;
    if(count gt 0) then begin
        ;
        for icount = 0, count-1 do begin
         comment_UC = STRUPCASE(comment[icount])       ; convert to upper case for comparisons
         ;
         if(STRPOS(comment_UC,"BD1") ne -1) then begin
          ;
          ;   found the "BD1" listing, determine the HWP position angle
          ;
          if(STRPOS(comment_UC,"HOME") ne -1) then summary[ifile].hwp_angle = 0 else begin
             ;
             ;  locate "goto" and ". done" to bracket HWP step number
             ;
             goto_location = STRPOS(comment_UC,"GOTO")
             done_location = STRPOS(comment_UC,". DONE")
             num = done_location - goto_location -5     ; number of characters
             p = 0
             if(num gt 0 and num lt 4) then begin
                HWP_step_string = STRMID(comment_UC,goto_location+5,num)
                reads, HWP_step_string, p
                summary[ifile].hwp_angle = p
             endif else begin
                print, "ERROR - couldn't parse HWP step number. Comment line = ",comment[icount]
             endelse
          endelse
         endif
         ;
       endfor
    endif else begin
        ;
        ;  found NO comments -- FITS header error
        ;
        print,"ERROR - Found no FITS COMMENTS lines for file ",files_in[ifile]
    endelse
    ;
    ;  next file
    ;
endfor
;
;   all files read in to summary structure, process
;
;  *** debug first ***
;
for ifile = 0, nfiles-1 do begin
    print, "File ",summary[ifile].filename," Object Name ",summary[ifile].object_name,"  HWP Steps = ",summary[ifile].hwp_angle
endfor
;
;   sort the structure by object name
;
obj_ind = SORT(summary[*].object_name)
;
;   for each unique object name, make a root file name
;
root = strarr(2000)
root[0] = summary[obj_ind[0]].object_name   ; load first name
nobjects = 0
for i = 1, nfiles-1 do begin
    ;
    ;  test this file name for match with current
    ;
    if(summary[obj_ind[i]].object_name ne root[nobjects]) then begin
       nobjects++
       root[nobjects] = summary[obj_ind[i]].object_name
    endif
endfor
nobjects++
;
;   now have array of root file names and number of objects
;
root = root[0:nobjects-1]
;
;   for each object name, collect all the relevant file entries
;
for iobj = 0, nobjects-1 do begin
    ;
    root_ind = WHERE(summary[*].object_name eq root[iobj])
    ;
    ;  move these to another structure
    ;
    nn = N_ELEMENTS(root_ind)
    object = REPLICATE(info,nn)
    object = summary[root_ind]     ; make a clean copy for just the one object
    ;
    ;  now, sort by HWP angle
    ;
    hwp_ind = SORT(object[*].hwp_angle)
    ;
    ;  extract hwp values in order
    ;
    n_hwp = N_ELEMENTS(hwp_ind)
    ;
    old_hwp = -1
    fileopen = 0
    nentry = 0
    for j = 0, n_hwp-1 do begin
       ;
       ; test for change in hwp value
       ;
       if(object[hwp_ind[j]].hwp_angle ne old_hwp) then begin
         ;
         ;    new hwp found, create file name, open file, write first entry
         ;
         if(fileopen eq 1) then begin
          if(nentry gt 0) then print,"  wrote ",nentry," filenames to output file"
          FREE_LUN, lu_obj       ; close any open file
         endif
         old_hwp = object[hwp_ind[j]].hwp_angle
         path = dir_out + STRTRIM(root[iobj]) + '_HWP_' + STRTRIM(STRING(format='(i)',old_hwp),2) + '.dat'
         ;
         print, "opening output file ",path
         ;
         openw, lu_obj, path, /GET_LUN
         fileopen = 1
         printf, lu_obj, object[hwp_ind[j]].filename   ; only item to save is filename
         nentry = 1
         ;
       endif else begin
         ;
         ;    ok, just another entry, add it to the file
         ;
         printf, lu_obj, object[hwp_ind[j]].filename
         nentry++
       endelse
    endfor
    ;
    ;  close last file
    ;
    FREE_LUN, lu_obj
    ;
    ;  next object
    ;
endfor
;
;   all objects done. quit.
;
;
print, "Prepare Polar Science for Coadd is Done!"
end