# This is a translation of the "refcat.tcl" script into Perl.
#   The TCL script was taking over half an hour to create reference catalogs
#   when given a really big reference catalog (1,500,000 stars), so
#   I decided to make a Perl script which can be used as an alternative.
#   It runs about four times faster on my Linux desktop.
#
#       usage:   refcat.pl  param_file
#
# MWR  4/18/2003
#  
#
# 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


use strict;

# some global variables
my($debug);
my($nimages);
my(@min_ra1, @max_ra1, @min_ra2, @max_ra2, @min_dec, @max_dec);


$debug = 0;

# we are going to need a set of (possibly) big arrays to contain
#    information on the limits in RA and Dec for all the images
#    taken during this run.  These arrays will hold that information.
#    We set a maximum size which OUGHT to be reasonable.
$nimages = 0;
$min_ra1[0] = 0;
$max_ra1[0] = 0;
$min_ra2[0] = 0;
$max_ra1[0] = 0;
$min_dec[0] = 0;
$max_dec[0] = 0;

###########################################################################
# 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
#    


# vars needed in the main loop
my($param_file);
my($line); my(@words);
my($reftype);
my($ref_file); my($image_file); my($output_file);

# parameters we read from the param file
my(%refcat_dir); my(%refcat_name); 
my(%ref_ra_col); my(%ref_dec_col);
my($imagelist_name); my($imagelist_ra_col); my($imagelist_dec_col);
my($fov);
my($output_dir); my(%output_name);



if ($#ARGV != 0) {
  printf stderr "usage: refcat.pl param_file\n";
  exit(1);
}
$param_file = $ARGV[0];
if ($debug > 0) {
  printf "param_file is $param_file \n";
}

##   # get information from the param file
  open(PARAM_FILE, "$param_file") || 
            die("refcat.pl: can't open file param_file $param_file");
  while (<PARAM_FILE>) {
    $line = " " . $_;
    chomp($line);
    @words = split(/\s+/, $line);
    if ($words[1] eq "astrom_refcat_dir") {
      $refcat_dir{"astrom"} = $words[2];
    }
    if ($words[1] eq "astrom_refcat_name") {
      $refcat_name{"astrom"} = $words[2];
    }
    if ($words[1] eq "photom_refcat_dir") {
      $refcat_dir{"photom"} = $words[2];
    }
    if ($words[1] eq "photom_refcat_name") {
      $refcat_name{"photom"} = $words[2];
    }
    if ($words[1] eq "astrom_ref_ra_col") {
      $ref_ra_col{"astrom"} = $words[2];
    }
    if ($words[1] eq "astrom_ref_dec_col") {
      $ref_dec_col{"astrom"} = $words[2];
    }
    if ($words[1] eq "photom_ref_ra_col") {
      $ref_ra_col{"photom"} = $words[2];
    }
    if ($words[1] eq "photom_ref_dec_col") {
      $ref_dec_col{"photom"} = $words[2];
    }
    if ($words[1] eq "imagelist_name") {
      $imagelist_name = $words[2];
    }
    if ($words[1] eq "imagelist_ra_col") {
      $imagelist_ra_col = $words[2];
    }
    if ($words[1] eq "imagelist_dec_col") {
      $imagelist_dec_col = $words[2];
    }
    if ($words[1] eq "fov") {
      $fov = $words[2];
    }
    if ($words[1] eq "output_dir") {
      $output_dir = $words[2];
    }
    if ($words[1] eq "astrom_mini_catalog") {
      $output_name{"astrom"} = $words[2];
    }
    if ($words[1] eq "photom_mini_catalog") {
      $output_name{"photom"} = $words[2];
    }
  }
  close(PARAM_FILE);

##   # we increase the field of view by a bit to take into account slop in
##   #   the pointing
    $fov *= 1.25;
    if ($debug > 0) {
      printf "fov is %f \n", $fov;
    }

##   # do exactly the same thing for both catalogs ....
    foreach $reftype ("astrom", "photom") {
      if ($debug > 0) {
        printf "top of loop for $reftype\n";
      }
      
##     # open the reference catalog
      $ref_file = $refcat_dir{$reftype} . "/" .$refcat_name{$reftype};
      open(REF_FILE, "$ref_file") || 
               die("refcat.pl: can't open refcat $ref_file");

##     # open the list of images
      $image_file = $output_dir . "/" . $imagelist_name;
      open(IMAGE_FILE, "$image_file") || 
               die("refcat.pl: can't open imagelist $image_file");
       

##     # open the output file, into which we're going to put the stars
##     #   which appear within the image boundaries
      $output_file = $output_dir . "/" . $output_name{$reftype};
      open(OUTPUT_FILE, ">$output_file") || 
               die("refcat.pl: can't open output file $output_file");

##     # fill the "min_ra1", etc. arrays with the limits of the observed fields
       if (create_limit_lists($imagelist_ra_col, $imagelist_dec_col, $fov) != 0) {
         printf STDERR "refcat.pl: create_limit_lists returns with error\n";
         exit(1);
       }
       if ($debug > 0) {
         printf "min_ra1[0] is %10.4f  \n", $min_ra1[0];
       }

##     # now walk through the reference catalog and select only lines containing
##     #   stars which fall within a field
       if (make_subset($ref_ra_col{$reftype}, $ref_dec_col{$reftype}) != 0) {
         printf stderr "refcat.pl: make_subset returns with error\n";
         exit(1);
       }
                   
##     # close all the files
       close(OUTPUT_FILE);
       close(IMAGE_FILE);
       close(REF_FILE);

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

# if all went well, exit with return value 0
exit 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
#
sub create_limit_lists {
    my($fileid, $ra_col, $dec_col, $fov);
    my($DEG2RAD);
    my($line); my(@words);
    my($img); 
    my($ra_cent); my($dec_cent);
    my($cosdec); my($ra_half); my($dec_half);
    my($minr, $maxr);
    my($inr1, $axr1, $inr2, $axr2);
    my($ind, $axd);

    $ra_col = $_[0];
    $dec_col = $_[1];
    $fov = $_[2];

    $min_ra1[0] = -999;
    $max_ra1[0] = -999;
    $min_ra2[0] = -999;
    $max_ra2[0] = -999;
    $min_dec[0] = -999;
    $max_dec[0] = -999;

    $DEG2RAD = 3.14159/180.0;

    if ($debug > 0) {
      printf "entering create_limits_lists"
    }

##   # walk through the list of images, getting the RA, Dec limits
     #    on each one
    $nimages = 0;
    while (<IMAGE_FILE>) {
      $line = " " . $_;
      chomp($line);
      @words = split(/\s+/, $line);

      # skip comments
      if (($words[0] =~ "^#") || ($words[1] =~ "^#")) {
        if ($debug > 0) {
          printf "  comment %s \n", $line;
        }
        next;
      }
   
      # skip lines which aren't images of type "object"
      if ($words[5] ne "object") {
        if ($debug > 0) {
          printf "  not object %s \n", $line;
        }
        next;
      }

      $img = $words[1];
      $ra_cent = $words[6];
      $dec_cent = $words[7];
      
    
     if ($debug > 0) {
       printf "  next image is $img  ra $ra_cent  dec $dec_cent \n";
     }

##     # we calculate the number of degrees of RA covered by this field
     $cosdec = cos($dec_cent*$DEG2RAD);
     $ra_half = ($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)
     $minr = $ra_cent - $ra_half;
     $maxr = $ra_cent + $ra_half;
     if ($minr < 0) {
       $inr1 = 0;
       $axr1 = 0;
       $inr2 = 360.0 + $minr;
       $axr2 = 360;
     } elsif ($maxr > 360) {
       $inr1 = 0;
       $axr1 = $maxr - 360.0;
       $inr2 = $minr;
       $axr2 = 360.0;
     }
     else {
       $inr1 = $minr;
       $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
       $inr2 = 400;
       $axr2 = 400;
     }

##     # now, we calculate the limits in Dec.
##     #   for convenience, place the Dec
##     #            ind, axd                min, max Dec
     $dec_half = $fov/2.0;
     $ind = $dec_cent - $dec_half;
     $axd = $dec_cent + $dec_half;

##     # now, we add these values to the growing arrays of limits
     $min_ra1[$nimages] = $inr1;
     $max_ra1[$nimages] = $axr1;
     $min_ra2[$nimages] = $inr2;
     $max_ra2[$nimages] = $axr2;
     $min_dec[$nimages] = $ind;
     $max_dec[$nimages] = $axd;
     if ($debug > 0) {
       printf "  limits are %6.1f %6.1f  (%6.1f %6.1f)   %6.1f %6.1f \n",
                $min_ra1[$nimages], $max_ra1[$nimages], 
                $min_ra2[$nimages], $max_ra2[$nimages], 
                $min_dec[$nimages], $max_dec[$nimages];
     }
     $nimages++;


    }
  
  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
#
sub make_subset {
    my($ra_col, $dec_col);
    my($nfields);
    my($i);
    my($line); my(@words);
    my($count); my($subcount);
    my($ra, $dec);
    my($field);
  

    $ra_col = $_[0];
    $dec_col = $_[1];

    if ($debug > 0) {
      printf "entering make_subset\n";
    }

    $nfields = $nimages;

    if ($debug > 0) {
      for ($i = 0; $i < $nimages; $i++) {
        printf "make_subset: RA %6.1f %6.1f (%6.1f %6.1f) Dec %6.1f %6.1f \n",
                 $min_ra1[$i], $max_ra1[$i], $min_ra2[$i], $max_ra2[$i],
                 $min_dec[$i], $max_dec[$i];
      }
    }


##   # we walk through the reference catalog, checking one star at a time
    $count = 0;
    $subcount = 0;
    while (<REF_FILE>) {
      $line = " " . $_;
      # don't chomp line, we'll need the newline later if we print it out
      @words = split(/\s+/, $line);
     
      $count++;
      if ($debug > 0) {
        if ($count % 1000 == 0) {
          printf "count $count \n";
        }
      }

##     # skip lines starting with the comment character '#'
      if (substr($words[1], 0, 1) eq "#") {
        next;
      }

##     # sanity check
      if (($#words < $ra_col) || ($#words < $dec_col)) {
        printf stderr "make_subset: skipping bad line ..$line..\n";
        next;
      }

      $ra = $words[$ra_col + 1];
      $dec = $words[$dec_col + 1];
      if ($debug > 2) {
         printf "  ra $ra  dec $dec \n";
      }
## 
       for ($field = 0; $field < $nfields; $field++) {
         if ( ($dec >= $min_dec[$field]) && ($dec <= $max_dec[$field]) ) {
           if ( (($ra >= $min_ra1[$field]) && ($ra <= $max_ra1[$field]))  ||
                (($ra >= $min_ra2[$field]) && ($ra <= $max_ra2[$field])) ) { 
 
             print OUTPUT_FILE $line;
             $subcount++;
             last;
           }
         }
       }

     }
## 
     if ($debug > 0) {
       printf "make_subset: select $subcount of $count stars \n";
     }

  return(0);

}