#
# $Name:  $
#
# $Log: FITS.pm,v $
# Revision 1.4  2002/02/19 14:48:10  robc
# Perltidy
#
# Revision 1.3  2001/11/23 16:54:01  robc
# Add FITS_get_pixel_size and FITS_remove_data
#
# Revision 1.2  2001/10/06 09:59:46  robc
# Add FITS_get_date_time and FITS_find_matched.
# Change FITS_get_ra_dec to use CRVAL1, CRVAL2 over RA, DEC since
#    Astro::Time::str2rad doesn't always seem to work correctly.
#
# Revision 1.1  2001/09/24 14:08:16  robc
# *** empty log message ***
#

use strict;
use CFITSIO qw( :longnames :constants );
use Astro::Time;
use Astro::Coord;
use Math::Trig qw( asec );
use Carp qw( cluck );

package FITS;
use Exporter();

use vars qw(
   @ISA
   @EXPORT
   %file_hash
   %file_valid_fits
   %file_calc
   );

@ISA = qw( Exporter );

@EXPORT = qw(
   FITS_find_matched
   FITS_get_az_elv
   FITS_get_ra_dec
   FITS_get_airmass
   FITS_get_zenith_distance
   FITS_get_date_time
   FITS_get_hour_angle
   FITS_get_all_keywords
   FITS_add_keyword
   FITS_find_fits
   FITS_get_key_value
   FITS_get_pixel_size
   FITS_get_size
   FITS_group_keyword
   FITS_remove_data
   FITS_key_exists
   );

sub FITS_find_matched( $$ )
   {
   my ( $dir, $max_diff ) = @_;

   my @fits_files      = FITS_find_fits( $dir );
   my $file_imgtyp_ref = FITS_group_keyword( 'IMAGETYP', @fits_files );
   my $file_filter_ref = FITS_group_keyword( 'FILTER', @fits_files );
   my %filter_time;
   my %file_time;
   my %test;
   my %reverse_filter = reverse %$file_filter_ref;
   my %matched;
   my $index = 0;
   my @times;
   my @filters;

   die "Error: Don't support other than 2 filters in a directory"
      if keys %reverse_filter != 2;

   foreach my $filter ( sort keys %reverse_filter )
      {
      foreach my $file ( keys %$file_imgtyp_ref )
         {
         next unless $file_imgtyp_ref->{ $file } =~ /^OBJECT/io;
         next unless $file_filter_ref->{ $file } =~ /$filter/;
         my @tmp = split /:/, @{ FITS_get_date_time( $file ) }[ 1 ];
         $file_time{ $filter }{ $file } =
            Astro::Time::hms2time( $tmp[ 0 ], $tmp[ 1 ], $tmp[ 2 ] );
         }
      my @data = sort values %{ $file_time{ $filter } };
      $filters[ $index ] = $filter;
      $times[ $index++ ] = \@data;
      }
   my %reverse_time_1 = reverse %{ $file_time{ $filters[ 0 ] } };
   my %reverse_time_2 = reverse %{ $file_time{ $filters[ 1 ] } };
   while ( scalar @{ $times[ 0 ] } > 0 )
      {
      my $time     = pop @{ $times[ 0 ] };
      my $min      = abs( $time - $times[ 1 ][ 0 ] );
      my $min_time = $times[ 1 ][ 0 ];
      for ( my $x = 1; $x < scalar @{ $times[ 1 ] }; ++$x )
         {
         my $diff = abs( $time - $times[ 1 ][ $x ] );
         if ( $diff < $min )
            {
            $min      = $diff;
            $min_time = $times[ 1 ][ $x ];
            }
         last if $min == 0;
         }
      if ( $min > $max_diff )
         {
         printf "No match for %s minimum time diff: %d seconds\n",
            $reverse_time_1{ $time }, $min * 24 * 3600;
         printf "\tfrom %s\n", $reverse_time_2{ $min_time };
         }
      else
         {
         $matched{ $reverse_time_1{ $time } } = $reverse_time_2{ $min_time };
         }
      }
   return \%matched;
   }

sub FITS_get_az_elv( $ )
   {
   my ( $file ) = @_;

   FITS_grab_header( $file );

   if ( !exists $file_calc{ $file }{ az } && $file_valid_fits{ $file } )
      {
      my ( $tra, $tdec ) = FITS_get_ra_dec( $file );
      my $dec = Astro::Time::rad2turn( $tdec );
      my $ha  = Astro::Time::rad2turn( FITS_get_hour_angle( $file ) );
      my $lat =
         Astro::Time::deg2turn( eval FITS_get_key_value( $file, 'LATITUDE' ) );
      my ( $taz, $telv ) = Astro::Coord::eqazel( $ha, $dec, $lat );
      $file_calc{ $file }{ az }  = Astro::Time::turn2rad( $taz );
      $file_calc{ $file }{ elv } = Astro::Time::turn2rad( $telv );
      }
   elsif ( !$file_valid_fits{ $file } )
      {
      die "Error: $file not a fits file";
      }

   return $file_calc{ $file }{ az }, $file_calc{ $file }{ elv };
   }

sub FITS_get_ra_dec( $ )
   {
   my ( $file ) = @_;

   FITS_grab_header( $file );

   if ( !exists $file_calc{ $file }{ ra } && $file_valid_fits{ $file } )
      {
      if ( !exists $file_hash{ $file }->{ CRVAL2 } )
         {
         $file_calc{ $file }{ ra } =
            Astro::Time::str2rad( eval FITS_get_key_value( $file, 'RA' ), 'H' );
         $file_calc{ $file }{ dec } =
            Astro::Time::str2rad( eval FITS_get_key_value( $file, 'DEC' ),
                                  'H' );
         }
      else
         {
         $file_calc{ $file }{ ra } =
            Astro::Time::deg2rad( eval FITS_get_key_value( $file, 'CRVAL2' ) );
         $file_calc{ $file }{ dec } =
            Astro::Time::deg2rad( eval FITS_get_key_value( $file, 'CRVAL1' ) );
         }
      }
   elsif ( !$file_valid_fits{ $file } )
      {
      die "Error: $file not a fits file";
      }

   return $file_calc{ $file }{ ra }, $file_calc{ $file }{ dec };
   }

sub FITS_get_airmass( $ )
   {
   my ( $file ) = @_;

   FITS_grab_header( $file );

   if ( !exists $file_calc{ $file }{ airmass } && $file_valid_fits{ $file } )
      {
      my $zenith = FITS_get_zenith_distance( $file );

      #
      # Calculate the Airmass from the zenith distance per
      # http://star-www.st-and.ac.uk/starlink/stardocs/sc6.htx/node15.html
      $file_calc{ $file }{ airmass } =
         $zenith - 0.0018167 * ( $zenith - 1 ) - 0.0028750 * ( $zenith - 1 )
         **2 - 0.0008083 * ( $zenith - 1 )**3;

      #
      # Another formula for calculating the airmass from the above page
      # Trig identies from
      # http://forum.swarthmore.edu/dr.math/faq/formulas/faq.trig.html
      my $za   = Math::Trig::asec( $zenith );
      my $sec2 = 1 + ( 1 - cos( 2 * $za ) ) / ( 1 + cos( 2 * $za ) );
      $file_calc{ $file }{ airmass2 } =
         $zenith * ( 1 - 0.0012 * ( $sec2 - 1 ) );
      }
   elsif ( !$file_valid_fits{ $file } )
      {
      die "Error: $file not a fits file";
      }

   return $file_calc{ $file }{ airmass }, $file_calc{ $file }{ airmass2 };
   }

sub FITS_get_zenith_distance( $ )
   {
   my ( $file ) = @_;
   my $zenith;

   FITS_grab_header( $file );

   if ( !exists $file_calc{ $file }{ zenith } && $file_valid_fits{ $file } )
      {
      my $hour_angle = FITS_get_hour_angle( $file );
      my $lat        =
         Astro::Time::deg2rad( eval FITS_get_key_value( $file, 'LATITUDE' ) );
      my ( $ra, $dec ) = FITS_get_ra_dec( $file );

      #
      # Zenith distance per
      # http://star-www.st-and.ac.uk/starlink/stardocs/sc6.htx/node32.html
      $zenith =
         ( 1 / ( sin( $lat ) * sin( $dec ) + cos( $lat ) * cos( $dec ) *
         cos( $hour_angle ) ) );

      #
      # Set up for repeat offenders
      $file_calc{ $file }{ zenith } = $zenith;
      }
   elsif ( !$file_valid_fits{ $file } )
      {
      die "Error: $file not a fits file";
      }

   return $file_calc{ $file }{ zenith };
   }

sub FITS_get_date_time( $ )
   {
   my ( $file ) = @_;

   FITS_grab_header( $file );

   if ( !exists $file_calc{ $file }{ date_time } && $file_valid_fits{ $file } )
      {
      my @data = split /T/, eval FITS_get_key_value( $file, 'DATE-OBS' );

      #
      # If only one entry, it's the date, so we need to get the time also
      if ( scalar @data == 1 )
         {
         push @data, eval FITS_get_key_value( $file, 'UT' );
         }
      $file_calc{ $file }{ date_time } = \@data;
      }
   elsif ( !$file_valid_fits{ $file } )
      {
      die "Error: $file not a fits file";
      }

   return $file_calc{ $file }{ date_time };
   }

sub FITS_get_hour_angle( $ )
   {
   my ( $file ) = @_;
   my $hour_angle = 0;

   FITS_grab_header( $file );

   if ( !exists $file_calc{ $file }{ hour_angle } && $file_valid_fits{ $file } )
      {
      my @data = split /T/, eval FITS_get_key_value( $file, 'DATE-OBS' );

      #
      # If only one entry, it's the date, so we need to get the time also
      if ( scalar @data == 1 )
         {
         push @data, eval FITS_get_key_value( $file, 'UT' );
         }
      my $long =
         Astro::Time::deg2turn( eval FITS_get_key_value( $file, 'LONGITUD' ) );
      my ( $ra, $dec ) = FITS_get_ra_dec( $file );
      my @date = split /-/, $data[ 0 ];
      my @time = split /:/, $data[ 1 ];
      my $ut = Astro::Time::hms2time( $time[ 0 ], $time[ 1 ], $time[ 2 ] );
      my $lst = Astro::Time::cal2lst(
                                      $date[ 2 ], $date[ 1 ], $date[ 0 ],
                                      $ut,        $long
                                      );

      #
      # Hour angle is lst - ra per
      # http://star-www.st-and.ac.uk/starlink/stardocs/sc6.htx/node32.html
      $hour_angle = Astro::Time::turn2rad( $lst ) - $ra;
      while ( $hour_angle < 0 )
         {
         $hour_angle += Math::Trig::pi;
         }

      #
      # Set up for repeat offenders
      $file_calc{ $file }{ hour_angle } = $hour_angle;
      }
   elsif ( !$file_valid_fits{ $file } )
      {
      die "Error: $file not a fits file";
      }

   return $file_calc{ $file }{ hour_angle };
   }

sub FITS_get_all_keywords( $ )
   {
   my ( $file ) = @_;

   FITS_grab_header( $file );

   return $file_valid_fits{ $file }, $file_hash{ $file };
   }

sub FITS_add_keyword( $$$ )
   {
   my ( $file, $keyword, $value ) = @_;

   FITS_grab_header( $file );
   if ( exists $file_hash{ $file }->{ $keyword } )
      {
      warn "Keyword \'$keyword\' already exists - use FITS_update_keyword";
      }
   else
      {
      my $status = 0;
      my $orig;
      my $fptr = CFITSIO::open_file( $file, CFITSIO::READWRITE, $status );
      die "Failed to open $file for READWRITE: $status\n" unless !$status;

      if ( $value =~ /^[-.\d]+$/o )
         {
         $fptr->write_key_flt( $keyword, $value,
                         -( length( $value ) - length( int( $value ) ) ), undef,
                         $status );
         }
      else
         {
         $fptr->write_key_str( $keyword, $value, undef, $status );
         }
      die "Failed to write key $keyword to $value: $status\n" unless !$status;
      $fptr->close_file( $status );
      die "Managed to screw up closing file $file: $status\n" unless !$status;
      }
   }

sub FITS_group_keyword( $@ )
   {
   my ( $keyword, @files ) = @_;
   my %groups;

   for ( @files )
      {
      FITS_grab_header( $_ );
      if ( $file_valid_fits{ $_ } )
         {
         if ( FITS_key_exists( $_, $keyword ) )
            {
            my $value = uc eval FITS_get_key_value( $_, $keyword );
            chomp $value;
            $value =~ s/\s+$//g;
            $groups{ $_ } = $value;
            }
         }
      }
   return \%groups;
   }

sub FITS_find_fits( $ )
   {
   my ( $dir ) = @_;
   my @valid;

   die "$dir doesn't exist" unless opendir( DIR, $dir );
   my @files = map $dir . "/" . $_, sort readdir( DIR );
   closedir( DIR );

   for ( @files )
      {
      FITS_grab_header( $_ );
      push @valid, $_ if $file_valid_fits{ $_ };
      }

   return @valid;
   }

sub FITS_get_key_value( $$ )
   {
   my ( $file, $key ) = @_;

   FITS_grab_header( $file );
   die "Error: $file doesn't contain key $key"
      unless $file_hash{ $file }->{ uc $key };

   return $file_hash{ $file }->{ uc $key };
   }

sub FITS_get_pixel_size( $ )
   {
   my $file = shift;

   FITS_grab_header( $file );
   if ( !$file_hash{ $file }->{ CDELT1 } and !$file_hash{ $file }->{ CDELT2 } )
      {
      warn "\nWarning: using default values for pixel size\n";
      $file_hash{ $file }->{ CDELT1 } = 0.002126;
      $file_hash{ $file }->{ CDELT2 } = -0.002126;
      }
   return ( $file_hash{ $file }->{ CDELT1 }, $file_hash{ $file }->{ CDELT2 } );
   }

sub FITS_get_size( $ )
   {
   my $file = shift;

   FITS_grab_header( $file );
   return ( $file_hash{ $file }->{ NAXIS1 }, $file_hash{ $file }->{ NAXIS2 } );
   }

sub FITS_remove_data( $ )
   {
   my $file = shift;

   delete $file_hash{ $file };
   delete $file_valid_fits{ $file };

   return;
   }

sub FITS_key_exists( $$ )
   {
   my ( $file, $key ) = @_;

   FITS_grab_header( $file );
   return exists $file_hash{ $file }->{ uc $key };
   }

#
# Functions not exported
#
sub FITS_grab_header( $ )
   {
   my $file   = shift;
   my $status = 0;
   my %temp;

   if ( not exists $file_hash{ $file } )
      {
      open( FILE, $file ) or die "Cannot open $file";
      my $start = <FILE>;
      if ( defined $start and $start =~ /(SIMPLE|XTENSION)/ )
         {
         ( $file_hash{ $file }, $status ) = CFITSIO::fits_read_header( $file );
         }
      else
         {
         $status = 1;
         }
      if ( $status )
         {
         $file_valid_fits{ $file } = 0;
         }
      else
         {
         $file_valid_fits{ $file } = 1;
         }
      }
   }