#!/usr/bin/perl
#
# reformat the SN catalog so that 
#   RA and Dec are in decimal degrees
#
# MWR 10/17/2007
use POSIX;

$debug = 0;
$DEGTORAD = (3.14159/180.0);

if (0 == 1) {
  $ra = 15.0;
  $dec = 10.0;

  ($glong, $glat) = radec_to_galactic($ra, $dec);
  printf " glong glat  %10.5f %10.5f \n", $glong, $glat;
}






while (<STDIN>) {
  $line = $_;

  if ($debug > 0) {
    printf "next line is ..%s.. \n", $line;
  }
  if ($line =~ /^#/) {
    next;
  }

  # get rid of some characters
  $line =~ s/\(/ /g;
  $line =~ s/\)/ /g;
  $line =~ s/z=/  /g;
 
  # get the name, replace spaces by underscores
  $name_str = substr($line, 0, 7);
  $name_str =~ s/ //g;
  if ($debug > 0) {
    printf " name_str is ..%s.. \n", $name_str;
  }

  # figure out RA, convert to decimal degrees
  $ra_str = substr($line, 19, 10);
  @words = split(/\s+/, $ra_str);
  $ra_h = $words[0];
  $ra_m = $words[1];
  $ra_s = $words[2];
  if ($debug > 0) {
    printf " RA  ..%s.. %d %d %5.2f ", $ra_str, $ra_h, $ra_m, $ra_s;
  }
  $ra_deg = $ra_h + ($ra_m/60.0) + ($ra_s/3600.0);
  $ra_deg *= 15.0;
  if ($debug > 0) {
    printf "  = %10.5f \n", $ra_deg;
  }
  
  # figure out Dec, convert to decimal degrees
  $dec_str = substr($line, 29, 11);
  @words = split(/\s+/, $dec_str);
  $dec_d = $words[0];
  $dec_m = $words[1];
  $dec_s = $words[2];
  if ($debug > 0) {
    printf " Dec  ..%s.. %d %d %5.2f ", $dec_str, $dec_d, $dec_m, $dec_s;
  }
  if ($dec_d >= 0) {
    $dec_sign = 1;
  }
  else {
    $dec_d = 0.0 - $dec_d;
    $dec_sign = -1.0;
  }
  $dec_deg = $dec_d + ($dec_m/60.0) + ($dec_s/3600.0);
  $dec_deg *= $dec_sign;
  if ($debug > 0) {
    printf "  = %10.5f \n", $dec_deg;
  }
  
  # calculate the galactic longitude and latitude for this object
  ($glong, $glat) = radec_to_galactic($ra_deg, $dec_deg);

  # now print out a copy of the original line, but with the
  #   new name, and with RA and Dec in decimal degrees,
  #   and with glong and glat inserted after RA, Dec
  $first = substr($line, 7, 12);
  $last = substr($line, 38, length($line) - 38);
  printf "%-8s %s %10.5f %10.5f %10.5f %10.5f %s", 
           $name_str, $first, $ra_deg, $dec_deg, $glong, $glat, $last;



}




exit 0;




#######################################################################
# Convert a position in (RA, Dec) -- decimal degrees -- 
#    into galactic coordinates (glong, glat) -- also in decimal degrees
#
# 
# usage:   (glong, glat)  =  radec_to_galactic ( ra, dec) 
#
sub radec_to_galactic {

  # coords of galactic North Pole in the RA, Dec system
  my($rag) = 192.85;
  my($decg) = 27.117;
  # don't know what this is
  my($pgl) = 123.0;
  my($pgq, $gpq, $gq, $glong, $glat);

  $ra = $_[0];
  $dec = $_[1];
 
  if ($debug > 0) {
    printf " ra %10.5f  dec %10.5f \n", $ra, $dec;
  }

  ($pgq, $gpq, $gq) = solve_triangle($ra - $rag, 90.0 - $decg, 
                            90.0 - $dec);

  $glong = $pgl - $pgq;
  while ($glong >= 360.0) {
    $glong -= 360; 
  }
  while ($glong < 0.0) {
    $glong += 360.0;
  }
  $glat = 90.0 - $gq;
  while ($glat >= 90.0) {
    $glat -= 180.0;
  }
  while ($glat < -90.0) {
    $glag += 180.0;
  }

  return($glong, $glat);

}




#################################################################
# Solve a spherical triangle given two sides and included angle.
#
# usage:  ($b, $c, $bc)  =  solve_triangle($a, $ab, $ac)
sub solve_triangle {
 


  $a = $_[0];
  $ab = $_[1];
  $ac = $_[2];

  $a *= $DEGTORAD;
  $ab *= $DEGTORAD;
  $ac *= $DEGTORAD;

  $bcnum = sqrt(   pow(cos($a/2.0)*sin($ab-$ac), 2.0 )
                 + pow(sin($a/2.0)*sin($ab+$ac), 2.0 )
                 + pow(sin($ab)*sin($ac)*sin($a), 2.0 ) );
  $bcdenom = cos($ab)*cos($ac) + sin($ab)*sin($ac)*cos($a);
  
  $bc = atan2($bcnum, $bcdenom);

  $b = atan2(sin($ac)*sin($a), cos($ac)*sin($ab) - cos($ab)*sin($ac)*cos($a) );
  $c = atan2(sin($ab)*sin($a), cos($ab)*sin($ac) - cos($ac)*sin($ab)*cos($a) );

  $bc /= $DEGTORAD;
  $b /= $DEGTORAD; 
  $c /= $DEGTORAD;

  return($b, $c, $bc);

}
