#!/usr/bin/perl
#
# Practice making Hammer-Aitoff plots
#

# modify these to suit the datafile we're going to plot
$infile = "euve_allsrc2.txt";
$ra_col = 3;
$dec_col = 4;
$outfile = "makeplots_out.ps";


$cmdfile = "makeplots.gnu";
$temp_outline_file = "QQ_ha_outline.dat";
# add this prefix to filenames when we create Hammer-Aitoff versions
$temp_prefix = "QQ_";
$ra_spacing = 30.0;
$dec_spacing = 30.0;

# should we make a graph using the horizontal range 0 to 360 degrees,
#    or -180 to +180 degrees?
$do_neg180_to_pos180 = 0;


make_aitoff_outline_datafile($temp_outline_file, $ra_spacing, $dec_spacing);

# 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.
$temp_datafile = sprintf "%s%s", $temp_prefix, $infile;
make_aitoff_datafile($infile, $ra_col, $dec_col, $temp_prefix);


# and now create a nice plot in Aitoff form
#    which has a grid for RA and Dec lines
#    and small dots for the objects from the given datafile
make_aitoff_plot($temp_datafile, $outfile, $temp_prefix);

exit(0);



#####################################################################
# Create a nice graph showing some (RA, Dec) positions in Hammer-Aitoff
#    projection.  
#
# Arguments:        datafile outfile  temp_prefix 
#
# where
#         datafile        is name of file with converted Aitoff coords
#         outfile         is name of Postscript/PNG files we create
#         temp_prefix     is prefix we use when making temp files 
#
# 
sub make_aitoff_plot {
 
  my($datafile);
  my($outfile);
  my($temp_prefix);
  my($cmd);
  my($retval, $png_file);

  $infile = $_[0];
  $outfile = $_[1];
  $temp_prefix = $_[2];

  
  open (CMDFILE, ">$cmdfile") || die("can't open $cmdfile for writing");
  printf CMDFILE "set output '$outfile' \n";
  printf CMDFILE "set term postscript color 'Helvetica,18' size 12,8 \n";
  printf CMDFILE "set nokey \n";
#  printf CMDFILE "set xlabel 'Mean magnitude' \n";
#  printf CMDFILE "set ylabel 'Internal scatter in magnitude' \n";
#  printf CMDFILE "set title 'Random subset of Mark IV, V-band' \n";
  
  if ($do_neg180_to_pos180 == 0) {
    $cmd = "plot [-10:370][-100:100] '$temp_outline_file' using 1:2 ";
  }
  else {
    $cmd = "plot [-190:190][-100:100] '$temp_outline_file' using 1:2 ";
  }
  $cmd .= "      with lines lc -1 lw 0.1 t '' ";
  $cmd .= " , ";
  $cmd .= " '$infile' using 1:2 with points ";
  $cmd .= "      pt 5 ps 0.3 lc 1 lw 2 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 PNG copy of any Postscript file
    $png_file = $outfile;
    $png_file =~ s/.ps/.png/;
    $cmd = "convert -flatten -background white -rotate 90 $outfile $png_file ";
    printf "cmd is ..$cmd.. \n";
    $retval = `$cmd`;

  return(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;

  # and now, if the user has asked for output in the range of
  # -180 to +180 degrees for the horizontal direction,  
  # change all "x" values so that the range changes 
  #     from     0   to   360  degrees
  #     to    -180   to  +180  degrees
  #
  if ($do_neg180_to_pos180 == 1) {
    if ($x > 180) {
      $x -= 360.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];

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

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

  close(NEW_DATAFILE);

}

