#!/usr/local/bin/perl
# usage: phaseit.pl filename name period [epoch]
#
# by Michael Richmond <richmond@stupendous.cis.rit.edu>
# and Michael Koppelman <lolife@bitstream.net>

use strict;
use POSIX;

# set type to one of 'postscript' or 'png'
# the postscript files are much prettier and can be converted to PNG or JPEG
my $type = "postscript";
#my $type = "png";

# tell us where your gnuplot is.
my $GNUPLOT = '/sw/bin/gnuplot';


### MAIN ###
my( $file, $star, $period, $epoch ) = @ARGV;

if( ! $ARGV[2] ) {
	die "usage: $0 filename name period [epoch]\n";
}

my $data;	# Array ref of hashrefs with jd, mag and err for each line from file
my $debug = 0;

my $x=0;
open( FILE, $file ) || die "Can't open file $file: $!\n";
while( my $line = <FILE> ) {
	next if( $line =~/^#/ );
	( $data->[$x]->{'jd'}, $data->[$x]->{'mag'}, $data->[$x]->{'err'} ) = split( /\s+/, $line );
	++$x;
}
close( FILE );

my @rtn = &make_phased_graph( $type, $star, $epoch, $period, $data );
if( ! $rtn[0] ) {
	print "Created $rtn[1]\n";
}
else {
	print "Failed.\n";
}
exit(0);
### end MAIN ###


sub make_phased_graph {
	my( $type, $name, $epoch_jd, $period, $data ) = @_;

	my( $term, $plotfile );
	if( $type eq 'postscript' ) {
		$term = 'postscript enhanced mono';
		$plotfile = "$name.$$.eps";
	}
	else {
		$term = 'png color';
		$plotfile = "$name.$$.png";
	}

	my($i);
	my($proc_retval);
	my($plot_cmd_file, $plot_data_file);
	my($retval, $cmd);

	$proc_retval = 1;
	$plot_cmd_file = "$name.$$.gnu";
	$plot_data_file = "$name.$$.dat";


	# if there are no data lines, we return with no graph
	if (scalar($data) < 1 ) {
		printf STDERR "make_phased_graph: no data to plot\n";
		return($proc_retval, undef);
	}

	# if the epoch_jd value is invalid, just use the first value
	if ( !$epoch_jd) {
		$epoch_jd = $data->[0]->{'jd'};
		print STDERR "Epoch defaulting to $epoch_jd\n";
	}

	# if the PERIOD is less than zero, then we can't make a phased graph;
	#		so return with error
	if ($period < 0) {
		printf STDERR "Can't have negative period, mate!\n";
		return($proc_retval, undef);
	}

	# calculate the phase for each observed magnitude, and put data into
	#	 a file with lines that look like this:
	#
	#					JD	 phase	phase-1	phase+1 mag err
	#
	if (create_phased_data($plot_data_file, $epoch_jd, $period, $data) != 0) {
		printf STDERR "make_phased_graph: create_phased_data returns with error";
		return($proc_retval, $plotfile);
	}

	my ($min_mag, $max_mag) = find_extreme_mags( $data );

	# open a file to hold GNUPLOT commands
	#	 and put into it all the command necessary to make a nice plot
	open(CMD_FILE, ">$plot_cmd_file") || die("make_phased_graph: can't open file $plot_cmd_file");
	printf CMD_FILE "set output '%s' \n", $plotfile;
	printf CMD_FILE "set term $term \n";
	printf CMD_FILE "set grid \n";
	printf CMD_FILE "set xlabel 'Phase' \n";
	printf CMD_FILE "set xtics -0.25,0.25,1.25 \n";
	printf CMD_FILE "set ylabel 'Magnitude' \n";
	printf CMD_FILE "set title '%s (Ephemeris = %.5f + %.5fE)\n", $name, $epoch_jd, $period;
	printf CMD_FILE "plot [-0.25:1.25][%f:%f] '%s' using 2:5:6 notitle with yerrorbars lt 1 pt 7 ", $max_mag, $min_mag, $plot_data_file;
	printf CMD_FILE " ,	'%s' using 3:5:6 notitle with yerrorbars lt 1 pt 7 ", $plot_data_file;
	printf CMD_FILE " ,	'%s' using 4:5:6 notitle with yerrorbars lt 1 pt 7 ", $plot_data_file;
	printf CMD_FILE "\n";
	close(CMD_FILE);
															 
	# execute the "gnuplot" program on the commands to create
	#		a file holding the graph
	$cmd = "$GNUPLOT $plot_cmd_file";
	$retval = `$cmd`;
	if ($? != 0) {
		printf STDERR "make_phased_graph: gnuplot returns with error \n";
		return($proc_retval, $plotfile);
	}

	# all is done
	$proc_retval = 0;

	return($proc_retval, $plotfile);
}

sub create_phased_data {
	my( $datafile_name, $epoch_jd, $period, $data ) = @_;

	my($i);
	my($retval);
	my($line, @words);
	my($jd, $mag, $err);
	my($phase, $x);
	my($vmag);

	# by default, we haven't succeeded
	$retval = 1;

	if ($debug > 0) {
		printf "create_phased_data: file ..%s.. epoch ..%12.5f.. period %12.5f \n", $datafile_name, $epoch_jd, $period;
	}
	# sanity checks
	if ($epoch_jd < 0) {
		printf STDERR "create_phased_data: bad epoch %12.5f \n", $epoch_jd;
		return($retval);
	}
	if ($period <= 0) {
		printf STDERR "create_phased_data: bad period 12.5f \n", $period;
		return($retval);
	}
	
	open(PHASED_DATA, ">$datafile_name") || die("create_phased_data: can't open $datafile_name");

	foreach my $d( @{$data} ) {
		$jd	= $d->{'jd'};
		$mag = $d->{'mag'};
		$err = $d->{'err'};

		$x = ($jd - $epoch_jd)/$period;
		$phase = $x - floor($x);
		if ($debug > 0) {
			 printf "	jd ..%s..	mag ..%s..", $jd, $mag;
			 printf " phase %8.3f \n", $phase;
		}

		# now we print a line into the output datafile
		printf PHASED_DATA " %12.5f %7.3f %7.3f %7.3f %7.3f %7.3f\n", $jd, $phase, $phase-1.0, $phase+1.0, $mag, $err;

	}
	close(PHASED_DATA);

	# if we get here, all is well!
	$retval = 0;
	return($retval);
}

sub find_extreme_mags {
	my( $data ) = @_;

	my($i);
	my($default_min_mag, $default_max_mag);
	my($bad_min_mag, $bad_max_mag);
	my($min_mag, $max_mag); 
	my($mag, $mag_column);
	my($line, @words);
	my($num_lines);

	# these will signal no real data if they persist
	$bad_min_mag = 100;
	$bad_max_mag = -100;

	# these are sensible defaults
	$default_min_mag = 6.0;
	$default_max_mag = 15.0;

	# find the min, max values in the data
	$min_mag = $bad_min_mag;
	$max_mag = $bad_max_mag;
	foreach my $d ( @{$data} ) {
		$mag = $d->{'mag'};
		if ($mag < $min_mag) {
			$min_mag = $mag;
		}
		if ($mag > $max_mag) {	
			$max_mag = $mag;
		}
	}
	if ($min_mag == $bad_min_mag) {
		$min_mag = $default_min_mag;
	}
	if ($max_mag == $bad_max_mag) {
		$max_mag = $default_max_mag;
	}

	return($min_mag, $max_mag);
}