#!/usr/bin/perl
#
# Make Hammer-Aitoff plots showing novae 
#          as a function of RA, Dec
#          as a function of galactic lat, long



# plot of nova positions (RA, Dec) 
$cmdfile = "make_nova_aitoff_plot.gnu";
$output_file = "nova_ra_dec_plot.ps";

  
  open (CMDFILE, ">$cmdfile") || die("can't open $cmdfile for writing");
  printf CMDFILE "set output '$output_file' \n";
  printf CMDFILE "set term postscript color \n";
  printf CMDFILE "set nogrid \n";
  printf CMDFILE "set key top left \n";
  printf CMDFILE "set xlabel 'Right Ascension (degrees) ' \n";
  printf CMDFILE "set ylabel 'Declination (degrees) ' \n";
  printf CMDFILE "set title 'Novae in the Milky Way (up to 1988) ' \n";
  
  $cmd = "plot [-20:380][-95:95] 'QQ_radec_2173_catalog.dat' ";
  $cmd .=  "      using 1:2 pt 7 ps 0.7 lw 3 t '' ";
  $cmd .= " , ";
  $cmd .= "  'QQ_ha_outline.dat' using 1:2 ";
  $cmd .= "      lt -1 lw 0.5 pt 6 ps 0.1 t ''  ";

  printf CMDFILE "$cmd \n";
  printf CMDFILE "set output \n";
  printf CMDFILE "quit \n";
  close(CMDFILE) ;
  
  $retval = `gnuplot < $cmdfile`;
  printf "retval is ..%s..\n", $retval;

  # create a GIF copy of any Postscript file
    $gif_file = $output_file;
    $gif_file =~ s/.ps/.gif/;
    $cmd = "convert -rotate 90 $output_file $gif_file ";
    printf "cmd is ..$cmd.. \n";
    $retval = `$cmd`;

if (1 == 1) {

# plot of nova positions (glong, glat) 
$cmdfile = "make_nova_aitoff_plot.gnu";
$output_file = "nova_glong_glat_plot.ps";

  
  open (CMDFILE, ">$cmdfile") || die("can't open $cmdfile for writing");
  printf CMDFILE "set output '$output_file' \n";
  printf CMDFILE "set term postscript color \n";
  printf CMDFILE "set nogrid \n";
  printf CMDFILE "set key top left \n";
  printf CMDFILE "set xlabel 'Galactic longitude (degrees) ' \n";
  printf CMDFILE "set ylabel 'Galactic latitude (degrees) ' \n";
  printf CMDFILE "set title 'Novae in the Milky Way (up to 1988) ' \n";
  printf CMDFILE "set xtics (\"180\" 0, \"270\" 90 ,   ";
  printf CMDFILE "           \"0\" 180, \"90\" 270 , \"180\" 360 ) \n";
  
  $cmd = "plot [-20:380][-95:95] 'QQ_glonglat_2173_catalog.dat' ";
  $cmd .=  "      using 1:2 pt 7 ps 0.7 lw 3 t '' ";
  $cmd .= " , ";
  $cmd .= "  'QQ_ha_outline.dat' using 1:2 ";
  $cmd .= "      lt -1 lw 0.5 pt 6 ps 0.1 t ''  ";

  printf CMDFILE "$cmd \n";
  printf CMDFILE "set output \n";
  printf CMDFILE "quit \n";
  close(CMDFILE) ;
  
  $retval = `gnuplot < $cmdfile`;
  printf "retval is ..%s..\n", $retval;

  # create a GIF copy of any Postscript file
    $gif_file = $output_file;
    $gif_file =~ s/.ps/.gif/;
    $cmd = "convert -rotate 90 $output_file $gif_file ";
    printf "cmd is ..$cmd.. \n";
    $retval = `$cmd`;

}



exit 0;


#####################################################################
# Create a temporary data file with a regular grid in (RA, Dec),
#    which has been converted into Hammer-Aitoff coords (x, y).
#    Thus, each line looks like
#
#              x    y
#
# we can then plot this datafile to make a nice outline of the sky.
#
# We break this line into stanzas, one stanza per line on the graph.
# We place two blank lines between each stanza, so that we can tell
# gnuplot to connect the dots within each stanza, but not between stanzas.
#
# Arguments:        datafile ra_spacing dec_spacing
#
# where
#         datafile        is name of file we create
#         ra_spacing      is spacing of RA grid lines in degrees
#         dec_spacing     is spacing of Dec grid lines in degrees
#
# 
sub make_aitoff_outline_datafile {
 
  my($datafile);
  my($ra_spacing, $dec_spacing);
  my($num_points);
  my($ra, $dec);
  my($i, $j, $k);
  my($x, $y);

  $datafile = $_[0];
  $ra_spacing = $_[1];
  $dec_spacing = $_[2];

  # how many points should we place along each grid line?
  $num_points = 50;

  open(DATAFILE, ">$datafile") || 
       die("make_aitoff_outline_datafile: can't open file $datafile");

  # first, we make the lines of constant RA
  for ($ra = 0; $ra <= 360.0; $ra += $ra_spacing) {
    $dec_incr = 180.0/$num_points;
    $dec = -90.0;
    for ($i = 0; $i <= $num_points; $i++) {

      # now, convert RA and Dec to Hammer-Aitoff (x, y)
      ($x, $y) = convert_radec_to_hammeraitoff($ra, $dec);
      printf DATAFILE " %8.4f %8.4f \n", $x, $y;
      $dec += $dec_incr;
    }
    printf DATAFILE "\n\n";
  }
  # next, we make the lines of constant Dec
  for ($dec = -90; $dec <= 90.0; $dec += $dec_spacing) {
    $ra_incr = 360.0/$num_points;
    $ra = 0.0;
    for ($i = 0; $i <= $num_points; $i++) {

      # now, convert RA and Dec to Hammer-Aitoff (x, y)
      ($x, $y) = convert_radec_to_hammeraitoff($ra, $dec);
      printf DATAFILE " %8.4f %8.4f \n", $x, $y;
      $ra += $ra_incr;
    }
    printf DATAFILE "\n\n";
  }

  close(DATAFILE);

}





####################################################################
# Given (RA, Dec), convert to (x, y) in Hammer-Aitoff projection
#
# The two arguments are RA and Dec in decimal degrees,
# and we return a two-element list for x, y.
# We scale (x, y) so that they range from
#
#            0   <   x   <   360      
#          -90   <   y   <    90      
#
# so that labels will correspond to (RA, Dec) in decimal degrees.
#
# MWR 7/22/2006

sub convert_radec_to_hammeraitoff {
  my($ra, $dec);
  my($tra, $tdec);
  my($x, $y, $z);
  my($DEGTORAD);
  my($pi);

  $ra = $_[0];
  $dec = $_[1];

  $pi = 3.14159;
  $DEGTORAD = $pi/180.0;

  # convert RA and Dec into radians
  $tra = ($ra*$DEGTORAD) - $pi; 
  $tdec = ($dec*$DEGTORAD);
 
  # calc auxiliary variable z
  $z = sqrt( 1.0 + cos($tdec)*cos($tra/2.0) );
  
  # calc Hammer-Aitoff normalized coords (x, y)
  $x = ( cos($tdec)*sin($tra/2.0) ) / $z;
  $y = sin($tdec) / $z;

  # now scale the normalized coords so that X ranges over 0 to 360,
  #  and y ranges over -90 to +90.
  $x = ($x*180.0) + 180.0;
  $y *= 90.0;

  return($x, $y);
}


####################################################################
# Given a datafile with (RA, Dec) in columns "ra_col", "dec_col",
#    create a new datafile which has only RA and Dec, converted
#    into Hammer-Aitoff format.
#
# Ignore lines which begin with a comment char "#"
#
# MWR 7/22/2006

make_aitoff_datafile($orig_datafile, $ra_col, $dec_col, $temp_prefix);

sub make_aitoff_datafile  {
  my($orig_datafile, $new_datafile);
  my($temp_prefix);
  my($ra, $dec);
  my($tra, $tdec);
  my($x, $y, $z);
  my($DEGTORAD);
  my($pi);

  $orig_datafile = $_[0];
  $ra_col = $_[1];
  $dec_col = $_[2];
  $temp_prefix = $_[3];

  $new_datafile = sprintf "%s%s", $temp_prefix, $orig_datafile;
  open(NEW_DATAFILE, ">$new_datafile") ||
           die("make_aitoff_datafile: can't open file $new_datafile");
  open(ORIG_DATAFILE, "$orig_datafile");
  while (<ORIG_DATAFILE>) {
    $line = $_;
    if ($line =~ /^#/) {
      next;
    }
    $line = " " . $line;
    chomp($line);
    @words = split(/\s+/, $line);
    $ra = @words[$ra_col];
    $dec = @words[$dec_col];

    # shift the "RA" (or longitude) value so that "RA=0" falls in the very
    #   middle of the graph, instead of at the left edge.
    #   We assume decimal degrees.
    $ra = $ra + 180.0;
    if ($ra >= 360.0) {
      $ra -= 360.0;
    }

    ($x, $y) = convert_radec_to_hammeraitoff($ra, $dec);

    printf NEW_DATAFILE " %lf %lf \n", $x, $y;
  }
  close(ORIG_DATAFILE);

  close(NEW_DATAFILE);

}

