#
# This TCL script reads a list of information about a set of 
#   TASS Mark IV images, and figures out where the camera
#   was pointing.  It then extracts a subset of stars from
#   a reference catalog, picking only those stars which
#   fall within the boundaries of some exposure.
#
# It writes a file to disk, which looks just like the reference
#   catalog, but contains only stars in the observed regions.
#
# MWR 10/24/2000
#    
# Modified "create_limit_lists" so that it uses the standard routines
#            select_input_list, single_image_info, get_image_value
#   instead of opening the "make_list.out" file directly.
#   This takes longer, but gives us more flexibility.  We can now
#   select only images with type "object".
# MWR 8/3/2001
#
# Modified "create_limit_lists" so that it only includes images
#   which have the "include" flag set to 1
# MWR 2/8/2003
#
# We now support two reference catalogs, one for astrometry and a second
#   for photometry.  We do exactly the same thing to each one: extract
#   a subset for local use into the output directory.  
#   This will allow us to cover star-poor regions of the sky, in which
#   we NEED to include stars with poor photometry in order to find
#   a good astrometric solution.
# MWR 4/15/2003
#
# Remove the "imagelist_name" param from the "refcat.param" file;
#   we can use the "get_list_name" function, which is a safer
#   and more portable way to get the information.
# MWR 4/24/2003


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

source param.tcl


###########################################################################
# Read information from two big catalogs:
#           one for astrometry, one for photometry.
#   Select only stars which are
#   with the boundaries of images taken on a particular night.
#   Spit out a subset of the the big catalogs into two small
#   "reference catalogs" we'll use, one for astrometry and 
#   one for photometry.
#
#
# RETURNS
#    0           if all goes well
#    1           if there's a problem
#    
proc refcat { param_file } {
  global debug

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

  # get information from the param file
  if { [file exists $param_file] != 1 } {
    puts stderr "refcat: can't open param file $param_file"
    return 1
  }
  set refcat_dir(astrom) [get_param $param_file "astrom_refcat_dir"]
  set refcat_name(astrom) [get_param $param_file "astrom_refcat_name"]
  set refcat_dir(photom) [get_param $param_file "photom_refcat_dir"]
  set refcat_name(photom) [get_param $param_file "photom_refcat_name"]
  set ref_ra_col(astrom) [get_param $param_file "astrom_ref_ra_col"]
  set ref_dec_col(astrom) [get_param $param_file "astrom_ref_dec_col"]
  set ref_ra_col(photom) [get_param $param_file "photom_ref_ra_col"]
  set ref_dec_col(photom) [get_param $param_file "photom_ref_dec_col"]
  set imagelist_ra_col [get_param $param_file "imagelist_ra_col"]
  set imagelist_dec_col [get_param $param_file "imagelist_dec_col"]
  set fov [get_param $param_file "fov"]
  set output_dir [get_param $param_file "output_dir"]
  set output_name(astrom) [get_param $param_file "astrom_mini_catalog"]
  set output_name(photom) [get_param $param_file "photom_mini_catalog"]

  # get the name of the "make_list.out" file via a special function
  set imagelist_name [get_list_name]

  # we increase the field of view by a bit to take into account slop in
  #   the pointing
  set fov [expr $fov*1.25]

   
  # do exactly the same thing for both catalogs ....
  foreach reftype { astrom photom } {

    # open the reference catalog
    set file $refcat_dir($reftype)/$refcat_name($reftype)
    if { [file exists $file] != 1 } {
      puts stderr "refcat: input file $file doesn't exist -- quitting"
      return 1
    }
    set ref_fileid [open $file "r"]
    if { [catch { set ref_fileid } ] == 1 } {
      puts stderr "refcat: error opening input file $file -- quitting"
      return 1
    }
  
    # open the list of images
    set file $output_dir/$imagelist_name
    if { [file exists $file] != 1 } {
      puts stderr "refcat: input file $file doesn't exist -- quitting"
      return 1
    }
    set list_fileid [open $file "r"]
    if { [catch { set list_fileid } ] == 1 } {
      puts stderr "refcat: error opening input file $file -- quitting"
      return 1
    }
  
    # open the output file, into which we're going to put the stars
    #   which appear within the image boundaries
    set file $output_dir/$output_name($reftype)
    set output_fileid [open $file "w"]
    if { [catch { set output_fileid } ] == 1 } {
      puts stderr "refcat: error opening output file $file -- quitting"
      return 1
    }
  
    # create a set of lists containing the limits of the observed fields
    if { [create_limit_lists $list_fileid $imagelist_ra_col $imagelist_dec_col \
                $fov min_ra1 max_ra1 min_ra2 max_ra2 min_dec max_dec] != 0 } {
      puts stderr "refcat: create_limit_lists returns with error -- quitting"
      return 1
    }
  
    # now walk through the reference catalog and select only lines containing
    #   stars which fall within a field
    if { [make_subset $ref_fileid $output_fileid \
                 $ref_ra_col($reftype) $ref_dec_col($reftype) \
                 $min_ra1 $max_ra1 $min_ra2 $max_ra2 $min_dec $max_dec] != 0 } {
      puts stderr "refcat: make_subset returns with error -- quitting"
      return 1
    }
  
  
    # close all the files
    close $ref_fileid
    close $list_fileid
    close $output_fileid

  # end of loop over the reference catalog type (astrom or photom)
  }  
 
  
  return 0
}



##############################################################################
# Given a fileid for the list of observed fields, the columns within
#   that file of RA and Dec values, and the size of the
#   field of view in degrees, figure out the limits of RA and Dec
#   for each observed field.  Create a set of lists:
#
#        min_ra1             (low) min RA of each field  (decimal degrees)
#        max_ra1             (low) max RA of each field  (decimal degrees)
#        min_ra2             (high) min RA of each field  (decimal degrees)
#        max_ra2             (high) max RA of each field  (decimal degrees)
#        min_dec             min Dec of each field  (decimal degrees)
#        max_dec             max Dec of each field  (decimal degrees)
#
# If a field does not wrap around the 0-360 degree boundary in RA,
#   then we set only the "min_ra1" and "max_ra1" values to the real
#   limits; we set the "min_ra2" and "max_ra2" values to bogus values
#   in that case.  But if an observed field does wrap around the
#   0-360 boundary, we set both sets of RA limits.
#
# RETURNS:
#    0              if all goes well
#    1              if not
#
proc create_limit_lists { fileid ra_col dec_col fov \
                                   min_ra1 max_ra1 min_ra2 max_ra2 \
                                   min_dec max_dec } {
  global debug 
  upvar $min_ra1 minr1 $max_ra1 maxr1 $min_ra2 minr2 $max_ra2 maxr2
  upvar $min_dec mind  $max_dec maxd

  set DEG2RAD [expr 3.14159/180.0]

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

  set minr1 {}
  set maxr1 {}
  set minr2 {}
  set maxr2 {}
  set mind  {}
  set maxd  {}

  # get the list of all object fields, for which we'll need reference stars
  set param_list [list [list "type" "object"] [list "include" "1"] ]
  set image_list [select_input_list $param_list]
  if { [llength $image_list] == 0 } {
    puts stderr "create_limits_list: no object fields?"
	 return 1
  }

  # we walk through the list of fields, adding the limits of each one 
  #   to the growing TCL lists "minr1", "maxr1", etc.
  foreach img $image_list {
    if { $debug > 0 } {
	   puts "  next image is $img"
	 }
	 if { [single_image_info $img info_list] != 0 } {
	   puts stderr "create_limits_list: single_image_info fails for $img"
		return 1
	 }

    # get central RA and Dec
    set ra_cent [get_image_value $info_list "ra"]
	 set dec_cent [get_image_value $info_list "dec"]
    if { $debug > 0 } {
	   puts "create_limit_lists: ra_cent $ra_cent  dec_cent $dec_cent"
	 }

    # we calculate the number of degrees of RA covered by this field
    set cosdec [expr cos($dec_cent*$DEG2RAD)]
    set ra_half [expr ($fov/2.0)/$cosdec]
    

    # now, we calculate the limits in RA -- must watch for wraparound
    #   for convenience, place the RA limits in variables
    #            inr1, axr1              min, max RA on low end
    #            inr2, axr2              min, max RA on high end (if needed)
    #
    set minr [expr $ra_cent - $ra_half]
    set maxr [expr $ra_cent + $ra_half]
    if { $minr < 0 } {
      set inr1 0
      set axr1 $maxr
      set inr2 [expr 360.0+$minr]
      set axr2 360
    } elseif { $maxr > 360 } {
      set inr1 0
      set axr1 [expr $maxr-360.0]
      set inr2 $minr
      set axr2 360
    } else {
      set inr1 $minr
      set axr1 $maxr
      # there is no need for a second set of RA values, so we just
      #   pick some values which will never include a real star's RA
      set inr2 400
      set axr2 400
    }

    # now, we calculate the limits in Dec.
    #   for convenience, place the Dec
    #            ind, axd                min, max Dec
    set dec_half [expr ($fov/2.0)]
    set ind [expr $dec_cent - $dec_half]
    set axd [expr $dec_cent + $dec_half]

    # now, we append these limits to the lists
    lappend minr1 $inr1
    lappend maxr1 $axr1
    lappend minr2 $inr2
    lappend maxr2 $axr2
    lappend mind  $ind
    lappend maxd  $axd

  }

  return 0
}



##############################################################################
# Given fileids for input reference catalog and output file,
#   and given a set of TCL lists containing field limits:
#
#        min_ra1             (low) min RA of each field  (decimal degrees)
#        max_ra1             (low) max RA of each field  (decimal degrees)
#        min_ra2             (high) min RA of each field  (decimal degrees)
#        max_ra2             (high) max RA of each field  (decimal degrees)
#        min_dec             min Dec of each field  (decimal degrees)
#        max_dec             max Dec of each field  (decimal degrees)
#
# we walk through the reference catalog.  For each star, we check to see
#   if it falls within at least one field.  
#
#         If yes, we copy the line containing that star to the output
#         If no,  we skip star
#
# RETURNS:
#    0              if all goes well
#    1              if not
#
proc make_subset { ref_fileid output_fileid  ra_col dec_col \
                     min_ra1 max_ra1 min_ra2 max_ra2 min_dec max_dec } {
  global debug 

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

  set nfields [llength $min_ra1]

  if { $debug > 2 } {
    set field 0
    while { $field < $nfields } {
      puts [format "make_subset: RA %6.1f %6.1f (%6.1f %6.1f) Dec %6.1f %6.1f " \
                             [lindex $min_ra1 $field] \
                             [lindex $max_ra1 $field] \
                             [lindex $min_ra2 $field] \
                             [lindex $max_ra2 $field] \
                             [lindex $min_dec $field] \
                             [lindex $max_dec $field]  ]
      incr field
    }
  }


  # we walk through the reference catalog, checking one star at a time
  set count 0
  set subcount 0
  while { [gets $ref_fileid line] != -1 } {

    incr count
    if { $debug > 1 } { 
      if { [expr $count%1000] == 0 } {
        puts "count $count"
      }
    }

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

    # sanity check
    if { ([llength $line] < $ra_col) || ([llength $line] < $dec_col) } {
      puts stderr "make_subset: skipping bad line ..$line.."
      continue
    }

    set ra [lindex $line $ra_col]
    set dec [lindex $line $dec_col]

    set field 0
    while { $field < $nfields } {

      if { (($ra >= [lindex $min_ra1 $field]) && \
            ($ra <= [lindex $max_ra1 $field]))      ||    \
           (($ra >= [lindex $min_ra2 $field]) && \
	    ($ra <= [lindex $max_ra2 $field])) } {
        if { ($dec >= [lindex $min_dec $field]) &&  \
	     ($dec <= [lindex $max_dec $field]) } {
          puts $output_fileid $line
          incr subcount
	  break
        }
      }

      incr field
    }

  }

  if { $debug > 1 } {
    puts "make_subset: select $subcount of $count stars"
  }

  return 0
}


namespace eval Refcat {
  namespace export init_astrom_refcat get_astrom_refcat_name \
                   init_photom_refcat get_photom_refcat_name
  variable Astrom_Refcat_Name {}
  variable Photom_Refcat_Name {}

#######################################################################
# Set the name of the mini-astrometric catalog we'll create from a big catalog.
#   The mini-catalog contains only those stars which appear in 
#   the area covered by a night's observations.
#
proc init_astrom_refcat { mini_catalog } {
  global debug
  variable Astrom_Refcat_Name 

  if { $debug > 0 } {
    puts "entering init_astrom_refcat "
  }
  set Astrom_Refcat_Name $mini_catalog

  return 0
}

######################################################################
# Return the name of the mini-astrometric catalog.
#
proc get_astrom_refcat_name { } {
  global debug
  variable Astrom_Refcat_Name

  if { $debug > 0 } {
    puts "entering get_astrom_refcat_name"
  }
  return $Astrom_Refcat_Name
}


#######################################################################
# Set the name of the mini-photometric catalog we'll create from a big catalog.
#   The mini-catalog contains only those stars which appear in 
#   the area covered by a night's observations.
#
proc init_photom_refcat { mini_catalog } {
  global debug
  variable Photom_Refcat_Name 

  if { $debug > 0 } {
    puts "entering init_photom_refcat "
  }
  set Photom_Refcat_Name $mini_catalog

  return 0
}

######################################################################
# Return the name of the mini-photometric catalog.
#
proc get_photom_refcat_name { } {
  global debug
  variable Photom_Refcat_Name

  if { $debug > 0 } {
    puts "entering get_photom_refcat_name"
  }
  return $Photom_Refcat_Name
}


  # end of namespace Refcat
}