#!/usr/bin/perl
#
# Make diagnostic plots showing the variation in several
#   parameters (sky brightness, FWHM, etc.) as a function
#   of JD during the night.
#
# Place the plots into .ps and .gif files
#
# Assumes 
#         a. that this is running in a "pff" directory
#         b. there is an "addheader.log" file in the .. dir
#         c. there are .coo and .pht files in .. dir
#         d. there are .pht files with headers in this dir
#
# MWR 7/31/2001

$debug = 1;

# all .coo files and .pht files have this base part to their name
$base_name = "../[a-z]*clear-*";

$night = "12-inch clear UT Aug 08 2014";
$minjd = 3000000.0;
$maxjd = 0;
$JDBASE = 2455000.0;

# read info from the "addheader.log" file into arrays
#   and set the globals "minjd" and "maxjd"
$logfile = "../addheader.log";
$numfiles = 0;
open(LOGFILE, "$logfile") || die("can't open file $logfile");
while (<LOGFILE>) {
  @words = split(/\s+/);
  $image[$numfiles] = $words[1];
  $object[$numfiles] = $words[2];
  $filter[$numfiles] = $words[3];
  $exptime[$numfiles] = $words[4];
  $date[$numfiles] = $words[5];
  $ut[$numfiles] = $words[6];

  # calculate the JD for this image
  $cmd = "~/bin/jd ut=$ut[$numfiles] date=$date[$numfiles]";
  if ($debug > 0) {
    printf "cmd is ..$cmd..\n";
  }
  $jd[$numfiles] = `$cmd`;
  chomp($jd[$numfiles]);
  $jd[$numfiles] -= $JDBASE;

  # this checks for invalid dates and ignores them
  if ($jd[$numfiles] == 0.0 - $JDBASE) {
    next;
  }
  printf "file %5d  name %20s  has JD %12.5f \n", $numfiles,
              $image[$numfiles], $jd[$numfiles];

  if ($jd[$numfiles] < $minjd) {
    $minjd = $jd[$numfiles];
  }
  if ($jd[$numfiles] > $maxjd) {
    $maxjd = $jd[$numfiles];
  }

  $numfiles++;
}
close(LOGFILE);

if ($debug > 0) {
  printf "minjd %12.5f maxjd %12.5f \n", $minjd, $maxjd;
}



# make a plot of number of objects vs. JD
# 
# step 1: make a data file with columns of values to plot
$cmdfile = "gnuplot.in";
$datafile = "gnuplot.dat";
open (DATAFILE, ">$datafile") || die("can't open file $datafile for output");
$coo_string = sprintf("%s*.coo", $base_name);
@coos = glob($coo_string);
foreach $coo (@coos) {
  # strip the leading ../ from the name
  $name = $coo;
  $name =~ s/\.\.\///;
  $name =~ s/\.coo//;
  $name =~ s/\_0/\.0/;
  $name =~ s/\_fit/\.fit/;
  printf "next coo is %s\n", $name;
  $index = -1;
  for ($i = 0; $i < $numfiles; $i++) {
    if ($name eq $image[$i]) {
      $index = $i;
	   last;
	 }
  }
  if ($index == -1) {
    printf "error: can't find image ..%s.. for coo ..%s..\n", $name, $coo;
	 next;
  }

  # count the number of lines in the file
  $nobj = 0;
  open(COOFILE, "$coo") || die ("can't open file $coo");
  while (<COOFILE>) {
    $nobj++;
  }
  close(COOFILE);
  printf DATAFILE "%12.5f %5d \n", $jd[$index], $nobj;
}
close(DATAFILE);


# step 2: create the GNUPLOT command file, and execute it
$psfile = "numobj.ps";
$term_options = "postscript color";
open (CMDFILE, ">$cmdfile") || die("can't open $cmdfile for writing");
printf CMDFILE "set output '$psfile' \n";
printf CMDFILE "set term $term_options \n";
printf CMDFILE "set grid \n";
printf CMDFILE "set nokey \n";
printf CMDFILE "set xlabel 'Julian Date - $JDBASE' \n";
printf CMDFILE "set ylabel 'Number of objects in frame' \n";
printf CMDFILE "set title '$night' \n";

$cmd = "plot [$minjd:$maxjd][0:200] '$datafile' using 1:2 pt 6 ps 0.6 lw 3";
printf CMDFILE "$cmd \n";
printf CMDFILE "set output \n";
printf CMDFILE "quit \n";
close(CMDFILE);

$retval = `gnuplot < $cmdfile`;
printf "retval for gnuplot of nobj is ..$retval.. \n";
$giffile = $psfile;
$giffile =~ s/.ps/.png/;
$cmd = "convert -rotate 90 $psfile $giffile";
printf "cmd is ..$cmd.. \n";
$ret = `$cmd`;






# make a plot of FWHM vs. JD
# 
# step 1: make a data file with columns of values to plot
$cmdfile = "gnuplot.in";
$datafile = "gnuplot.dat";
open (DATAFILE, ">$datafile") || die("can't open file $datafile for output");
$coo_string = sprintf("%s*.coo", $base_name);
@coos = glob($coo_string);
foreach $coo (@coos) {
  # strip the leading ../ from the name
  $name = $coo;
  $name =~ s/\.\.\///;
  $name =~ s/\.coo//;
  $name =~ s/\_0/\.0/;
  $name =~ s/\_fit/\.fit/;
  printf "next coo is %s\n", $name;
  $index = -1;
  for ($i = 0; $i < $numfiles; $i++) {
    if ($name eq $image[$i]) {
      $index = $i;
	   last;
	 }
  }
  if ($index == -1) {
    printf "error: can't find image ..%s.. for coo ..%s..\n", $name, $coo;
	 next;
  }

  $fwhm = get_mean($coo, 5);
  printf DATAFILE "%12.5f %6.2f \n", $jd[$index], $fwhm;
}
close(DATAFILE);


# step 2: create the GNUPLOT command file, and execute it
$psfile = "fwhm.ps";
$term_options = "postscript color";
open (CMDFILE, ">$cmdfile") || die("can't open $cmdfile for writing");
printf CMDFILE "set output '$psfile' \n";
printf CMDFILE "set term $term_options \n";
printf CMDFILE "set grid \n";
printf CMDFILE "set nokey \n";
printf CMDFILE "set xlabel 'Julian Date - $JDBASE' \n";
printf CMDFILE "set ylabel 'average FWHM (pix)' \n";
printf CMDFILE "set title '$night' \n";

$cmd = "plot [$minjd:$maxjd][0:6] '$datafile' using 1:2 pt 6 ps 0.6 lw 3";
printf CMDFILE "$cmd \n";
printf CMDFILE "set output \n";
printf CMDFILE "quit \n";
close(CMDFILE);

$retval = `gnuplot < $cmdfile`;
printf "retval for gnuplot of nobj is ..$retval.. \n";
$giffile = $psfile;
$giffile =~ s/.ps/.png/;
$cmd = "convert -rotate 90 $psfile $giffile";
printf "cmd is ..$cmd.. \n";
$ret = `$cmd`;




# make a plot of sky value vs. JD
# 
# step 1: make a data file with columns of values to plot
$cmdfile = "gnuplot.in";
$datafile = "gnuplot.dat";
open (DATAFILE, ">$datafile") || die("can't open file $datafile for output");
$pht_string = sprintf("%s*.pht", $base_name);
@phts = glob($pht_string);
foreach $pht (@phts) {
  # strip the leading ../ from the name
  $name = $pht;
  $name =~ s/\.\.\///;
  $name =~ s/\.pht//;
  $name =~ s/\_0/\.0/;
  $name =~ s/\_fit/\.fit/;
  printf "next pht is %s\n", $name;
  $index = -1;
  for ($i = 0; $i < $numfiles; $i++) {
    if ($name eq $image[$i]) {
      $index = $i;
	   last;
	 }
  }
  if ($index == -1) {
    printf "error: can't find image for pht ..%s..\n", $pht;
	 next;
  }

  $sky = get_mean($pht, 4);
  # we calculate the LOGARITHM base 10 of the sky value
  if ($sky <= 0) {
    $skylog = 0;
  } else {
    $skylog = log($sky)/log(10.0);
  }
# no, don't use logarithm
  $skylog = $sky;
  printf DATAFILE "%12.5f %7.2f %8.4f\n", $jd[$index], $sky, $skylog;
}
close(DATAFILE);


# step 2: create the GNUPLOT command file, and execute it
$psfile = "sky.ps";
$term_options = "postscript color";
open (CMDFILE, ">$cmdfile") || die("can't open $cmdfile for writing");
printf CMDFILE "set output '$psfile' \n";
printf CMDFILE "set term $term_options \n";
printf CMDFILE "set grid \n";
printf CMDFILE "set nokey \n";
printf CMDFILE "set xlabel 'Julian Date - $JDBASE ' \n";
printf CMDFILE "set ylabel 'background sky (DN)' \n";
printf CMDFILE "set title '$night' \n";

#$cmd = "plot [$minjd:$maxjd][1.0:4.4] '$datafile' using 1:3 pt 6 ps 0.6 lw 3";
$cmd = "plot [$minjd:$maxjd][] '$datafile' using 1:3 pt 6 ps 0.6 lw 3";
printf CMDFILE "$cmd \n";
printf CMDFILE "set output \n";
printf CMDFILE "quit \n";
close(CMDFILE);

$retval = `gnuplot < $cmdfile`;
printf "retval for gnuplot of nobj is ..$retval.. \n";
$giffile = $psfile;
$giffile =~ s/.ps/.png/;
$cmd = "convert -rotate 90 $psfile $giffile";
printf "cmd is ..$cmd.. \n";
$ret = `$cmd`;






exit 0;




#########################################################################
# PROCEDURE: get_mean
#
# DESCRIPTION: open the given file, read numbers from the given 
#              column, calculate the mean value.  
#
# RETURNS:
#      the mean value
#
sub get_mean {
  
  my($filename, 
     $column,
	  $value,
	  $sum,
	  $num,
	  $mean,
	  $i);

  $filename = $_[0];
  $column = $_[1];

  $num = 0;
  $sum = 0;
  open(INFILE, $filename) || die ("get_mean: can't open file $filename");
  while (<INFILE>) {
    @words = split(/\s+/);
    $value = $words[$column];
	 if ($debug > 1) {
	   printf "  get_mean: next value is ..%s..\n", $value;
	 }
	 $sum += $value;
	 $num++;
  }
  close(INFILE);

  if ($num > 0) {
    $mean = $sum/$num;
  } else {
    $mean = 0.0;
  }

  return($mean);
}