package      Astangles;
require      Exporter;
@ISA       = qw(Exporter);
@EXPORT    = qw(median radec_dist dec2hms dms2s);
#@EXPORT_OK = qw();

# Subroutines for handling funny astronomical angle conventions, and doing
# calculations with the results.

# Adapted by Rob Reid from:

###############################################################################
#
# align_evt: Align event files prior to merging, based on detected X-ray sources
# Copyright (C) 2001  Tom Aldcroft
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 2
# of the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
#$VERSION = '$Id: align_evt,v 1.1.1.1 2002/05/21 18:47:34 aldcroft Exp $';  # ';
#
# http://cxc.harvard.edu/cal/ASPECT/align_evt/align_evt
#
#############################################################################

# I took out the X ray stuff (we don't have event files anyway) and turned the
# remainder into a library.

use POSIX;

# Set up some constants
my $d2r = 3.14159265358979/180.;
my $d2a = 3600.0;

sub median {
  my @a = sort {$a <=> $b} @_;
  my $mid = int(@a / 2);
  return (@a % 2) ? $a[$mid] : ($a[$mid-1]+$a[$mid])/2.0;
}

# Returns the distance in degrees between (a1, d1) and (a2, d2), which are both
# in degrees.
sub radec_dist {
  my ($a1, $d1, $a2, $d2) = @_;

  return(0.0) if ($a1==$a2 && $d1==$d2);

  return acos( cos($d1*$d2r)*cos($d2*$d2r) * cos(($a1-$a2)*$d2r) +
	       sin($d1*$d2r)*sin($d2*$d2r)) / $d2r;
}

# Converts 1d 2m 3s, or 1D2'3", to 3723.0 (i.e. a bare number) in arcseconds.
# Written by Rob Reid, but inspired by dec2hms below.
sub dms2s {
  $_ = join ' ', @_;
  s/[dms'"]/ /gi;
  @arg = split;
  my $as = 0.0;
  foreach my $d (@arg){
    $as *= 60.0;
    $as += $d;
  }
  return $as;
}

# Converts between sexigesimal (HMS) and decimal RA and Dec.  The
# direction of conversion is given by the number and form of inputs
# Returns two-element array of (RA, Dec) in either case.
sub dec2hms {
  $_ = join ' ', @_;
  s/[,:dhms'"]/ /gi;
  @arg = split;

  if (@arg == 2) {
    my $ra = shift;
    my $dec = shift;
    my ($rah, $ram, $ras);
    my ($dec_sign, $decd, $decm, $decs);
    my ($ra_hms, $dec_hms);

    $ra += 360.0 if ($ra < 0);
    $ra /= 15.;
    $rah = floor($ra);
    $ram = floor(($ra - $rah) * 60.);
    $ras = ($ra - $rah - $ram / 60.) * 60. * 60.;

    $dec_sign = ($dec < 0);
    $dec = abs($dec);
    $decd = floor($dec);
    $decm = floor(($dec - $decd ) * 60.);
    $decs = ($dec - $decd - $decm / 60) * 60. * 60.;

    $ra_hms = sprintf "%02dh%02dm%06.3fs", $rah, $ram, $ras;
    $dec_hms = sprintf "%s%02dd%02dm%05.2fs", $dec_sign ? '-' : ' ',
      $decd, $decm, $decs;

    return ($ra_hms,$dec_hms);
  } elsif (@arg == 6) {
    my ($rah, $ram, $ras, $decd, $decm, $decs) = @arg;
#    print "rah: $rah, ram: $ram, ras: $ras, decd: $decd, decm: $decm, decs: $decs\n";
    $ra = 15.0 * ($rah + $ram / 60.0 + $ras / 3600.0);
    $dec = abs($decd) + $decm / 60.0 + $decs / 3600.0;
    $dec = -$dec if ($decd < 0.0);
    return (sprintf("%12.7f", $ra), sprintf("%12.6f", $dec));
  } else {
    warn "dec2hms: Error -- enter either 6 or 2 arguments\n";
  }
}

# dec2hms toggles between hms and decimal degrees.  This makes sure that its
# output is in hms, d'" whether or not its input is already in that format or
# in decimal degrees.  Expects a single string as input.
sub ensurehms {
  unless($_[0] =~ /[Hh]/) {
    return dec2hms(@_);
  }
  return @_;
}

# dec2hms toggles between hms and decimal degrees.  This makes sure that its
# output is in decimal degrees whether or not its input is already in that
# format or in hms, d'".  Expects a single string as input.
sub ensuredecdeg {
  if($_[0] =~ /[Hh]/){
    return dec2hms(@_);
  }
  return @_;
}

sub xcorr_sources {
  $evt2 = shift;
  $src2 = $src2{$evt2};
  @celldet = get_cell_detect($src2);

  foreach (@celldet) {
    ($ra, $dec, $counts, $snr) = split;
    next unless ($snr >= $par{snr});

    foreach $cat (@cat) {
      $delta_rad = 3600 * radec_dist($ra, $dec, $cat->{ra}, $cat->{dec});
      next unless ($delta_rad < $par{match_rad});
      push @{$match{$src2}}, { ra => $ra,
			       dec => $dec,
			       snr => $snr,
			       d_dec => ($cat->{dec} - $dec) * $d2a,
			       d_ra  => ($cat->{ra} - $ra) * cos($cat->{dec}*$d2r) * $d2a,
			     };
      last;
    }
  }

  @d_ra = ();
  @d_dec= ();
  print "\nFollowing sources in $evt2 match baseline catalog:\n" if ($par{loud});
  printf("%10s %10s %6s %6s %6s\n", "RA", "Dec", "D_RA", "D_Dec", "SNR") if ($par{loud});
  foreach $m (@{$match{$src2}}) {
    printf("%10.5f %10.5f %6.2f %6.2f %6.2f\n", $m->{ra}, $m->{dec},
	   $m->{d_ra}, $m->{d_dec}, $m->{snr}) if ($par{loud});
    push @d_ra, $m->{d_ra};
    push @d_dec, $m->{d_dec};
  }
  print "\n" if ($par{loud});
    
  return (@d_ra) ? (median(@d_ra), median(@d_dec)) : (0.0, 0.0)
}