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

$cmdfile = "makeplot.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;



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.
make_aitoff_datafile("snlist_noz.dat", 5, 6, $temp_prefix);
#make_aitoff_datafile("my_2173_catalog.dat", 5, 6, $temp_prefix);


exit(0);



# plot of V-band scatter(s) versus mean mag

$output_file = "sigma_vs_mag_v.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 grid \n";
  printf CMDFILE "set key top left \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";
  
  $cmd = "plot [5:16][-0.22:0.50] 'subset.dat' using 8:11 ps 0.7 t 'raw' ";
  $cmd .= " , ";
  $cmd .= " 'subset.dat' using 8:(\$12-0.20) ps 0.7 t 'ensemble' ";
  $cmd .= " , ";
  $cmd .= " 0.0 lt -1 t ''  ";
  $cmd .= " , ";
  $cmd .= " -0.20 lt -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];

if (1 == 1) {
    # shift the "RA" (or longitude) value by 180 degrees
    #   so that 
    #   We assume decimal degrees.
    $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);

}

