# This script finds stars in a set of cleaned Mark IV frames
#   and places lists of stars (and photometry) into files
#   with extensions .coo and .pht
#   
#      clean image            H3R1586.695   
#      output of 'stars'      H3R1586.695.coo
#      output of 'phot'       H3R1586.695.pht
#      
#
# Added "gain" and "readnoise" parameters for photometry.
#   MWR 1/21/2001
#
# Added "xvista_dir" to all calls to XVista programs
#   --  MWR 2/9/2001
# 
# Removed "input_dir" from params.  We look in the "output_dir"
#   for cleaned images, and we also place our lists of stars there.
#   --  MWR 3/5/2001
#
# Fixed bug in "run_phot" -- was passing gain to phot program 
#   with arg "gain=xxx", should have been using arg "adu=xxx".
#   --  MWR 3/7/2001
#
# Improved "exec" call to wip in "make_fwhm_plot", and file names
#   created by that routine.
#   Thanks to Paul Roe.
#   --  MWR 6/14/2001
#
# Modified code to 
#      - count the stars in the image after finding them
#      - set the "nstar" field in the image list
#      - compare with acceptable limit on min number of stars
#      - only run the "phot" step if the image had acceptable number of stars
#   --  MWR 2/8/2003
#
#
# Removed "sanity check" that sky be greater than one; since we now
#   may subtract the background before finding stars, this isn't
#   going to be true all the time.
#   --  MWR 3/15/2003   (actually checked in 4/1/2003)
#
# Changed most of the parameters in the "stars.param" file to be
#   keyed lists, so that we can have separate values for the V and I
#   cameras.  For example, if the two cameras have different readnoise,
#   we'll now properly calculate its contribution correctly for each
#   camera.
#   --  MWR 4/11/2003
#
# When we call the "sky" program to calculate sky values, we now
#   use the "get_skyparams" function to read the parameters associated
#   with "sky", rather than reading those parameters from our own
#   "stars.param" file.  This allows us to share the single set
#   of values defined in the "sky.param" file.
# Also, the "sky" program may print a WARNING message to stderr.
#   We now use "catch" to prevent that warning message from
#   terminating the pipeline execution.  If the warning is 
#   caught, we mark the image as bad and return with error code
#   (but at least other frames can still be reduced).
#   --  MWR 5/15/2003
#
# Added support for "empscatter" parameter in "stars.param" file.
#   If empscatter  ==  1, then we pass "empscatter" option to phot program
#                            (noise determined from scatter in sky values)
#   If empscatter  ==  0, then we don't 
#                            (noise determined from sky value, not scatter)
#   --  MWR 3/26/2005


if { [catch { set debug } ] == 1 } {
  set debug 2
}

# contains source for the "get_param" proc
source param.tcl



#############################################################################
# Given a param file name, and a cleaned image, 
#   parse the param file for a few control values and then
#   run the XVista programs
#   
#       stars                to find stars in the images
#       phot                 to calculate their brightness
#       
# We assume files are in directories
#
#       output_dir           the corrected image files
#       output_dir           the .coo and .pht files these routines create
#
# Returns
#    0              if all goes well
#    1              if there's an error
#
proc stars { param_file clean_image } {
  global debug

  if { [file exists $param_file] != 1 } {
    puts stderr "stars: can't find param file $param_file"
    return 1
  }

  # the XVista programs can be found in a special directory
  set xvista_dir [get_directory "xvista_dir"]
  if { $xvista_dir == {} } {
    puts stderr "stars: get_directory fails for xvista_dir"
    return 1
  }

  set output_dir       [get_param $param_file "output_dir"]

  # parameters for calculating the sky value and sigma
  if { [get_skyparams skybin sky_min_accept sky_max_accept] != 0 } {
    puts stderr "stars: get_skyparams fails"
    return 1
  }

  # parameters for the "stars" program
  set minsig           [get_param $param_file "minsig"]
  set minfwhm          [get_param $param_file "minfwhm"]
  set maxfwhm          [get_param $param_file "maxfwhm"]
  set minround         [get_param $param_file "minround"]
  set maxround         [get_param $param_file "maxround"]
  set minsharp         [get_param $param_file "minsharp"]
  set maxsharp         [get_param $param_file "maxsharp"]
  set maxstars         [get_param $param_file "maxstars"]
  set movemax          [get_param $param_file "movemax"]
 
  # parameters for the "check_stars" procedure
  set min_stars        [get_param $param_file "min_stars"]

  # parameters for the "phot" program
  set apers            [get_param $param_file "apers"]
  set sky_inner        [get_param $param_file "sky_inner"]
  set sky_outer        [get_param $param_file "sky_outer"]
  set readnoise        [get_param $param_file "readnoise"]
  set gain             [get_param $param_file "gain"]
  set saturate         [get_param $param_file "saturate"]
  set empscatter       [get_param $param_file "empscatter"]


  if { $debug > 1 } {
    puts "stars: about to work on file $clean_image"
  }

  # get information about the image
  if { [single_image_info $clean_image image_info_list] != 0 } {
    puts stderr "stars: single_image_info fails for image $clean_image"
   return 1
  }
  # if this is not an "object" image, don't process it!  Just return
  #   with code 0, but do print a warning message
  set type [get_image_value $image_info_list "type"]
  if { [string compare $type "object"] != 0 } {
    if { $debug > 0 } {
      puts "stars: image $clean_image is not 'object', so skip it"
    }
    return 0
  }
  # we need to know the filter so that we can apply the appropriate mask
  #   file in 'run_phot', below
  set filter [get_image_value $image_info_list "filter"]

    # do everything necessary to find stars
    if { [run_stars $output_dir $xvista_dir \
                       $clean_image $filter $skybin \
                       $sky_min_accept $sky_max_accept \
                       $minsig \
                       $minfwhm $maxfwhm \
                       $minround $maxround \
                       $minsharp $maxsharp $maxstars $movemax] != 0 } {
      puts stderr "stars: run_stars fails on file $clean_image"
      puts stderr "quitting now "
      return 1
    }

    # Check to see if the number of stars found in the frame is above
    #   the minimum acceptable level.  If not, set the "include"
    #   flag for this image to zero.
    if { [check_stars $clean_image $min_stars] != 0 } {
      puts stderr "stars: check_stars returns with error"
      return 1
    }

    # make a diagnostic plot showing variation in FWHM across frame
    if { [make_fwhm_plot $output_dir $clean_image] != 0 } {
      puts stderr "stars: make_fwhm_plot returns with error"
      return 1
    }

    # If the frame is still acceptable, 
    #   do everything necessary to measure properties of the stars
    if { [image_include_check $clean_image] == 1 } {
      if { [run_phot $output_dir $xvista_dir \
                         $clean_image $apers \
                         $sky_inner $sky_outer \
		         $readnoise $gain $saturate $empscatter $filter] != 0 } {
        puts stderr "stars: run_phot fails on file $clean_image"
        puts stderr "quitting now "
        return 1
      }
    }

  return 0
}


#############################################################################
# Prepare to run the 'stars' program on the given image, then do so
# 
#   - find the sky value and sky sigma
#   - create an appropriate output file name, by appending ".coo"
#          to the image name
#   - run the 'stars' command
#   - count the number of stars found, and set the "nstar" value
#          for the image in the image list
#
# We assume
#
#          output_dir            holds the images
#          output_dir            will hold the .coo file we'll produce
#          xvista_dir            holds the XVista programs
#
#   
# Return 
#    0               if all goes well
#    1               if there's a problem
#    
proc run_stars { output_dir xvista_dir \
                                 file filter skybin \
                                 sky_min_accept sky_max_accept \
                                 minsig minfwhm maxfwhm \
                                 minround maxround \
                                 minsharp maxsharp maxstars movemax } {
  global debug

  if { $debug > 1 } {
    puts "entering run_stars"
  }

  set infile $output_dir/$file
  if { [file exists $infile] != 1 } {
    puts stderr "run_stars: can't find file $infile"
    return 1
  }

  # find the sky value and sky sigma value
  #   if the "sky" program calculates suspicious results, it will
  #      print a warning to stderr, and we will not get the results
  #      put into "str".  In that case, we have no choice but to
  #      mark the image as "bad" and return with an error code.
  #
  if { [catch { set str [exec $xvista_dir/sky $infile bin=$skybin \
                         min=$sky_min_accept max=$sky_max_accept] } ] != 0 } {
    puts stderr "run_stars: exec'ing sky caused exception for $infile "
    puts stderr "run_stars: marking $infile as 'bad' "
    if { [designate_bad $infile] != 0 } {
      puts stderr "run_stars: designate_bad fails on $infile"
      return 1
    }
    return 1
  }
  set ll [lindex $str 0]
  set ll [split $ll "="]
  set skyval [lindex $ll 1]
  set ll [lindex $str 1]
  set ll [split $ll "="]
  set skysig [lindex $ll 1]

  if { $debug > 1 } {
    puts "run_stars: file $file has sky $skyval  skysig $skysig"
  }


  # some sanity checks
  if { $skysig < 1 } {
    puts stderr "run_stars: skysig value for $file is $skysig, less than one"
    puts stderr "run_stars: quitting now "
    return 1
  }

  if { $debug > 1 } {
    puts "run_stars: file $file has sky $skyval  skysig $skysig"
  }

  # now get the elements of parameters which depend on the filter
  set minsig [get_image_value $minsig $filter]
  set minfwhm [get_image_value $minfwhm $filter]
  set maxfwhm [get_image_value $maxfwhm $filter]
  set minround [get_image_value $minround $filter]
  set maxround [get_image_value $maxround $filter]
  set minsharp [get_image_value $minsharp $filter]
  set maxsharp [get_image_value $maxsharp $filter]
  set movemax [get_image_value $movemax $filter]
  set maxstars [get_image_value $maxstars $filter]
  if { $debug > 2 } {
    puts "run_stars: minsig $minsig  minfwhm $minfwhm  maxfwhm $maxfwhm "
    puts "run_stars: minround $minround  maxround $maxround  "
    puts "run_stars: minsharp $minsharp  maxsharp $maxsharp  "
    puts "run_stars: movemax $movemax  maxstars $maxstars  "
  }


  # create the output file name
  set outfile $output_dir/${file}.coo

  # now run the 'stars' program
  set retcode [ catch { exec $xvista_dir/stars $infile outfile=$outfile \
                   sky=$skyval skysig=$skysig minsig=$minsig \
                   minfwhm=$minfwhm maxfwhm=$maxfwhm \
                   minround=$minround maxround=$maxround \
                   minsharp=$minsharp maxsharp=$maxsharp \
                   maxstars=$maxstars movemax=$movemax quiet } ]

  if { $retcode != 0 } {
    puts stderr "run_stars: exec stars returns with error code $retcode"
    return 1
  }

  # Figure out the average fwhm of stars we found, and set that value in
  #   the list file for this image
  if { [find_avg_fwhm $outfile avg_fwhm] != 0 } {
    puts stderr "run_stars: find_avg_fwhm fails for file $outfile "
    puts stderr "   use value 0.0 for fwhm, continue pipeline "
  }
  set klist [list "fwhm" $avg_fwhm]
  if { $debug > 0 } {
    puts "run_stars: file $infile has avg FWHM $avg_fwhm"
  }
  if { [replace_image_value $file $klist] != 0 } {
    puts stderr "run_stars: replace_image_value for FWHM fails for $file"
    return 1
  }

  # How many stars did we find?  Just count the lines in the .coo file
  if { [catch { set ret [exec wc $outfile] } ] != 0 } {
    puts stderr "run_stars: exec wc to count stars in $outfile fails"
    return 1
  }
  set nstar [lindex $ret 0]
  set klist [list "nstar" $nstar]
  if { $debug > 0 } {
    puts "run_stars: file $infile has $nstar stars"
  }
  if { [replace_image_value $file $klist] != 0 } {
    puts stderr "run_stars: replace_image_value fails for $file"
    return 1
  }


  return 0
}


#############################################################################
# PROCEDURE: check_stars
#
# Given an image name and a keyed list of acceptable minimum numbers
#   of stars, like this:
#
#         { V 1000 } { I 800 }
#
#   determine the filter of this image, and check to see if it has
#   more than the minimum number of stars.  If not, set its include
#   flag to zero.
#
# RETURNS
#    0        if all goes well
#    1        if not
#
proc check_stars { filename min_star_list } {
  global debug

  if { $debug > 0 } {
    puts "entering check_stars"
  }


  # get information about this image
  if { [single_image_info $filename image_info_list] != 0 } {
    puts "check_stars: single_image_info fails for image $filename"
    return 1
  }

  # pick out the filter
  set filter [get_image_value $image_info_list "filter"]
  if { $debug > 3 } {
    puts "check_stars: filter for $filename is $filter"
  }

  # get the minimum acceptable value for number of stars for this filter
  set this_min_nstar [get_image_value $min_star_list $filter]
  if { $debug > 3 } {
    puts "check_stars: min value in filter $filter is $this_min_nstar"
  }

  # Now, compare to the actual number of stars found
  set nstar [get_image_value $image_info_list "nstar"]
  if { $nstar < $this_min_nstar } {
    if { $debug > 0 } {
      puts "check_stars: image $filename has nstar $nstar < $this_min_nstar"
      puts "check_stars:    designating $filename as bad"
    }
    if { [designate_bad $filename] != 0 } {
      puts stderr "check_stars: designate_bad fails on $filename"
      return 1
    }
  }

  return 0
}



#############################################################################
# Prepare to run the 'phot' program on the given image, then do so
# 
#   - create an appropriate output file name, by appending ".pht"
#          to the image name
#   - run the 'phot' command
#   
# We assume that files are located as follows:
#
#       output_dir         contains the corrected images
#       output_dir         contains the .coo file (already created)
#                          will contain the .pht file we produce here
#       xvista_dir            holds the XVista programs
#
#
# Return 
#    0               if all goes well
#    1               if there's a problem
#    
proc run_phot { output_dir xvista_dir \
                           file apers sky_inner sky_outer \
                           readnoise gain saturate empscatter filter } {
  global debug

  if { $debug > 1 } {
    puts "entering run_phot"
  }

  set imagefile $output_dir/$file
  if { [file exists $imagefile] != 1 } {
    puts stderr "run_phot: can't find file $imagefile"
    return 1
  }

  # check to see if a mask file exists; if it does, we'll use it
  set maskfile [master_mask_name $filter]
  if { [file exists $maskfile] == 0 } {
    puts "run_phot: can't find mask file $maskfile, so will use no mask"
    set use_mask 0
  } else {
    set use_mask 1
  }


  # some sanity checks
  if { $sky_inner < 1 } {
    puts stderr "run_phot: sky_inner value for $file is $sky_inner, less than one"
    puts stderr "run_phot: quitting now "
    return 1
  }
  if { $sky_outer < 1 } {
    puts stderr "run_phot: sky_outer value for $file is $sky_outer, less than one"
    puts stderr "run_phot: quitting now "
    return 1
  }
  if { $readnoise < 0 } {
    puts stderr "run_phot: readnoise value for $file is $readnoise, less than zero"
    puts stderr "run_phot: quitting now"
    return 1
  }
  if { $gain <= 0 } {
    puts stderr "run_phot: gain value for $file is $gain, must be > zero"
    puts stderr "run_phot: quitting now"
    return 1
  }

  # now we deal with the parameters, which are lists keyed by filter.
  # We need to choose the proper value for the current filter
  #   before we build the command line.
  #
  # The "apers" parameter is especially tricky, since it may have
  #   a list of elements for each filter....  We build a string
  #   which contains all the apertures.
  set aper_str {}
  set this_aper_list [get_image_value $apers $filter]
  foreach ap $this_aper_list {
    lappend aper_str aper=$ap
  }
  if { $debug > 1 } {
    puts "run_phot: file $file has apers string $aper_str "
  }

  # get the other elements of parameters which depend on the filter
  set sky_inner [get_image_value $sky_inner $filter]
  set sky_outer [get_image_value $sky_outer $filter]
  set readnoise [get_image_value $readnoise $filter]
  set gain [get_image_value $gain $filter]
  set saturate [get_image_value $saturate $filter]
  if { $debug > 2 } {
    puts "run_phot: sky_inner $sky_inner  sky_outer $sky_outer "
    puts "run_phot: readnoise $readnoise  gain $gain   saturate $saturate "
  }

  # the "empscatter" argument is a little different from the others
  #    if it equals 1,  then we add "empscatter" keyword to the command line
  #                         of the "phot" executable.  This means that we 
  #                         will calc sky noise empirically, looking at
  #                         the scatter from mean/median in sky annulus
  #    if it equals 0,  then we do not add the "empscatter" keyword
  #                         to the command line.  In this case, the program
  #                         will calc sky noise based on number of electrons
  #                         it believes came from sky + readnoise + dark.
  if { $empscatter == 1 } {
    set empscatter_val "empscatter"
  } else {
    set empscatter_val " "
  }


  # create the input and output file name
  set infile  $output_dir/${file}.coo
  set outfile $output_dir/${file}.pht

  # now run the 'phot' program
  if { $use_mask == 0 } {
    if { $debug > 1 } {
      puts "run_phot: exec $xvista_dir/phot $imagefile \
                     infile=$infile outfile=$outfile \
                     $aper_str sky_inner=$sky_inner sky_outer=$sky_outer \
   		     readnoise=$readnoise adu=$gain \
                     $empscatter_val \
                     saturate=$saturate "
    }
    set retcode [ catch { exec $xvista_dir/phot $imagefile \
                     infile=$infile outfile=$outfile \
                     $aper_str sky_inner=$sky_inner sky_outer=$sky_outer \
  		     readnoise=$readnoise adu=$gain \
                     $empscatter_val \
                     saturate=$saturate } ]
  } else {
    if { $debug > 1 } {
      puts "run_phot: exec $xvista_dir/phot $imagefile \
                     infile=$infile outfile=$outfile \
		     maskfile=$maskfile \
                     $aper_str sky_inner=$sky_inner sky_outer=$sky_outer \
   		     readnoise=$readnoise adu=$gain \
                     $empscatter_val \
                     saturate=$saturate "
    }
    set retcode [ catch { exec $xvista_dir/phot $imagefile \
                     infile=$infile outfile=$outfile \
		     maskfile=$maskfile \
                     $aper_str sky_inner=$sky_inner sky_outer=$sky_outer \
  		     readnoise=$readnoise adu=$gain \
                     $empscatter_val \
                     saturate=$saturate } ]
  }

  if { $retcode != 0 } {
    puts stderr "run_phot: exec phot returns with error code $retcode"
    return 1
  }

  return 0
}



##############################################################################
# Given the name of an image, look at the ".coo" file associated with
#   it.  Make a plot which indicates variation in FWHM across the field
# 
# To make plots, we use the 'wip' program available from 
#   University of Maryland:
#
#      http://bima.astro.umd.edu/wip/
#
# If the program doesn't exist, we print an error message, but continue
#   program execution
#   
# usage:   make_fwhm_plot      output_dir image_name
# 
# Returns
#   0         if all goes well 
#                    (or if WIP doesn't exist)
#   1         if there's an error
#
proc make_fwhm_plot { output_dir image } {
  global debug
  
  if { $debug > 1 } {
    puts "entering make_fwhm_plot"
  }

  # we use the program "wip" to make plots; if the program isn't installed
  #   on the system, don't try to make a plot (but don't terminate the
  #   entire script execution, either)
  if { [catch { exec wip -q -d /dev/null -e exit }] != 0 } {
    puts stderr "make_fwhm_plot: can't execute wip, so no plots will be made"
    return 0
  }

  set coo_name $output_dir/${image}.coo
  set output_file $output_dir/${image}.fwhm_plot.ps
  set wip_file $output_dir/make_fwhm_plot.wip
  set ngrid 10
  set num_cells [expr $ngrid*$ngrid]
  set ccd_size 2030
  set cell_width [expr int($ccd_size/$ngrid)]

  if { [file exists $coo_name] != 1 } {
    puts stderr "make_fwhm_plot: file $coo_name doesn't exist"
    return 1
  }

  # initialize two arrays, 
  # 
  #    sum_array        sum of FWHM in each grid cell
  #    num_array        number of stars in each grid cell
  set cell 0
  while { $cell < $num_cells } {
    set sum_array($cell) 0
    set num_array($cell) 0
    incr cell
  }
  
  # run through the stars, calculating the mean FWHM inside
  #   each square on a grid within the image
  set fid [open $coo_name "r"]
  while { [gets $fid line] != -1 } {
    if { [llength $line] != 7 } {
      continue
    }
    set xc [lindex $line 1]
    set yc [lindex $line 2]
    set fwhm [lindex $line 4]

    set x_cell [expr int($xc/$cell_width)]
    set y_cell [expr int($yc/$cell_width)]

    set cell_index [expr ($x_cell*$ngrid)+$y_cell]
   
    set num_array($cell_index) [expr $num_array($cell_index)+1]
    set sum_array($cell_index) [expr $sum_array($cell_index)+$fwhm]

    if { $debug > 3 } {
      puts [format "   %6.0f %6.0f %5.2f  %4d %4d  %5d   %5d %7.2f " \
            $xc $yc $fwhm $x_cell $y_cell $cell_index \
            $num_array($cell_index) $sum_array($cell_index)]
    }
  }

  set fout [open $wip_file "w"]
  puts $fout "device $output_file/CPS"
  puts $fout "env 0 1 0 1 1"
  puts $fout "move 0.1 1.1"
  puts $fout "label FWHM:  orange 1.75-2.00    dark blue 2.75-3.00"
  puts $fout "fill 1"
  set x_cell 0 
  while { $x_cell < $ngrid } { 

    # the 'x' values denote rows, which run from 0 at top to 1 at bottom
    set yplot_max [expr 1.0-$x_cell*(1.0/$ngrid)]
    set yplot_min [expr $yplot_max-(1.0/$ngrid)]
   
    set y_cell 0
    while { $y_cell < $ngrid } {

      # the 'y' values denote cols, which run from 0 at left to 1 at right
      set xplot_min [expr $y_cell*(1.0/$ngrid)]
      set xplot_max [expr $xplot_min+(1.0/$ngrid)]

      set cell_index [expr ($x_cell*$ngrid)+$y_cell]
      if { $num_array($cell_index) < 1 } {
        set mean 0
      } else {
        set mean [expr $sum_array($cell_index)/$num_array($cell_index)]
      }

      # select the proper color here 
      set color [select_color $mean]
      puts $fout "color $color"
      puts $fout [format "rect  %6.3f %6.3f %6.3f %6.3f " \
                       $xplot_min $xplot_max $yplot_min $yplot_max]

      if { $debug > 3 } {
        puts [format "   x %4d  y %4d   num %5d  mean %6.2f %3d " \
                       $x_cell $y_cell $num_array($cell_index) $mean $color]
      }

      incr y_cell
    }

    incr x_cell
  }
  puts $fout "color 1"
  close $fout

  exec wip < $wip_file 
  exec rm $wip_file

  puts "look for FWHM plot in $output_file"

  return 0
}


#############################################################################
# Given the mean FWHM (in pixels) of stars in some portion of the field,
#   return a color for filling in that portion of the field.
#   
# See help for the "wip" plotting program.  Colors are:
# 
#        n    Color                     n     Color
#            --------------------------------------------------------------
#        0    Black (background)        8     Orange (Red + Yellow)
#        1    White (foreground)        9     Green + Yellow
#        2    Red                      10     Green + Cyan
#        3    Green                    11     Blue + Cyan
#        4    Blue                     12     Blue + Magenta
#        5    Cyan (Blue + Green)      13     Red + Magenta
#        6    Magenta (Red + Blue)     14     Dark Gray
#        7    Yellow (Red + Green)     15     Light Gray
#                                                                    
# Returns
#    the appropriate color
#    
proc select_color { fwhm } {
  global debug
 
  if { $fwhm < 1.0 } {
    return 0
  }
  if { $fwhm < 1.25 } {
    return 1
  }
  if { $fwhm < 1.50 } {
    return 13
  }
  if { $fwhm < 1.75 } {
    return 2
  }
  if { $fwhm < 2.00 } {
    return 8
  }
  if { $fwhm < 2.25 } {
    return 7
  }
  if { $fwhm < 2.50 } { 
    return 3
  }
  if { $fwhm < 2.75 } {
    return 5
  }
  if { $fwhm < 3.0 } {
    return 4
  }
  if { $fwhm < 3.25 } {
    return 12
  }
  if { $fwhm < 3.50 } {
    return 15
  }
  return 14

}




###########################################################################
# PROCEDURE: find_avg_fwhm
#
# DESCRIPTION: given a .coo file, calculate the average FWHM value
#              for stars.  Place the avg FWHM into the second argument.
#
# RETURNS:
#   0           if all goes well
#   1           if we can't calculate avg FWHM
#
proc find_avg_fwhm { coo_file output_avg_fwhm } {
  global debug
  upvar $output_avg_fwhm avg_fwhm

  if { [file exists $coo_file] != 1 } {
    puts stderr "find_avg_fwhm: file $coo_file doesn't exist"
    return 1
  }

  set sum_fwhm 0
  set num_stars 0
  set avg_fwhm 0
  
  # run through the stars, calculating the mean FWHM inside
  #   each square on a grid within the image
  set fid [open $coo_file "r"]
  while { [gets $fid line] != -1 } {
    if { [llength $line] != 7 } {
      continue
    }
    set fwhm [lindex $line 4]

    set sum_fwhm [expr $sum_fwhm+$fwhm]
    incr num_stars

    if { $debug > 3 } {
      puts [format "   this fwhm %5.2f  sum now %9.2f  num %6d  " \
            $fwhm $sum_fwhm $num_stars]
    }
  }
  close $fid

  if { $num_stars > 0 } {
    set avg_fwhm [expr $sum_fwhm/$num_stars]
  }
  # make the average FWHM look nice; much easier for human to read
  set avg_fwhm [format "%.2f" $avg_fwhm]

  return 0
}