#!/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);
}