# This script performs astrometry on a set of ".pht" files based on
#   Mark IV images.  It compares the detected sources to Tycho stars,
#   finds a transformation, and transforms TASS detections to 
#   (RA, Dec) coords.
#
# Removed the "input_dir" parameter, since this step takes as input
#   files produced by earlier steps of the pipeline -- which are
#   in the "output_dir".  So, this routine 
#       - expects to find its input in "output_dir"
#       - AND place its own output in "output_dir"
#   A bit confusing, but that's the way it is.
# Also added new parameters to control the "match" program.
#     MWR  1/21/2001
#
# Now expect that the reference catalog appears in the "output_dir",
#   rather than in a separate directory.
#     MWR  2/8/2001
#      
# Fixed "create_astrom_subset" to handle cases in which the possible
#   reference star coords wrap around RA = 0 = 360.
#     MWR  3/5/2001
#
# Add "match_scale" param to astrom.param file, send it as an arg
#   to the "run_match" routine.
#     MWR  4/7/2001
#   
# Improvements to "make_error_plot" to give each plot a different name,
#   and handle the case better when wip doesn't exist.
#   Thanks for Paul Roe.
#     MWR  6/14/2001
#
# Modified the section which extracts "delta_dec" from TRANS, because
#   new match v0.5 has extra args at end of TRANS.
# Also, insert special check in "run_match" to make sure that we 
#   don't include the "match_transonly" variable in argument list 
#   to match program if it is NULL.  NULL arguments cause the 
#   match program to halt immediately.
#     MWR 4/7/2002
#
# Changed param "astrom_mag" to "astrom_V_mag" and "astrom_I_mag",
#   so that we use the appropriate magnitudes to match up against detected
#   stars in different filters.
#     MWR 5/6/2002
#
# Modify "create_astrom_subset" to account properly for increase in
#   possible RA values when Dec is far from equator
# Also, modify the code which scans the TRANS for "delta_dec" value
#   so that it expects the proper number of elements; the TRANS now
#   has two extra elements, the "sx" and "sy" values suggested by
#   Blakeslee at JHU.  These elements are generated by match v0.6
#   and later.
#     MWR 11/22/2002
#
# Add new proc "replace_header_values", which puts the newly measured
#   RA and DEC values into the FITS header of the image file; specifically,
#   it modifies the FITS header values for
#        CRVAL1, CRVAL2, CRPIX1, CRPIX2
#     MWR 03/15/2003    
#
# Now use "get_astrom_refcat_name" instead of "get_refcat_name", since
#   we now support two different catalogs.
#     MWR 04/18/2003


if { [catch { set debug } ] != 0 } {
  set debug 2
}

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



#############################################################################
# Given a param file with input values, and the name of one reduced
#   image for which we have already done photometry, 
#   
#       - read and parse the values in the param file
#
#       - create a file of Tycho stars which are in the field of the
#            given TASS image
#
#       - create a subset of TASS measurements for matching
#
#       - run the "match" program 
#       
#       
# Returns
#    0              if all goes well
#    1              if there's an error
#
proc astrom { param_file this_image } {
  global debug

  set DEGTORAD 0.017453293

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

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

  set output_dir       [get_param $param_file "output_dir"]
  set temp_file        [get_param $param_file "temp_file"]

  # parameters for matching up
  set ra_offset        [get_param $param_file "ra_offset"]
  set dec_offset       [get_param $param_file "dec_offset"]
  set iter_center      [get_param $param_file "iter_center"]
  set tass_fov         [get_param $param_file "tass_fov"]
  set tass_npix        [get_param $param_file "tass_npix"]
  set axis_rowpos      [get_param $param_file "axis_rowpos"]
  set axis_colpos      [get_param $param_file "axis_colpos"]
  set tass_nstar       [get_param $param_file "tass_nstar"]
  set tass_subset      [get_param $param_file "tass_subset"]
  set tass_x           [get_param $param_file "tass_x"]
  set tass_y           [get_param $param_file "tass_y"]
  set tass_mag         [get_param $param_file "tass_mag"]
  set astrom_nstar     [get_param $param_file "astrom_nstar"]
  set astrom_infield   [get_param $param_file "astrom_infield"]
  set astrom_subset    [get_param $param_file "astrom_subset"]
  set astrom_x         [get_param $param_file "astrom_x"]
  set astrom_y         [get_param $param_file "astrom_y"]
  set astrom_V_mag     [get_param $param_file "astrom_V_mag"]
  set astrom_I_mag     [get_param $param_file "astrom_I_mag"]
  set match_exe        [get_param $param_file "match_exe"]
  set match_trirad     [get_param $param_file "match_trirad"]
  set match_nobj       [get_param $param_file "match_nobj"]
  set match_matchrad   [get_param $param_file "match_matchrad"]
  set match_order      [get_param $param_file "match_order"]
  set match_transonly  [get_param $param_file "match_transonly"]
  set match_recalc     [get_param $param_file "match_recalc"]
  set match_maxiter    [get_param $param_file "match_maxiter"]
  set match_haltsigma  [get_param $param_file "match_haltsigma"]
  set match_scale      [get_param $param_file "match_scale"]
  set apply_match_exe  [get_param $param_file "apply_match_exe"]
  set photom_match_exe [get_param $param_file "photom_match_exe"]
  set putsym_exe       [get_param $param_file "putsym_exe"]

  # we'll need to know the directory in which the various external
  #   executable programs are located.
  #   Prepend this directory to the name of the executables
  set match_dir [get_directory "match_dir"]
  set match_exe $match_dir/$match_exe
  set apply_match_exe $match_dir/$apply_match_exe
  set photom_match_exe $match_dir/$photom_match_exe

  # we also need to know the names of the "putsym" executable is
  #   located, and the name of the executable itself
  set bait_dir [get_directory "bait_dir"]
  set putsym_exe $bait_dir/$putsym_exe


  # we get the name of the mini-catalog with astrometric reference stars
  #   by calling a special function, rather than using a parameter,
  #   because the name is set with a parameter in ANOTHER param file
  set astrom_catalog [get_astrom_refcat_name]

  # get information about the image
  if { [single_image_info $this_image image_info_list] != 0 } {
    puts stderr "astrom: single_image_info fails for image $this_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 "astrom: image $this_image is not 'object', so skip it"
    }
    return 0
  }

  # we need to know the filter so that we can pick the appropriate 
  #   stars from the astrometric reference catalog in "create_astrom_subset"
  set filter [get_image_value $image_info_list "filter"]
  if { $filter == "V" } {
    set astrom_mag $astrom_V_mag
  } elseif { $filter == "I" } {
    set astrom_mag $astrom_I_mag 
  } else {
    if { $debug > 0 } {
      puts "astrom: image $image: no ref mag for filter $filter, so skip it"
    }
    return 0
  }

  set pht_file [format "%s/%s.pht" $output_dir $this_image]
  set ast_file [format "%s/%s.ast" $output_dir $this_image]

  # a sanity check
  if { ($astrom_nstar < $match_nobj) || ($tass_nstar < $match_nobj) } {
    puts stderr "astrom: nstar(s) smaller than match_nobj?"
  }
 
  # figure out the RA and Dec for this image
  if { [get_ra_dec $this_image $ra_offset $dec_offset \
                                      tass_ra tass_dec ] != 0 } {
    puts stderr "astrom: get_ra_dec fails for $this_image"
    return 1
  }


  # create a subset of TASS stars, for matching purposes
  if { [create_tass_subset $pht_file $output_dir $tass_nstar $tass_subset \
                               $axis_rowpos $axis_colpos \
                               $tass_x $tass_y \
                               $tass_npix $tass_fov $temp_file ] != 0 } {
    puts stderr "astrom: create_tass_subset fails on $this_image"
    return 1
  }

  # now, we run a loop one or two times in order to determine the
  #   central RA and Dec of the TASS image precisely.  What we do
  #   is to extract Tycho stars based on the nominal (RA, Dec) from
  #   the FITS header, run a 'match', and use the translational 
  #   portion of the TRANS to improve the image's central (RA, Dec);
  #   we repeat at most "iter_center" times

  set delta_ra 0
  set delta_dec 0
  set iter -1
  while { $iter < $iter_center } {

    # update our estimate of the position of the center of the image
    set tass_ra  [expr $tass_ra  + ($delta_ra/$DEGTORAD)]
    set tass_dec [expr $tass_dec + ($delta_dec/$DEGTORAD)]
    if { $debug > 1 } {
      puts [format "   new RA %9.4f   new Dec %9.4f " $tass_ra $tass_dec]
    } 

    # create the Tycho catalog data files -- both of them
    if { [create_astrom_subset $output_dir $tass_ra $tass_dec $tass_fov \
                              $astrom_catalog \
			      $astrom_x $astrom_y $astrom_mag \
                              $astrom_infield $astrom_subset \
                              $astrom_nstar $temp_file ] != 0 } {
      puts stderr "astrom: create_astrom_subset fails for $this_image"
      return 1
    }


    # run the matching program on the small subsets
    # Note that if the 'match' program fails to find a good match,
    #   we do print an error message, but return 0, not 1;
    #   this permits execution to continue on other fields.
    if { [run_match $output_dir $tass_subset  $tass_x  $tass_y  $tass_mag \
                    $astrom_subset $astrom_x $astrom_y $astrom_mag \
                    $match_exe $match_trirad $match_nobj $match_matchrad \
                    $match_transonly $match_recalc $match_order \
		    $match_maxiter $match_haltsigma $match_scale \
                    trans_structure ] != 0 } {
      puts stderr "astrom: run_match fails on file $this_image"
      return 0
    }

    if { [decode_trans $trans_structure] != 0 } {
      puts stderr "astrom: decode_trans returns with error"
    }

    set str [lindex $trans_structure 1]
    if { [scan $str "a=%f" delta_ra] != 1 } {
      puts stderr "astrom: can't scan for delta_ra from TRANS $trans_structure"
      return 1
    }
    # for the Dec, we need to check different coeffs if this is a linear
    #   or quadratic TRANS 
    #   This code works with match v 0.6, which added two new elements
    #   to the end of the TRANS structure (the "sx" and "sy" values)
    if { [llength $trans_structure] == 12 } {
      set str [lindex $trans_structure 4]
      set ret [scan $str "d=%f" delta_dec]
    } elseif { [llength $trans_structure] == 18 } {
      set str [lindex $trans_structure 7]
      set ret [scan $str "g=%f" delta_dec]
    } elseif { [llength $trans_structure] == 22 } {
      set str [lindex $trans_structure 9]
      set ret [scan $str "i=%f" delta_dec]
    } else {
      set ret 0
    }
    if { $ret != 1 } {
      puts stderr "astrom: can't scan for delta_dec from TRANS $trans_structure"
      return 1
    }

    incr iter
  }   
  # end of loop to determine central RA and Dec precisely

  # update the RA and Dec items for this image in the image list file
  #   Make sure we format the coordinates nicely: 
  #   decimal degrees, reported to 2 decimal places.
  set nice_ra [format "%.2f" $tass_ra]
  set klist [list "ra" $nice_ra ]
  if { [replace_image_value $this_image $klist] != 0 } {
    puts stderr "astrom: can't replace RA value for image $this_image"
    return 1
  }
  set nice_dec [format "%+.2f" $tass_dec]
  set klist [list "dec" $nice_dec]
  if { [replace_image_value $this_image $klist] != 0 } {
    puts stderr "astrom: can't replace Dec value for image $this_image"
    return 1
  }

  # update the RA and Dec values in the FITS header of the image file
  #   Note that
  #       CRVAL1        is new Dec value, in decimal degrees
  #       CRVAL2        is new RA  value, in decimal degrees
  #       CRPIX1        is col=Dec pixel position of center of optical axis
  #                          (i.e. get from "axis_colpos" param)
  #       CRPIX2        is row=RA  pixel position of center of optical axis
  #                          (i.e. get from "axis_rowpos" param)
  set longer_ra [format "%+.4f" $tass_ra]
  set longer_dec [format "%+.4f" $tass_dec]
  if { [replace_header_values $output_dir $this_image $putsym_exe \
                        $longer_ra $longer_dec \
                        $axis_rowpos $axis_colpos] != 0 } {
    puts stderr "astrom: can't replace FITS header values for image $this_image"
    return 1
  }
  

  # apply the match to the entire set of TASS detections
  #
  # first, we create a file with the 'scaled' versions of all TASS
  #   detections (not just the brightest few)
  set number_in_subset -1
  if { [create_tass_subset $pht_file $output_dir \
                               $number_in_subset $tass_subset \
                               $axis_rowpos $axis_colpos \
                               $tass_x $tass_y \
                               $tass_npix $tass_fov $temp_file ] != 0 } {
    puts stderr "astrom: create_tass_subset fails on file $this_image"
    return 1
  }

  # and now we apply the TRANS to these scaled detections
  #   (we put the TRANSformed detections back into the $tass_subset file)
  if { [apply_trans $output_dir $apply_match_exe $tass_subset $match_order  \
                    $trans_structure $tass_ra $tass_dec $ast_file ] != 0 } {
    puts stderr "astrom: apply_trans fails on file $this_image"
    return 1
  }

  # make a plot showing the errors in the astrometry over the field
  if { [make_error_plot $photom_match_exe $astrom_infield $ast_file \
                        $tass_ra $tass_dec $tass_fov ] != 0 } { 
    puts stderr "astrom: make_error_plot fails"
    return 1
  }


      

  # delete some of the temporary files we've created in the output directory
  foreach f [list $tass_subset $astrom_subset $astrom_infield] {
    if { $debug > 1 } {
      puts "checking for temp file $output_dir/$f"
    }
    if { [file exists $output_dir/$f] == 1 } {
#      exec /bin/rm -f $output_dir/$f
    }
  }
  # these files created by the "match" program
  foreach f { matched.mtA matched.mtB matched.unA matched.unB } {
    if { [file exists $output_dir/$f] == 1 } {
#      exec /bin/rm -f $output_dir/$f
    }
  }
  


  return 0
}


#############################################################################
# Create two simple ASCII data files, containing 
# 
#      1. all the stars from the astrometric catalog which fall within 
#              the TASS field of view, with coords in (RA, Dec)
#              
#      2. the brightest 'astrom_nstar' stars in the field of view,
#              with projected coords in the tangent plane (xi, eta)
#
# All lines within the big astrometric catalog must have
#
#           RA  values in the column given by     "astrom_x" arg
#           Dec values in the column given by     "astrom_y" arg
#           mag values in the column given by     "astrom_mag" arg
# 
# Ignore any line with a "mag" of 0.0.  That's not a real value.
#
# We transform the 'astrom_subset' entries in the following way:
# 
#      - project the stars onto a tangent plane, centered at 
#                 (tass_ra, tass_dec)
#      - calculate plate coordinates (xi, eta) to replace (RA, Dec)
#
# We leave all the other information in the astrometric catalog alone, but
# replace the (RA, Dec) values with these new (xi, eta) values
#      
# Give the datafiles the name specified by the "astrom_infield" 
#   and "astrom_subset" arguments.
# 
# Return
#   0              if all goes well
#   1              if there's a problem
#   
proc create_astrom_subset { output_dir tass_ra tass_dec tass_fov \
                            astrom_catalog \
			    astrom_x astrom_y astrom_mag \
                            astrom_infield astrom_subset \
                            astrom_nstar temp_file } {
  global debug

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

  if { $debug > 1 } {
    puts [format "create_astrom_subset: RA %9.4f  Dec %9.4f " \
                  $tass_ra $tass_dec]
    puts [format "create_astrom_subset: refcat is %s/%s " \
                  $output_dir $astrom_catalog]
  }

                 
  if { [file exists $output_dir/$astrom_catalog] != 1 } {
    puts "create_astrom_subset: can't find $output_dir/$astrom_catalog"
    return 1
  }

  set astrom_infield [format "%s/%s" $output_dir $astrom_infield]
  set astrom_subset  [format "%s/%s" $output_dir $astrom_subset]
  set temp_file      [format "%s/%s" $output_dir $temp_file]

  # create the output files
  set infield_fout [open $astrom_infield "w"]
  set subset_fout [open $astrom_subset "w"]
  if { $debug > 1 } {
    puts "create_astrom_subset: creating files $astrom_infield and $astrom_subset"
  }

  # figure out the limits in RA and Dec for the subset
  #   Must allow for increase in range of RA with Dec.
  #   Make simple check for Dec > 89, use fixed width in RA if so
  set DEGTORAD [expr 3.14159/180.0]
  set half [expr $tass_fov/2.0]
  if { $tass_dec > 89.0 } {
    set half_ra [expr $half*10]
  } else {
    set half_ra [expr $half/cos($tass_dec*$DEGTORAD)]
  }
  set min_ra [expr $tass_ra-$half_ra]
  set max_ra [expr $tass_ra+$half_ra]
  set min_dec [expr $tass_dec-$half]
  set max_dec [expr $tass_dec+$half]

  # we use these in case we need to make a second pass to pick out
  #   stars on the other side of the RA = 360 line, in case
  #   the field in question wraps around RA = 0 = 360.
  set number_of_passes 1
  set min_ra2 0
  set max_ra2 0

  # check to see if we need to make a second pass
  if { ($min_ra < 0) || ($max_ra > 360) } {
    puts "create_astrom_subset: warning: RA wraps around zero!"
    set number_of_passes 2
    if { $min_ra < 0 } {
      set min_ra2 [expr $min_ra + 360.0]
      set max_ra2 360.0
    } 
    if { $max_ra > 360.0 } {
      set min_ra2 0.0
      set max_ra2 [expr $max_ra - 360.0]
    }
  }

  # now read through the data file, and pick out the stars inside 
  #   the TASS field of view
  set fid [open $output_dir/$astrom_catalog "r"]
  set count 0
  set pass 0
  
  # we make at most 2 passes:
  #   pass 0 is usually the only one: 
  #        use min_ra, max_ra
  #   pass 1 occurs if there's wraparound in RA
  #        set min_ra, max_ra to the min_ra2, max_ra2 values
  #
  while { $pass < $number_of_passes } {

    if { $debug > 1 } {
      puts "create_astrom_subset: pass $pass, $min_ra $max_ra"
    }

    while { [gets $fid line] != -1 } {

      # skip lines which begin with the "#" comment character
      if { [regexp "^#" $line] == 1 } {
        continue
      }

      set ra [lindex $line $astrom_x]
      set dec [lindex $line $astrom_y]
      set mag [lindex $line $astrom_mag]

      if { $debug > 2 } {
        puts [format "    ra %9.5f  dec %9.5f  mag %6.3f " $ra $dec $mag ]
      }
     
      if { ($ra < $min_ra) || ($ra > $max_ra) } {
        continue
      }
      if { ($dec < $min_dec) || ($dec > $max_dec) } {
        continue
      }
      if { ($mag <= 0) } {
        continue
      }

      # if we get here, this is a good star.
      #   Put it out, as-is, with coords still in (RA, Dec) into the 
      #   astrom_infield file
      puts $infield_fout $line

      #   And transform its coordinates into the plate coords (xi, eta)
      #   for the astrom_subset
      set proj_line [project_tycho $line $tass_ra $tass_dec \
                                        $astrom_x $astrom_y $astrom_mag]
      if { $debug > 2 } {
        puts "      orig line: $line"
        puts "      proj line: $proj_line"
      }
      puts $subset_fout $proj_line
  
      incr count

    }

    # get ready for another pass (if we need to make one)
    incr pass
    set min_ra $min_ra2
    set max_ra $max_ra2
    seek $fid 0 start
  }

  close $infield_fout
  close $subset_fout
  close $fid

  if { $debug > 1 } {
    puts "create_astrom_subset: found $count good stars"
  }

  # now select the brightest $astrom_nstar of these for the subset file
  set sort_col $astrom_mag
  exec sort -n +$sort_col $astrom_subset > $temp_file 
  exec head -$astrom_nstar < $temp_file > $astrom_subset 
  exec rm $temp_file 

  return 0
}



#############################################################################
# Create a simple ASCII data file, containing a subset of the 
#   stars detected by TASS.  We pick only the brightest 'tass_nstar'
#   stars.
#
# usage: create_tass_subset  tass_file output_dir tass_nstar tass_subset  \ 
#                             axis_rowpos axis_colpos \
#                             tass_xcol tass_ycol \ 
#                             central_ra central_dec chip_pix tass_fov \ 
#                             temp_file 
#
# where
#         tass_file          is a list of stars detected in image,
#                               output of "phot",
#                               with pixel-based coords (x, y)
#
#         output_dir         is a directory we use for holding temporary
#                               files, as well as the final output files
#                               
#         tass_nstar         if >= 0, include only the brightest 'tass_nstar'
#                               detections in the output
#                            if   -1, include all detections in output
#                            
#         tass_subset        is name of output file we'll create
#         
#         axis_rowpos        is the row position on the chip of the optical
#                                axis
#                                
#         axis_colpos        is the col position on the chip of the optical
#                                axis
#
#         tass_xcol          is the column in 'tass_file' with 'x' coord
#         
#         tass_ycol          is the column in 'tass_file' with 'y' coord
#         
#         chip_pix           is the rough number of pixels along one
#                                 dimension of the image (i.e. if chip
#                                 is 512x512, then 'chip_pix' is 512)
#
#         tass_fov           is rough size of field of view, in degrees
#         
#         temp_file          is name we may use for temporary data file
#         
#                                 
#
#
# All lines within the TASS subset should look something like this:
# 
#    1    5.86  416.61   2761  21.00   13.349  0.012 
#    2    8.70  656.01   2756  20.07   15.022  0.052 
#    3    7.27 1232.15   2742  23.47   11.551 -0.003 
#
# where  
#       col 1   is the TASS ID
#       col 2   is the row
#       col 3   is the col
#       col 4   is the local sky
#       col 5   is the local sky sigma
#       col 6   is the instrumental mag from first aperture
#       col 7   is the mag uncertainty from first aperture
#
# There may be several apertures ....
#
# Give the datafile the name specified by the "tass_subset" argument.
#
# We rescale the TASS coords in the following manner, so that they
#   match in size (and direction) the projected Tycho plate coords
#   
#                TASS row - (axis_rowpos)
#        x  =   -------------------------   
#                 (chip_pix/(fov*DEGTORAD))
#
#                TASS col - (axis_colpos)
#        y  =   -------------------------   
#                 (chip_pix/(fov*DEGTORAD))
#
# The resulting scaled coordinates are close to radians, with
#   (0, 0) at the center of the optical axis (not necessarily the
#   center of the image).  This makes matching
#   against the projected Tycho coords easy.
#
#                   
# 
# Return
#   0              if all goes well
#   1              if there's a problem
#   
proc create_tass_subset { tass_file output_dir tass_nstar tass_subset \
                              axis_rowpos axis_colpos \
                              tass_xcol tass_ycol chip_pix tass_fov \
                              temp_file } {
  global debug
  set DEGTORAD 0.017453293

  if { $debug > 0 } {
    puts "entering create_tass_subset"
  }
  if { $debug > 1 } {
    puts "  args are $tass_file $output_dir $tass_nstar $tass_subset "
    puts "           $axis_rowpos $axis_colpos "
    puts "           $tass_xcol $tass_ycol $chip_pix $tass_fov"
  } 

  if { [file exists $tass_file] != 1 } {
    puts "create_tass_subset: can't find file $tass_file"
    return 1
  }

  set temp_file $output_dir/$temp_file
  set tass_subset $output_dir/$tass_subset

  # now sort the stars in the TASS detection file, and place the
  #   brightest 'tass_nstar' ones into 'tass_subset'
  #   (should be able to do this in one step, but it causes problems)
  if { $debug > 1 } {
    puts "create_tass_subset: creating file $tass_subset"
  }
  exec sort -n +5 $tass_file > $temp_file
  if { $tass_nstar >= 0 } {
    exec head -$tass_nstar < $temp_file > $tass_subset
  } else {
    exec cp $temp_file $tass_subset
  }
  exec rm $temp_file

  # rescale the 'row' and 'col' values in the TASS subset, so that 
  #   they match more closely the Tycho (RA, Dec) values
  set half_chip [expr $chip_pix/2]
  set scale_factor [expr $chip_pix/($tass_fov*$DEGTORAD)]
 
  if { $debug > 2 } {
    puts "create_tass_subset: half_chip $half_chip   scale_factor $scale_factor"
  }

  set fid [open $tass_subset "r"]
  set fout [open $temp_file "w"]
  while { [gets $fid line] != -1 } {

    set len [llength $line]

    set oldx [lindex $line $tass_xcol]
    set newx [format "%12.7e" [expr (($oldx-$axis_rowpos)/$scale_factor)] ]

    set oldy [lindex $line $tass_ycol]
    set newy [format "%12.7e" [expr (($oldy-$axis_colpos)/$scale_factor)] ]

    set field 0
    set newline {}
    while { $field < $len } {
      if { $field == $tass_xcol } {
        lappend newline $newx
      } elseif { $field == $tass_ycol } {
        lappend newline $newy
      } else {
        lappend newline [lindex $line $field]
      }
      incr field
    }
 
    if { $debug > 2 } {
      puts [format " old line $line \n"]
      puts [format " new line $newline \n"]
    }

    puts $fout $newline
  }
  close $fid
  close $fout

  exec mv $temp_file $tass_subset


  return 0
}




#############################################################################
# Run the 'match' program on the subsets of the TASS and Tycho catalogs.
#
# Place the TRANS structure into the 'trans_structure' argument on output
# 
# Return
#   0              if all goes well
#   1              if there's a problem
#   
proc run_match { output_dir tass_subset  tass_x  tass_y  tass_mag \
                 astrom_subset astrom_x astrom_y astrom_mag \
                 match_exe match_trirad match_nobj match_matchrad \
                 match_transonly match_recalc match_order \
		 match_maxiter match_haltsigma match_scale \
                 trans_structure } {
  upvar $trans_structure trans_struct
  global debug
  
 
  if { $debug > 1 } {
    puts "entering run_match"
  }

  set tass_subset [format "%s/%s" $output_dir $tass_subset]
  set astrom_subset [format "%s/%s" $output_dir $astrom_subset]

  if { [file exists $tass_subset] != 1 } {
    puts "run_match: can't find file $tass_subset"
    return 1
  }
  if { [file exists $astrom_subset] != 1 } {
    puts "run_match: can't find file $astrom_subset"
    return 1
  }

  if { $debug > 1 } {
    puts "run_match: about to run the matching program"
  }

  if { $debug > 1 } {
    puts "exec $match_exe $tass_subset  $tass_x  $tass_y  $tass_mag "
    puts "          $astrom_subset $astrom_x $astrom_y $astrom_mag "
    puts "          $match_order trirad=$match_trirad nobj=$match_nobj "
    puts "          matchrad=$match_matchrad $match_transonly $match_recalc"
    puts "          max_iter=$match_maxiter halt_sigma=$match_haltsigma"
    puts "          scale=$match_scale "
  }

if { 0 } {
exec ${match_exe}.verbose $tass_subset  \
                  $tass_x  $tass_y  $tass_mag \
                  $astrom_subset $astrom_x $astrom_y $astrom_mag  \
                  $match_order trirad=$match_trirad nobj=$match_nobj \
                  matchrad=$match_matchrad $match_transonly \
		  max_iter=$match_maxiter halt_sigma=$match_haltsigma \
                  scale=$match_scale $match_recalc > match.verbose
}


  # we cause an error if we include an empty argument.  the "match_transonly"
  #   parameter will sometimes be left blank, so we have to check for
  #   it and have special cases for the calls to "match".
  #
  # there must be a better way to do this ...
  #   MWR 4/7/2002
  if { $match_transonly == {} } {
    set retcode [catch { set trans_struct [exec $match_exe $tass_subset  \
                  $tass_x  $tass_y  $tass_mag \
                  $astrom_subset $astrom_x $astrom_y $astrom_mag  \
                  $match_order trirad=$match_trirad nobj=$match_nobj \
                  matchrad=$match_matchrad \
		            max_iter=$match_maxiter halt_sigma=$match_haltsigma \
                  scale=$match_scale $match_recalc]  } ] 
  } else {
    set retcode [catch { set trans_struct [exec $match_exe $tass_subset  \
                  $tass_x  $tass_y  $tass_mag \
                  $astrom_subset $astrom_x $astrom_y $astrom_mag  \
                  $match_order trirad=$match_trirad nobj=$match_nobj \
                  matchrad=$match_matchrad $match_transonly \
		            max_iter=$match_maxiter halt_sigma=$match_haltsigma \
                  scale=$match_scale $match_recalc]  } ] 
  }

  if { $retcode != 0 } {
    puts stderr "run_match: exec $match_exe returns with retcode $retcode"
    return 1
  }

  if { $debug > 1 } {
    puts "run_match: retcode is $retcode"
    puts "run_match: trans is $trans_struct"
  }


  return 0
}



#############################################################################
# Given an image name, figure out the RA and Dec, and place them
#   into the output args
#   
# usage: get_ra_dec      image_name   ra_offset dec_offset  ra   dec
#  
# where
# 
#             image_name           is the name of the (corrected) image
#
#             ra_offset            amount to add to the FITS RA value
#                                       in order to calc true RA
#                                       of center of image (degrees)
#             
#             dec_offset           amount to add to the FITS Dec value
#                                       in order to calc true Dec
#                                       of center of image (degrees)
#             
#             ra                   on output, will contain the RA
#                                       of the image (decimal degrees)
#                                       
#             dec                  on output, will contain the Dec 
#                                       of the image (decimal degrees)
#                                       
# Returns
#   0                 if all goes well
#   1                 if there's a problem
#   
proc get_ra_dec { image_name ra_offset dec_offset ra dec } {
  upvar $ra ra_val $dec dec_val
  global debug

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

  if { [single_image_info $image_name pair_list] != 0 } {
    puts stderr "get_ra_dec: single_image_info returns with error"
    return 1
  }

  if { $debug > 1 } {
    puts "  get_ra_dec: pair_list is $pair_list"
  }

  set ra_val  [get_image_value $pair_list "ra"]
  set dec_val [get_image_value $pair_list "dec"]

  set ra_val  [expr $ra_val + $ra_offset]
  set dec_val [expr $dec_val + $dec_offset]

  if { $debug > 1 } {
    puts "  get_ra_dec: ra is $ra_val   dec is $dec_val"
  }

  return 0
}


############################################################################
# Given the name of a file containing all the (scaled) TASS detections
#   for some image, and given a TRANS structure, 
#   
#       1. apply the TRANS structure to all the detections;
#                   this yields (xi, eta) coords in the tangent plane
#                   which is centered on the given central_ra, central_dec
#                   
#       2. deproject the (xi, eta) coords into (RA, Dec), using the
#                   central_ra and central_dec
#
#   Put the resulting, transformed, (RA, Dec) values for all detections 
#   into the file specified by the 'output_file' arg
#   
# usage:   apply_trans       output_dir apply_match_exe \
#                            tass_detections match_order \
#                            trans_structure central_ra central_dec \
#                            output_file
#
# where
#            output_dir              is a directory in which the data files
#                                         all exist (both input file
#                                         "tass_detections" and output file
#                                         "output_file")
#
#            apply_match_exe         is the full path of a Unix program
#                                         which applies TRANS structures
#                                         to datafiles
#                                         
#            tass_detections         is a file containing the (scaled)
#                                         coords of all objects detected
#                                         in a TASS image
#                                         
#            trans_structure         is a string containing the TRANS 
#                                         coeffs, returned by 'match'
#
#            match_order             describes the order of the transformation
#                                         of the plate solution:
#                                              linear
#                                              quadratic
#                                              cubic
#                                         
#            central_ra              is the central RA of the field,
#                                         in decimal degrees
#                                         
#            central_dec             is the central Dec of the field
#                                         in decimal degrees
#
#            output_file             is name of a file into which we place
#                                         the results
# 
# Returns
#    0               if all goes well
#    1               if there's a problem
#
proc apply_trans { output_dir apply_match_exe tass_detections match_order \
                                    trans_structure \
                                    central_ra central_dec output_file } {
  global debug

  if { $debug > 0 } {
    puts "entering apply_trans, tass_detections is $tass_detections"
  }

  # we need to prepend the output directory name to the "tass_detections"
  #   file; but the "output_file" already has the output directory in it
  set tass_detections [format "%s/%s" $output_dir $tass_detections]

  # sanity check
  if { [file exists $apply_match_exe] != 1 } { 
    puts stderr "apply_trans: can't find program $apply_match_exe"
    return 1
  }
  if { [file exists $tass_detections] != 1 } { 
    puts stderr "apply_trans: can't find file $tass_detections"
    return 1
  }

  # first, we decode the TRANS structure, placing the 16 coefficients
  #   into vars "a, b, c, d, e, f ... m, n, o, p".  If the TRANS has
  #   a low order (linear or quadratic), then just read the elements
  #   it does contain, and ignore the rest
  if { [string compare $match_order "linear"] == 0 } {
    if { [scan $trans_structure "%s a=%f b=%f c=%f d=%f e=%f f=%f" \
                           junk a b c d e f ] != 7 } {
      puts stderr "apply_trans: the following linear TRANS is invalid"
      puts stderr "   $trans_structure "
      return 1
    }
  } elseif { [string compare $match_order "quadratic"] == 0 } {
    if { [scan $trans_structure "%s a=%f b=%f c=%f d=%f e=%f f=%f g=%f h=%f i=%f j=%f k=%f l=%f" \
                           junk a b c d e f g h i j k l] != 13 } {
      puts stderr "apply_trans: the following quadratic TRANS is invalid"
      puts stderr "   $trans_structure "
      return 1
    }
  } elseif { [string compare $match_order "cubic"] == 0 } {
    if { [scan $trans_structure "%s a=%f b=%f c=%f d=%f e=%f f=%f g=%f h=%f i=%f j=%f k=%f l=%f m=%f n=%f o=%f p=%f" \
                           junk a b c d e f g h i j k l m n o p ] != 17 } {
      puts stderr "apply_trans: the following TRANS is invalid"
      puts stderr "   $trans_structure "
      return 1
    }
  } else {
    puts stderr "apply_trans: invalid match_order $match_order"
    return 1
  }

  # now, we call the "apply_match" command to do the transformations 
  #   for us
  set cmd ""
  set cmd [format "exec $apply_match_exe $tass_detections 1 2 "]
  set cmd [format "%s %s" $cmd "$central_ra $central_dec "]
  if { [string compare $match_order "linear"] == 0 } {
    set cmd [format "%s %s" $cmd "linear"]
    set cmd [format "%s %s" $cmd "$a $b $c $d $e $f "]
  } elseif { [string compare $match_order "quadratic"] == 0 } {
    set cmd [format "%s %s " $cmd "quadratic"]
    set cmd [format "%s %s " $cmd "$a $b $c $d $e $f $g $h $i $j $k $l "]
  } elseif { [string compare $match_order "cubic"] == 0 } {
    set cmd [format "%s %s " $cmd "cubic"]
    set cmd [format "%s %s " $cmd "$a $b $c $d $e $f $g $h $i $j $k $l $m $n $o $p"]
  } else {
    puts stderr "apply_trans: invalid match_order $match_order"
    return 1
  }
  set cmd [format "%s %s " $cmd "outfile=$output_file "]
  puts "cmd is $cmd"

  if { [catch { eval $cmd } ] != 0 } {
    puts stderr "apply_trans: $apply_match_exe returns with error"
    return 1
  }
  

  return 0
}


############################################################################
# Given the name of a file containing all the (scaled) TASS detections
#   for some image, and given a TRANS structure, 
#   
#       1. apply the TRANS structure to all the detections;
#                   this yields (xi, eta) coords in the tangent plane
#                   which is centered on the given central_ra, central_dec
#                   
#       2. deproject the (xi, eta) coords into (RA, Dec), using the
#                   central_ra and central_dec
#
#   Put the resulting, transformed, (RA, Dec) values for all detections 
#   back into the original file.
#   
# usage:   apply_trans       apply_match_exe tass_detections trans_structure \ 
#                            central_ra  central_dec
#
# where
#            apply_match_exe         is the full path of a Unix program
#                                         which applies TRANS structures
#                                         to datafiles
#                                         
#            tass_detections         is a file containing the (scaled)
#                                         coords of all objects detected
#                                         in a TASS image
#                                         
#            trans_structure         is a string containing the TRANS 
#                                         coeffs, returned by 'match'
#                                         
#            central_ra              is the central RA of the field,
#                                         in decimal degrees
#                                         
#            central_dec             is the central Dec of the field
#                                         in decimal degrees
#
#            temp_file               is name of a file we can use for
#                                         temp storage of results
# 
# Returns
#    0               if all goes well
#    1               if there's a problem
#
proc apply_trans6 { apply_match_exe tass_detections trans_structure \
                                    central_ra central_dec temp_file } {
  global debug

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

  # sanity check
  if { [file exists $apply_match_exe] != 1 } { 
    puts stderr "apply_trans6: can't find program $apply_match_exe"
    return 1
  }
  if { [file exists $tass_detections] != 1 } { 
    puts stderr "apply_trans6: can't find file $tass_detections"
    return 1
  }

  # first, we decode the TRANS structure, placing the 12 coefficients
  #   into vars "a, b, c, d, e, f, g, h, i, j, k, l"
  if { [scan $trans_structure "%s a=%f b=%f c=%f d=%f e=%f f=%f g=%f h=%f i=%f j=%f k=%f l=%f" \
                           junk a b c d e f g h i j k l] != 13 } {
    puts stderr "apply_trans6: the following TRANS is invalid"
    puts stderr "   $trans_structure "
    return 1
  }
puts "trans values are $a $b $c $d $e $f .. $g $h $i $j $k $l"

  # now, we call the "apply_match" command to do the transformations 
  #   for us
  puts "exec $apply_match_exe $tass_detections 1 2 "
  puts "                          $central_ra $central_dec "
  puts "                          $a $b $c $d $e $f $g $h $i $j $k $l"
  puts "                          outfile=$temp_file "
  if { [catch { exec $apply_match_exe $tass_detections 1 2 \
                                  $central_ra $central_dec \
                                  $a $b $c $d $e $f $g $h $i $j $k $l \
                                  outfile=$temp_file } ] != 0 } {
    puts stderr "apply_trans6: $apply_match returns with error"
    return 1
  }

if { 0 } {
  exec mv $temp_file $tass_detections 
}
puts "tass_detections in $tass_detections   output in $temp_file"

  return 0
}


############################################################################
# Given the name of a file containing all the (scaled) TASS detections
#   for some image, and given a TRANS structure, 
#   
#       1. apply the TRANS structure to all the detections;
#                   this yields (xi, eta) coords in the tangent plane
#                   which is centered on the given central_ra, central_dec
#                   
#       2. deproject the (xi, eta) coords into (RA, Dec), using the
#                   central_ra and central_dec
#
#   Put the resulting, transformed, (RA, Dec) values for all detections 
#   back into the original file.
#   
# usage:   apply_trans8      apply_match_exe tass_detections trans_structure \ 
#                            central_ra  central_dec
#
# where
#            apply_match_exe         is the full path of a Unix program
#                                         which applies TRANS structures
#                                         to datafiles
#                                         
#            tass_detections         is a file containing the (scaled)
#                                         coords of all objects detected
#                                         in a TASS image
#                                         
#            trans_structure         is a string containing the TRANS 
#                                         coeffs, returned by 'match'
#                                         
#            central_ra              is the central RA of the field,
#                                         in decimal degrees
#                                         
#            central_dec             is the central Dec of the field
#                                         in decimal degrees
#
#            temp_file               is name of a file we can use for
#                                         temp storage of results
# 
# Returns
#    0               if all goes well
#    1               if there's a problem
#
proc apply_trans8 { apply_match_exe tass_detections trans_structure \
                                    central_ra central_dec temp_file } {
  global debug

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

  # sanity check
  if { [file exists $apply_match_exe] != 1 } { 
    puts stderr "apply_trans8: can't find program $apply_match_exe"
    return 1
  }
  if { [file exists $tass_detections] != 1 } { 
    puts stderr "apply_trans8: can't find file $tass_detections"
    return 1
  }

  # first, we decode the TRANS structure, placing the 16 coefficients
  #   into vars "a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p"
  if { [scan $trans_structure "%s a=%f b=%f c=%f d=%f e=%f f=%f g=%f h=%f i=%f j=%f k=%f l=%f m=%f n=%f o=%f p=%f" \
                           junk a b c d e f g h i j k l m n o p] != 17 } {
    puts stderr "apply_trans8: the following TRANS is invalid"
    puts stderr "   $trans_structure "
    return 1
  }
puts "trans values are $a $b $c $d $e $f $g $h .. $i $j $k $l $m $n $o $p"

  # now, we call the "apply_match" command to do the transformations 
  #   for us
  puts "exec $apply_match_exe $tass_detections 1 2 "
  puts "                   $central_ra $central_dec "
  puts "                   $a $b $c $d $e $f $g $h $i $j $k $l $m $n $o $p"
  puts "                   outfile=$temp_file "
  if { [catch { exec $apply_match_exe $tass_detections 1 2 \
                           $central_ra $central_dec \
                           $a $b $c $d $e $f $g $h $i $j $k $l $m $n $o $p\
                           outfile=$temp_file } ] != 0 } {
    puts stderr "apply_trans8: $apply_match returns with error"
    return 1
  }

if { 0 } {
  exec mv $temp_file $tass_detections 
}
puts "tass_detections in $tass_detections   output in $temp_file"

  return 0
}



############################################################################
# Given a line of information from an astrometric catalog, in which
#
# where  
#       col 'astrom_x'      is the RA (decimal degrees)
#       col 'astrom_y'      is the Dec (decimal degrees)
#       col 'astrom_mag'    is a magnitude we'll use for rough sorting
#
# calculate the "plate coordinates" of each star, on a tangent plane
#   which is centered at the given RA and Dec.
#   
# usage:  project_tycho       tycho_line central_ra central_dec
#                                        astrom_x  astrom_y  astrom_mag
# 
# where
#            tycho_line       is a line of information from Tycho subset
#            
#            central_ra       is central RA of tangent plane (degrees)
#            
#            central_dec      is central Dec of tangent plane (degrees)
#
#            astrom_x         is the column containing RA values in catalog
#
#            astrom_y         is the column containing Dec values in catalog
#
#            astrom_mag       is the column containing mag values in catalog
#            
# Returns
#   a formatted line with projected coords     if all goes well
#   terminates execution                       if something goes wrong
#   
proc project_tycho { tycho_line central_ra central_dec \
                             astrom_x astrom_y astrom_mag } {
  global debug
  set DEGTORAD 0.017453293

  if { $debug > 2 } {
    puts "entering project_tycho"
  }

  # there must be at least 'last_col' columns in each line of catalog
  set last_col $astrom_x
  if { $astrom_y > $last_col } {
    set last_col $astrom_y
  }
  if { $astrom_mag > $last_col } {
    set last_col $astrom_mag
  }

  # sanity checks
  set llen [llength $tycho_line]
  if { $llen < $last_col } {
    puts stderr "project_tycho: bad line $tycho_line"
    puts stderr "quitting now"
    exit 1
  }

  if { $debug > 2 } {
    puts "   input line $tycho_line"
  }

  set ra [lindex $tycho_line $astrom_x]
  set dec [lindex $tycho_line $astrom_y]
  set mag [lindex $tycho_line $astrom_mag]


  # to convert from (RA, Dec) to project plate coords (xi, eta), 
  #    use the complex formulae shown in Olkin et al., 
  #    PASP 108, 202 (1996).
  #    (this is better done in C for speed, I'm sure)
  set cosdec0 [expr cos($DEGTORAD*$central_dec)]
  set sindec0 [expr sin($DEGTORAD*$central_dec)]
  set cosdec  [expr cos($DEGTORAD*$dec)]
  set sindec  [expr sin($DEGTORAD*$dec)]
if { 1 } {
  set sin_delta_ra [expr sin($DEGTORAD*($ra - $central_ra))]
  set cos_delta_ra [expr cos($DEGTORAD*($ra - $central_ra))]
} else {
  set sin_delta_ra [expr sin($DEGTORAD*($central_ra - $ra))]
  set cos_delta_ra [expr cos($DEGTORAD*($central_ra - $ra))]
}

  # calculate xi, the pseudo-RA coordinate 
  set xx [expr $cosdec*$sin_delta_ra]
  set yy [expr $sindec0*$sindec + $cosdec0*$cosdec*$cos_delta_ra]
  set xi [expr $xx/$yy]
        
  # calculate eta, the pseudo-Dec coordinate 
  set xx [expr $cosdec0*$sindec - $sindec0*$cosdec*$cos_delta_ra]
  set eta [expr $xx/$yy]


  # now create the output line
  set ret_line " "
  set word_index 0
  while { $word_index < $llen } {
    set word [lindex $tycho_line $word_index]
    if { $word_index == $astrom_x } {
      set ret_line [format "%s %12.7e " $ret_line $xi ]
    } elseif { $word_index == $astrom_y } {
      set ret_line [format "%s %12.7e " $ret_line $eta ]
    } else {
      set ret_line [format "%s %s " $ret_line $word ]
    }
    incr word_index
  }

  if { $debug > 2 } {
    puts "  output line $ret_line"
  }
  return $ret_line
}


############################################################################
# Given a line of information from a TASS photometry file, like this:
#
#    1   0.002341 -0.001103  349   7.87   15.740  0.039 
#    2  -0.010353  0.022387  344   7.80   16.415  0.068 
#
#
#       col 1   is the TASS ID number in the image
#       col 2   is the eta (pseudo-RA) coord in the projection
#       col 3   is the xi  (pseudo-Dec) coord in the projection
#       col 4   is the local sky value
#       col 5   is the local sky sigma
#       col 6   is the instrumental mag
#       col 7   is the estimated uncertainty in the instrumental mag
#
# calculate the (RA, Dec) coords of each star, by de-projecting the
#   (xi, eta) "plate coordinates" which are centered at the given RA and Dec.
#   
# usage:  deproject_tass      tass_line central_ra central_dec
# 
# where
#            tass_line        is a line of information from TASS photometry
#                                  file, with (xi, eta) coords
#            
#            central_ra       is central RA of tangent plane (degrees)
#            
#            central_dec      is central Dec of tangent plane (degrees)
#            
# Returns
#   a formatted line with de-projected coords  if all goes well
#   terminates execution                       if something goes wrong
#   
proc deproject_tass { tass_line central_ra central_dec } {
  global debug
  set DEGTORAD 0.017453293

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

  # sanity checks
  if { [llength $tass_line] != 7 } {
    puts stderr "deproject_tass: bad line $tass_line"
    puts stderr "quitting now"
    exit 1
  }

  if { $debug > 1 } {
    puts "   input line $tass_line"
  }

  if { [scan $tass_line "%d %f %f %f %f %f %f" \
                             tass_id xi eta sky skysig mag magsig ] != 7 } {
    puts stderr "deproject_tass: bad line $tass_line"
    puts stderr "quitting now"
    exit 1
  }

  # to convert from plate coords (xi, eta) to spherical coords (RA, Dec)
  #    use the complex formulae shown in Olkin et al., 
  #    PASP 108, 202 (1996).
  #    (this is better done in C for speed, I'm sure)
  set cosdec0 [expr cos($DEGTORAD*$central_dec)]
  set sindec0 [expr sin($DEGTORAD*$central_dec)]

  # calculate RA
  set xx $xi
  set yy [expr ($cosdec0 - $eta*$sindec0)]
  set delta_ra [expr atan($xx/$yy)]
  set ra [expr $central_ra - $delta_ra/$DEGTORAD)]

  # calculate Dec 
  set xx [expr ($sindec0 + $eta*$cosdec0)*(cos($delta_ra))]
  set zz [expr ($xx/$yy)]
  set dec [expr atan($zz)/$DEGTORAD]


  # now create the output line
  set ret_line [format "%d %10.5f %10.5f %5d %6.2f %7.3f %7.3f " \
                      $tass_id $ra $dec $sky $skysig $mag $magsig]

  if { $debug > 1 } {
    puts "  output line $ret_line"
  }
  return $ret_line
}



#############################################################################
# Make a plot which shows distortion over the field of view, by matching
#   up detected stars with Tycho stars.
#   
# Uses the "wip" program to make the Postscript plot.  If the "wip" program
#   isn't available, just print a warning message and return immediately
#   with code 0.
#
# Returns
#   0          if all goes well (or if "wip" doesn't exist, 
#                                  and we return without making any plots)
#   1          if there's a real problem
# 
proc make_error_plot { photom_match_exe tycho_file tass_file \
                       tass_ra tass_dec tass_fov } {
  global debug

  set DEGTORAD 0.017453293
  set dist_factor 0.1

  # divide the field of view into a grid with this many cells in 
  #   each direction
  set ngrid  10

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

  # if the "wip" program doesn't exist, then don't both to make
  #   the plotting input file; just print a warning message and return
  #   immediately, with code 0
  if { [catch { exec wip -q -d /dev/null -e exit }] != 0 } {
    puts stderr "make_error_plot: can't execute wip, so make no plots"
    return 0
  }

  # place output of "photom_match" here
  set output_file photom_match.out

  # place data to be used in making plot here
  set plot_file make_error_plot.dat

  # place input script for WIP here
  set wip_file  make_error_plot.wip

  # wip output file
  set wip_ps_file  ${tass_file}.error_plot.ps

  # maximum dist between stars to be considered a match (arcsec)
  set match_radius  10       

  if { $debug > 0 } {
    puts "    exec $photom_match_exe $tass_file 2 3 2000 v=6 "
    puts "         $tycho_file 2 3 2000 v=4  radius=$match_radius "
    puts "         > $output_file  "
  }
  set retcode  [catch { exec $photom_match_exe $tass_file 2 3 2000 v=6 \
                              $tycho_file 2 3 2000 v=4  radius=$match_radius \
                              > $output_file } ] 
  if { $retcode != 0 } {
    puts stderr "make_error_plot fails, retcode $retcode "
    return 1
  }

  # calculate 
  #    a. min, max values of RA and Dec (in decimal degrees)
  #    b. boundaries of the grid cells in RA and Dec
  set plot_radius [expr (($tass_fov*1.1)/2)]
  set min_ra [expr $tass_ra-$plot_radius]
  set max_ra [expr $tass_ra+$plot_radius]
  set min_dec [expr $tass_dec-$plot_radius]
  set max_dec [expr $tass_dec+$plot_radius]

  set grid_elem_size [expr $plot_radius/($ngrid*0.5)]
  
  # we need three arrays for each grid cell, to hold
  #   grid_array_num         number of points in the cell
  #   grid_array_ravec       sum of RA distances in the cell
  #   grid_array_decvec      sum of Dec distances in cell
  #
  # each array contains $ngrid*$ngrid elements, and we calculate an
  #   index like this:
  #   
  #       index = (xpos*$ngrid) + ypos
  #       
  #   where 'xpos' and 'ypos' are the positions within the grid
  set index 0
  set num_cells [expr $ngrid*$ngrid]
  while { $index < $num_cells } {
    set grid_array_num($index) 0
    set grid_array_ravec($index) 0
    set grid_array_decvec($index) 0
    incr index
  }
 

  # now, read the data in the output file and modify it slightly so that
  #   we have the information we need to make a vector plot.
  set fid [open $output_file "r"]
  while { [gets $fid line] != -1 } {

    if { [regexp \#.* $line ] == 1 } {
      continue
    }

    set ra [lindex $line 0]
    set dec [lindex $line 1]
    set dra [lindex $line 2]
    set ddec [lindex $line 3]
    set dist [lindex $line 4]
    set raw_mag [lindex $line 5]
    set cat_mag [lindex $line 6]

    # figure out which grid element this star inhabits
    set grid_x [expr int(floor(($ra-$min_ra)/$grid_elem_size))]
    set grid_y [expr int(floor(($dec-$min_dec)/$grid_elem_size))]
    set grid_index [expr int((($grid_x*$ngrid) + $grid_y))]

    # add these values to the running sums in "grid_array" arrays
    set grid_array_num($grid_index) [expr $grid_array_num($grid_index)+1]
    set grid_array_ravec($grid_index) \
                   [expr $grid_array_ravec($grid_index)+$dra]
    set grid_array_decvec($grid_index) \
                   [expr $grid_array_decvec($grid_index)+$ddec]

    if { $debug > 3 } {
    puts [format "grid %4d now has %9.4f %9.4f %6.2f " $grid_index \
                  $grid_array_num($grid_index) \
                  $grid_array_ravec($grid_index) \
                  $grid_array_decvec($grid_index) ]
    }

  }
  close $fid

  # now, finally, we can write the average values for distance and 
  #   angle in each grid cell into the plotting file
  set plot_fid [open $plot_file "w"]
  set grid_index 0
  while { $grid_index < $num_cells } {

    set n [expr int($grid_array_num($grid_index))]
    set x_pos [expr int(floor($grid_index/$ngrid))]
    set y_pos [expr int(floor($grid_index - $x_pos*$ngrid))]
    if { $debug > 3 } {
      puts [format "   n %6s  x %6s y %6s " $n $x_pos $y_pos]
      puts [format "   n %6d  x %6d y %6d " $n $x_pos $y_pos]
    }
 
    set ra  [expr $min_ra  + (($x_pos + 0.5)*$grid_elem_size)]
    set dec [expr $min_dec + (($y_pos + 0.5)*$grid_elem_size)]
   

    if { $n == 0 } {
      puts $plot_fid [format "%9.4f %9.4f %9.4f %9.4f " $ra $dec 0 0]

    } else {

      set avg_ra [expr $grid_array_ravec($grid_index)/$n]
      set avg_dec [expr $grid_array_decvec($grid_index)/$n]
      if { $debug > 3 } {
        puts [format "    avg_ra %9.4f  avg_dec %9.4f " $avg_ra $avg_dec]
      }

      # calculate the average distance between TASS and Tycho positions,
      #   scale the "dist" from arcseconds to degrees, for plotting
      set dist [expr sqrt($avg_ra*$avg_ra + $avg_dec*$avg_dec)]
      set scaled_dist [expr $dist*$dist_factor]

      # calculate the angle of the error, in degrees, ccw from pos x
      set angle [expr atan2($avg_dec, $avg_ra)]
      set angle [expr $angle/$DEGTORAD]
      puts $plot_fid [format "%9.4f %9.4f %9.4f %9.4f " \
                                 $ra $dec $scaled_dist $angle]
      if { $debug > 3 } {
        puts [format "        dist %9.4f  scaled %9.4f  angle %6.1f" \
                            $dist $scaled_dist $angle]
      }
    }
    
    incr grid_index
  }
  close $plot_fid

  # now we make an input script for the WIP plotting package
  set fid [open $wip_file "w"]
  puts $fid "device $wip_ps_file/PS"
  puts $fid "data $plot_file"
  puts $fid "env $min_ra $max_ra $min_dec $max_dec 1"
  puts $fid "xc 1"
  puts $fid "yc 2"
  puts $fid "pc 3"
  puts $fid "ecol 4"
  puts $fid "expand 0.4"
  puts $fid "vector 70 0.7"
  close $fid

  # call WIP to make the postscript file
  exec wip < $wip_file

  # get rid of temp files we've created
  if { 0 } {
    exec rm $wip_file
    exec rm $plot_file
  }


  return 0
}


#############################################################################
# Given a linear TRANS structure (list of its elements, actually),
#   decode it into scale and rotation angle.
#   
# Returns
#   0          if all goes well
#   1          if there's a problem
#
proc decode_trans { trans } {
  global debug

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

  set DEGTORAD 0.017453293

  if { [llength $trans] != 7 } {
    puts stderr "decode_trans: invalid linear TRANS "
    return 1
  }
  
  if { [scan [lindex $trans 1] "a=%f" a] != 1 } {
    puts stderr [format "decode_trans: invalid A value %s" [lindex $trans 1] ]
    return 1
  }
  if { [scan [lindex $trans 2] "b=%f" b] != 1 } {
    puts stderr [format "decode_trans: invalid B value %s" [lindex $trans 2] ]
    return 1
  }
  if { [scan [lindex $trans 3] "c=%f" c] != 1 } {
    puts stderr [format "decode_trans: invalid C value %s" [lindex $trans 3] ]
    return 1
  }
  if { [scan [lindex $trans 4] "d=%f" d] != 1 } {
    puts stderr [format "decode_trans: invalid D value %s" [lindex $trans 4] ]
    return 1
  }
  if { [scan [lindex $trans 5] "e=%f" e] != 1 } {
    puts stderr [format "decode_trans: invalid E value %s" [lindex $trans 5] ]
    return 1
  }
  if { [scan [lindex $trans 6] "f=%f" f] != 1 } {
    puts stderr [format "decode_trans: invalid F value %s" [lindex $trans 6] ]
    return 1
  }
 
  # now it's safe to calculate the needed factors; we'll take the average
  #   of the moduli of b and f, and c and e, to get slightly higher
  #   accuracy than if we just took b and f.
  set avg_bf [expr (abs($b)+abs($f))/2.0]
  set avg_ce [expr (abs($c)+abs($e))/2.0]
  set scal [expr sqrt(($avg_bf*$avg_bf)+($avg_ce*$avg_ce))]

  set pa [expr atan2($avg_bf, $avg_ce)/$DEGTORAD]

  puts [format "    scale %9.5f    pa %7.2f " $scal $pa]
 
  return 0
}



###########################################################################
# PROCEDURE: replace_header_values
#
# DESCRIPTION:
#   Update the RA and Dec values in the FITS header of the image file
#   Note that we will try to change
#
#         CRVAL1        is new Dec value, in decimal degrees
#         CRVAL2        is new RA  value, in decimal degrees
#         CRPIX1        is col=Dec pixel position of center of optical axis
#                            (i.e. get from "axis_colpos" param)
#         CRPIX2        is row=RA  pixel position of center of optical axis
#                            (i.e. get from "axis_rowpos" param)
#
# RETURNS
#     0          if all goes well
#     1          if an error occurs

proc replace_header_values { output_dir this_image putsym_exe \
                                       nice_ra nice_dec \
                                       axis_rowpos axis_colpos } {
  
  global debug

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

  # update the RA and Dec values
  set cmd "$putsym_exe file=$output_dir/$this_image "
  set cmd "$cmd CRVAL1=$nice_dec CRVAL2=$nice_ra "
  if { $debug > 1 } {
    puts "cmd is ..$cmd.. "
  }
  if { [ catch { eval exec $cmd } ] != 0 } {
    puts stderr "replace_header_values: putsym CRVAL1 CRVAL2 fails for file $this_image"
    return 1
  }

  # update the CRPIX1 and CRPIX2 values
  set cmd "$putsym_exe file=$output_dir/$this_image "
  set cmd "$cmd CRPIX1=$axis_colpos CRPIX2=$axis_rowpos "
  if { $debug > 1 } {
    puts "cmd is ..$cmd.. "
  }
  if { [ catch { eval exec $cmd } ] != 0 } {
    puts stderr "replace_header_values: putsym CRPIX1 CRPIX2 fails for file $this_image"
    return 1
  }


  return 0
}