#!/usr/bin/perl
#
# Finds the best match for the given star/raw data files
#
# $Name:  $
#
# $Log: matchStars.pl,v $
# Revision 1.5  2002/01/29 04:55:33  robc
# General re-write so matching is 'deterministic'.
# Add support for AppConfig config file parsing.
# Remove mounds of old code.
# Perltidy.
#
# Revision 1.4  2001/12/09 01:36:08  robc
# Change to use pixel scale from fits header rather than hard coded values
#
# Revision 1.3  2001/11/14 04:01:15  robc
# Perltidy code
#
# Revision 1.2  2001/10/06 10:15:34  robc
# General hacking to try to get it working correctly.
#
# Revision 1.1  2001/09/23 17:31:41  robc
# *** empty log message ***
#

use warnings;
use strict;
use lib( $ENV{ TASSIV_FITS_PM } );

use CFITSIO qw( :longnames :constants );
use IO::Socket;
use Math::NumberCruncher;
use FITS;
use Astro::Time;
use DB_File;
use AppConfig qw( :argcount :expand );

autoflush STDOUT;

my $ID = q$Id: matchStars.pl,v 1.5 2002/01/29 04:55:33 robc Exp $;
my $version = join ( ' ', ( split ( ' ', $ID ) )[ 1 .. 3 ] );
$version =~ s/,v\b//;
$version =~ s/(\S+)$/($1)/;

sub find_best( $$$$$$$$ );
sub initialize_star_catalog( $$$ );
sub find_subset( $$$$$$ );
sub calculate_perc( $$$$ );
sub setup_args( $$ );

die "Please specify the fit file to use\n" unless my $fit_file = shift;

my $CONFIG_FILE = '/tass/src/tassiv_reduce/tassiv_reduce.cfg';

my $c = AppConfig->new(
   {
      CREATE => 1,
      GLOBAL => {
                  EXPAND   => EXPAND_VAR,
                  ARGCOUNT => ARGCOUNT_ONE
      }
   }
   );
$c = setup_args( $c, $CONFIG_FILE );

my $star_file = $fit_file . ".raw";
my $dest_file = $fit_file . ".mat";

die "$star_file doesn\'t exist" unless -e $star_file;
die "$dest_file already exists" if -e $dest_file;
print "$fit_file...\n";

my ( $file_ra, $file_dec ) = FITS_get_ra_dec( $fit_file );
$file_ra  = rad2deg( $file_ra );
$file_dec = rad2deg( $file_dec );
my ( $size_x, $size_y )  = FITS_get_size( $fit_file );
my ( $ra_cov, $dec_cov ) = FITS_get_pixel_size( $fit_file );
$ra_cov  = abs( $ra_cov * $size_x );
$dec_cov = abs( $dec_cov * $size_y );

my $start_ra  = $file_ra - $ra_cov * $c->search_ra;
my $ra_inc    = 2 * $ra_cov * $c->search_ra / $c->div_ra;
my $start_dec = $file_dec - $dec_cov * $c->search_dec;
my $dec_inc   = 2 * $dec_cov * $c->search_dec / $c->div_dec;

my $match_star = 'match_star.list';
my $match_cat  = 'match_cat.list';
my $temp_file  = $fit_file . '.temp';

my $cmd = 'sort -n -k ' . ( $c->mag_col + 1 ) . " $star_file | head -";
$cmd .= $c->keep_process . " > $match_star";
`$cmd`;

my ( $band ) = uc eval FITS_get_key_value( $fit_file, 'FILTER' );
my ( $real_ra, $real_dec ) = ( $file_ra, $file_dec );
my ( $x, $y, $per );
my $done = 0;

my ( %star_ra, %star_dec, @star_catalog );
$star_ra{ min }  = $start_ra - $ra_cov;
$star_ra{ max }  = $file_ra + $ra_cov * $c->search_ra + $ra_cov;
$star_dec{ min } = $start_dec - $dec_cov;
$star_dec{ max } = $file_dec + $dec_cov * $c->search_dec + $dec_cov;

@star_catalog = initialize_star_catalog( $band, \%star_ra, \%star_dec );

$c->div_ra, $c->div_dec;

for ( my $x = 0; $x < $c->iterations; ++$x )
   {
   ( $real_ra, $real_dec, $per ) =
      find_best( $start_ra, $start_dec, $ra_cov, $dec_cov, $ra_inc, $dec_inc,
                 $c->div_ra, $c->div_dec );
   last if $per < $c->keep_percent;
   $ra_inc  /= $c->inc_reduction;
   $dec_inc /= $c->inc_reduction;
   $start_ra  = $real_ra - $c->div_ra * $ra_inc / 2;
   $start_dec = $real_dec - $c->div_dec * $dec_inc / 2;
   }

if ( $per > $c->keep_percent )
   {
   printf( "Finding final match using RA: %6.2f (%6.2f) DEC: %6.2f (%6.2f)\n",
           $real_ra, $file_ra, $real_dec, $file_dec );

   $cmd = 'sort -n -k ' . ( $c->mag_col + 1 ) . " $star_file | head -";
   $cmd .= $c->keep_final . " > $match_star";
   `$cmd`;

   find_subset( $real_ra, $real_dec, $ra_cov, $dec_cov, $match_cat,
                $c->keep_final );

   my $cmd = "match $match_star";
   $cmd .= ' ' . $c->x_col;
   $cmd .= ' ' . $c->y_col;
   $cmd .= ' ' . $c->mag_col;
   $cmd .= " $match_cat 1 2 3";
   $cmd .= ' linear nobjs=' . $c->keep_final;
   $cmd .= ' triad=0.0015 matchrad=0.0000242';

   #      $cmd .= " recalc max_iter=3";
   #      $cmd .= " halt_sigma=1.0e-11 scale=1.0";

   my $best_output = `$cmd`;

   `wc -l matched.mtA` =~ /(\d+)/;
   my $match_count = $1;

   $best_output =~ /a\=([-.\d]+)\s+
                    b\=([-.\d]+)\s+
                    c\=([-.\d]+)\s+
                    d\=([-.\d]+)\s+
                    e\=([-.\d]+)\s+
                    f\=([-.\d]+)\s+
                   /x;
   my $a = $1;
   my $b = $2;
   my $c = $3;
   my $d = $4;
   my $e = $5;
   my $f = $6;

   my $sum1 = abs( $b ) + abs( $f );
   my $sub1 = abs( abs( $b ) - abs( $f ) );
   my $sum2 = abs( $c ) + abs( $e );
   my $sub2 = abs( abs( $c ) - abs( $e ) );

   if ( $sum1 == 0 or $sum2 == 0 )
      {
      print "Very bad match\n";
      }
   else
      {
      my $per1 = 100 * $sub1 / $sum1;
      my $per2 = 100 * $sub2 / $sum2;
      if ( ( $per1 > 10 ) || ( $per2 > .2 ) )
         {
         printf "Bad match %5.1f%% & %5.1f%%\n", $per1, $per2;
         }
      else
         {
         printf "Good match %5.1f%% & %5.1f%%\n", $per1, $per2;

         tie my %thash, "DB_File", $dest_file
            or die "Error: cannot tie hash to $dest_file";
         $thash{ type }   = "Linear";
         $thash{ a }      = $a;
         $thash{ b }      = $b;
         $thash{ c }      = $c;
         $thash{ d }      = $d;
         $thash{ e }      = $e;
         $thash{ f }      = $f;
         $thash{ ra }     = $real_ra;
         $thash{ dec }    = $real_dec;
         $thash{ count }  = $match_count;
         $thash{ to_cat } = 1;
         untie %thash;
         }
      }
   }
else
   {
   open( FILE, ">> matchStars.bad" );
   printf FILE "$star_file matched %6.2f%%\n", $per;
   close( FILE );
   printf( "Not enough stars matched %5.2f%%\n\n", $per );
   }

sub find_best( $$$$$$$$ )
   {
   my (
        $start_ra, $start_dec, $ra_cov, $dec_cov,
        $ra_inc,   $dec_inc,   $end_x,  $end_y
      )
      = @_;

   my ( $x, $y, $percent, $output, $use_ra, $use_dec );
   my ( $best_percent, $bra, $bdec ) = ( 0, $start_ra, $start_dec );

   printf(
           "%6.2f->%6.2fx%4.2f %5.2f->%5.2fx%4.2f ",
           $start_ra,
           $start_ra + $ra_inc * $end_x,
           $ra_inc,
           $start_dec,
           $start_dec + $dec_inc * $end_y,
           $dec_inc
           );

   for ( $x = 0; $x < $end_x; ++$x )
      {
      $use_ra = $start_ra + $x * $ra_inc;
      for ( $y = 0; $y < $end_y; ++$y )
         {
         $use_dec = $start_dec + $y * $dec_inc;

         ( $percent, $output ) =
            calculate_perc( $use_ra, $use_dec, $ra_cov, $dec_cov );

         if ( ( $percent > $best_percent ) )
            {
            $best_percent = $percent;
            $bra          = $use_ra;
            $bdec         = $use_dec;
            }
         }
      printf( "." );
      }
   printf( " %6.2f%% %6.2f %6.2f\n", $best_percent, $bra, $bdec );

   return ( $bra, $bdec, $best_percent );
   }

sub initialize_star_catalog( $$$ )
   {
   my ( $band, $rah, $dech ) = @_;
   my ( @catalog, @data, $cmd );

   my $ra   = $rah->{ min } + ( $rah->{ max } - $rah->{ min } ) / 2;
   my $dec  = $dech->{ min } + ( $dech->{ max } - $dech->{ min } ) / 2;
   my $dra  = ( $rah->{ max } - $ra );
   my $ddec = ( $dech->{ max } - $dec );

   print "Initializing star catalog for:\n";
   printf "\tRA: %6.2f +-%5.2f DEC: %6.2f +-%5.2f...\n", $ra, $dra, $dec, $ddec;

   $cmd = "findStars.pl $band $ra $dec $dra $ddec |";

   open FILE, "> test.out";

   open( DATA, $cmd );
   while ( <DATA> )
      {
      @data = split ( /\s+/, $_ );
      push @catalog, [ @data ];
      print FILE join ' ', @data, "\n";
      }
   close FILE;

   @catalog = sort( { $$a[ 3 ] <=> $$b[ 3 ] } @catalog );

   return ( @catalog );
   }

sub find_subset( $$$$$$ )
   {
   my ( $ra, $dec, $ra_cov, $dec_cov, $file, $keep ) = @_;
   my $star_count = 0;

   my @ra_min  = $ra - $ra_cov / 2;
   my @ra_max  = $ra + $ra_cov / 2;
   my $dec_min = $dec - $dec_cov / 2;
   my $dec_max = $dec + $dec_cov / 2;

   if ( ( $ra_min[ 0 ] < $star_ra{ min } ) || ( $ra_max[ 0 ] > $star_ra{ max } )
        || ( $dec_min < $star_dec{ min } ) || ( $dec_max > $star_dec{ max } ) )
      {
      print( "\n$ra_min[ 0 ] < $star_ra{ min }\n" );
      print( "$ra_max[ 0 ] > $star_ra{ max }\n" );
      print( "$dec_min < $star_dec{ min }\n" );
      print( "$dec_max > $star_dec{ max }\n" );
      die "\nError: didn't initialize catalog with enough stars\n";
      }

   my $split_range = 0;

   if ( $ra_min[ 0 ] < 0 )
      {
      $ra_min[ 1 ] = $ra_min[ 0 ] + 360;
      $ra_min[ 0 ] = 0;
      $ra_max[ 1 ] = 360;
      $split_range = 1;
      }
   elsif ( $ra_max[ 0 ] > 360 )
      {
      $ra_min[ 1 ] = 0;
      $ra_max[ 1 ] = $ra_max[ 0 ] - 360;
      $ra_max[ 0 ] = 360;
      $split_range = 1;
      }

   open( STAR, "> $temp_file" );

   if ( $split_range )
      {
      foreach my $star ( @star_catalog )
         {
         if ( ( $$star[ 2 ] >= $dec_min ) && ( $$star[ 2 ] <= $dec_max ) )
            {
            if ( ( ( $$star[ 1 ] >= $ra_min[ 0 ] )
               && ( $$star[ 1 ] <= $ra_max[ 0 ] ) )
               || ( ( $$star[ 1 ] >= $ra_min[ 1 ] )
                    && ( $$star[ 1 ] <= $ra_max[ 1 ] ) ) )
               {
               print STAR join ( "\t", @$star ) . "\n";
               last if ( $keep && ( ++$star_count >= $keep ) );
               }
            }
         }
      }
   else
      {
      foreach my $star ( @star_catalog )
         {
         if ( ( $$star[ 2 ] >= $dec_min ) && ( $$star[ 2 ] <= $dec_max ) )
            {
            if ( ( $$star[ 1 ] >= $ra_min[ 0 ] )
                 && ( $$star[ 1 ] <= $ra_max[ 0 ] ) )
               {
               ++$star_count;
               print STAR join ( "\t", @$star ) . "\n";
               last if ( $keep && ( ++$star_count >= $keep ) );
               }
            }
         }
      }
   close( STAR );

   `project_coords $temp_file 1 2 $ra $dec outfile=$file`;
   unlink( $temp_file );

   return ( $star_count );
   }

sub calculate_perc( $$$$ )
   {
   my ( $ra, $dec, $ra_cov, $dec_cov ) = @_;
   my ( $match_count, $output ) = ( 0 );

   if (
       find_subset( $ra, $dec, $ra_cov, $dec_cov, $match_cat, $c->keep_process )
      )
      {
      my $cmd = "match $match_star";
      $cmd .= ' ' . $c->x_col;
      $cmd .= ' ' . $c->y_col;
      $cmd .= ' ' . $c->mag_col;
      $cmd .= " $match_cat 1 2 3";
      $cmd .= ' linear nobjs=' . $c->keep_process;
      $cmd .= ' triad=0.0015 matchrad=0.0000242';
      `$cmd`;
      `wc -l matched.mtA` =~ /(\d+)/;
      $match_count = $1;

      #print "--$match_count--\n";
      #print "$cmd\n";
      #exit;
      }
   return ( 100 * $match_count / $c->keep_process );
   }

sub setup_args( $$ )
   {
   my ( $c, $file ) = @_;
   $c->define(
               Match_Stars_search_ra => {
                                          DEFAULT => 1,
                                          ALIAS   => 'search_ra'
               },
               Match_Stars_search_dec => {
                                           DEFAULT => 1,
                                           ALIAS   => 'search_dec'
               },
               Match_Stars_div_ra => {
                                       DEFAULT => 5,
                                       ALIAS   => 'div_ra'
               },
               Match_Stars_div_dec => {
                                        DEFAULT => 5,
                                        ALIAS   => 'div_dec'
               },
               Match_Stars_inc_reduction => {
                                              DEFAULT => 3.5,
                                              ALIAS   => 'inc_reduction'
               },
               Match_Stars_iterations => {
                                           DEFAULT => 3,
                                           ALIAS   => 'iterations'
               },
               Match_Stars_keep_percent => {
                                             DEFAULT => 19,
                                             ALIAS   => 'keep_percent'
               },
               Match_Stars_keep_process => {
                                             DEFAULT => 61,
                                             ALIAS   => 'keep_process'
               },
               Match_Stars_keep_final => {
                                           DEFAULT => 91,
                                           ALIAS   => 'keep_final'
               },
               Match_Stars_mag_column => {
                                           DEFAULT => 2,
                                           ALIAS   => 'mag_col'
               },
               Match_Stars_x_column => {
                                         DEFAULT => 0,
                                         ALIAS   => 'x_col'
               },
               Match_Stars_y_column => {
                                         DEFAULT => 1,
                                         ALIAS   => 'y_col'
               },
               );
   $c->file( $file );

   $c->define( 'help|h!' );
   $c->define( 'version|v!' );
   $c->define( 'force|f!' );

   $c->args( \@ARGV );

   return $c;
   }