#!/usr/bin/perl -w

use strict;
use Getopt::Long;
use POSIX qw(log10 ceil);
use PGPLOT;
use Astro::Time;

my $xcol = 1;
my @ycol = ();
my @label = ();
my $skip = 0;
my $header = undef;
my $logx = 0;
my $logy = 0;
my $size = 1;
my $symbol = 5;
my $line = 0;
my $dashed = 0;
my $width = undef;
my $time = 0;
my $sectime = 0;
my $histo = 0;
my $cumul = 0;
my $panel = 0;
my $delta = 0;
my $scale = 1;
my $xscale = 1;
my $xoffset = undef;
my $yoffset = 0;
my $ramp = 0;
my $nox = 0;
my $xxmin = undef;
my $xxmax = undef;
my $yymin = undef;
my $yymax = undef;
my $smooth = undef;
my @mark = ();
my @xmark = ();
my $dev = '/xs';
my $title = undef;
my $xlab = undef;
my $ylab = undef;
my $colour = 1;
my $nocolour = 0;
my $invert = 0;
my $help = 0;
my $xdiv = undef;
my $ydiv = undef;
my $maxfilt = undef;
my $minfilt = undef;
my $xfilt = undef;

GetOptions('xaxis=i'=>\$xcol, 'x=i'=>\$xcol, 'ycol=s'=>\@ycol,
	   'col=s'=>\@ycol, 'skip=i'=>\$skip, 'header'=>\$header,
	   'logx'=>\$logx, 'logy'=>\$logy, 'xmin=s'=>\$xxmin,
	   'xmax=s'=>\$xxmax, 'dev=s'=>\$dev, 'ymin=s'=>\$yymin,
	   'ymax=s'=>\$yymax, 'line'=>\$line, 'time'=>\$time,
	   'sectime'=>\$sectime, 'nocolour'=>\$nocolour,
	   'nocolor'=>\$nocolour, 'help'=>\$help,
	   'kumulative'=>\$cumul, 'scale=f'=>\$scale, 'nox'=>\$nox,
	   'xscale=f'=>\$xscale, 'smooth=i'=>\$smooth, 'xdiv=f'=>\$xdiv,
	   'title=s'=>\$title, 'xlab=s'=>\$xlab, 'ylab=s'=>\$ylab,
	   'mark=f'=>\@mark, 'xmark=f'=>\@xmark, 'xoffset=s'=>\$xoffset, 'label=s'=>\@label,
	   'dashed'=>\$dashed, 'div=f'=>\$ydiv, 'ydiv=f'=>\$ydiv,
	   'maxfilt=f'=>\$maxfilt, 'minfilt=f'=>\$minfilt, 'xfilt=f'=>\$xfilt, 'histo'=>\$histo,
	   'size=f'=>\$size, 'symbol=i'=>\$symbol, 'delta'=>\$delta,
	   'yoffset=s'=>\$yoffset, 'panel'=>\$panel, 'invert'=>\$invert,
	   'width=f'=>\$width);

$colour = 0 if $nocolour;

if ($help) {
print<<EOF
 Usage: quickplot.pl [options] <file>
 
 Options:
   xaxis <col>     Column with x axis data (default first column)
   col/ycol <col>  Column with y axis data (default second column). 
                   Can be specified multiple times to plot more than
 	           one column.
   dev             Pgplot device to use
   header          Skip the first "skip" rows
   skip <n>        Number of rows to skip if -header options used (default 1)
   logx            Take log10 of the x-axis data
   logy            Take log10 of the y-axis data
   xmin <min>      Force minimum x axis (default autoscale)
   xmax <max>      Force maximum x axis (default autoscale)
   ymin <min>      Force minimum y axis (default autoscale)
   ymax <max>      Force maximum x axis (default autoscale)
   maxfilt <val>   Ignore any values > val
   minfilt <val>   Ignore any values < val
   xfilt <val>     Ignore any lines with xvalule > val
   line            Plot lines rather than points
   dashed          Plot dashed lines
   nocolour        Don't use colour
   kumulative      Plot cumulative plot(?)
   delta           Plot change in values
   time            Interprete x axis as formatted time 
   sectime         Plot x axis in time, assuming values = time
   scale           Scale all y values by this amount
   xscale          Scale x axis values by this amount
   xoffset         Add this number to all "x" values, after scaling
   offset          Add this number to all "y" values, after scaling
   smooth <n>      Plot (as a line) an average of n values
   mark <value>    Draw a horizontal line at the given value. Multiple
                   values allowed (pass -mark multiple times)
   xmark <value>   Draw a vertical line at the given value. Multiple
                   values allowed (pass -xmark multiple times)
   title <value>   Title for plot
   xlab            Label for X-axis
   ylab            Label for Y-axis
   help            This message

EOF
  ;
exit 0;
}

my @colours = (2,3,4,5,6,7);
my @styles = (1, 2, 3, 4, 5);
my $ncol = scalar(@colours);
my $nstyle = scalar(@styles);

$panel = 0 if $histo;

$nox=1 if ($cumul || $histo);
if ($nox) {
  $xcol = undef;
} else {
  $xcol--;
}
my $i=0;
while ($i<@ycol) {
  if ($ycol[$i] =~ /^(\d+):(\d+)\+(\d+)$/) {
    my @new = ();
    for (my $j=$1; $j<=$2; $j+=$3) {push @new, $j};
    splice @ycol, $i, 1, @new;
    $i += scalar(@new);
  } elsif ($ycol[$i] =~ /^(\d+):(\d+)$/) {
    splice @ycol, $i, 1, ($1..$2);
    $i += abs($2-$1)+1;
  } elsif ($ycol[$i] =~ /^\d+$/) {
    $i++;
  } else {
    die "Did not understand column $ycol[$i]\n";
  }
}
foreach (@ycol) {
  $_--;
}
@ycol = (0) if (@ycol==0);

if ($xdiv) {
  $xscale /= $xdiv;
}

if ($ydiv) {
  $scale /= $ydiv;
}

if ($header && $skip==0) {
  $skip=1;
}

if ($time) {
  if (defined $xxmin ) {
    if ($xxmin =~ /^(\d+)\/(\S+)$/) {
      $xxmin = str2turn($2, 'H')  + $1;
      $xxmin *= 24 * 60 * 60; # day to seconds
    } else {
      $xxmin = str2turn($xxmin, 'H') * 24 * 60 * 60; # String to seconds
    }
  }

  if (defined $xxmax ) {
    if ($xxmax =~ /^(\d+)\/(\S+)$/) {
      $xxmax = str2turn($2, 'H')  + $1;
      $xxmax *= 24 * 60 * 60; # day to seconds
    } else {
      $xxmax = str2turn($xxmax, 'H') * 24 * 60 * 60; # String to seconds
    }
  }
} 


if (!defined $xoffset) {
  $xoffset = 0;
} elsif ($time || $sectime) {
  $xoffset = str2turn($xoffset, 'H') * 24*60*60;
}

if (!defined $yoffset) {
  $yoffset = 0;
}

# Skip header section
while ($skip>0 && defined(<>)) {
  $skip--;
}

die "EOF reached\n" if ($skip>0);

my (@x, @y, $x, $y, @ally);
my @line = ();
my $n = 0;
my $yramp=0;
while (<>) {
  # Clean up line
  s/^\s+//;
  s/\s+$//;
  next if ($_ eq '');

  @line = split;
  my $badval = 0;
  if (defined $xcol) {
    if ($time) {
      $x = $line[$xcol];
    } else {
      $x = ($line[$xcol]*$xscale) + $xoffset;
    }
  } else {
    $x = ($n+1)*$xscale + $xoffset;
  }
  if (defined $xfilt && $x > $xfilt) {
    $badval = 1;
    last;
  }
  @y = ();
  foreach (@ycol) {
    if (!defined($line[$_])) {
      $badval = 1;
      last;
    }
    $y = $line[$_]*$scale - $yoffset - $yramp;
    
    if (defined $maxfilt && $y>$maxfilt) {
      $badval = 1;
      last;
    }
    if (defined $minfilt && $y<$minfilt) {
      $badval = 1;
      last;
    }
    push @y, $y;
  }
  $yramp+=$ramp;

  if (defined($x) && !$badval) {
    if ($time) {
      if ($x =~ /^(\d+)\/(\S+)$/) {
	$x = str2turn($2, 'H')  + $1;
	$x *= 24 * 60 * 60; # day to seconds
	$x += $xoffset;
      } else {
	$x = str2turn($x, 'H') * 24 * 60 * 60; # String to seconds
      }
    }

    push @x, $x;
    for (my $i=0; $i<@y; $i++) {
      push @{$ally[$i]}, $y[$i];
    }
    $n++;
  } else {
    warn "Skipping $_\n";
  }
}

if ($logx) {
  foreach (@x) {
    $_ = log10($_);
  }
}

if ($logy) {
  foreach (@ally) {
    foreach (@$_) {
      $_ = log10($_);
    }
  }
}


die "\nNot enough points read to plot\n" if (scalar @x < 2);

$n = (@x > @ally) ? @x : @ally;


if ($cumul) {
  foreach (@x) {
    $_ = $_/$n*100;
  }

  foreach (@ally) {
    @$_ = sort {$a <=> $b} @$_;
  }
}

if ($delta) {
  foreach (@ally) {
    for (my $i=1; $i<@$_; $i++) {
      $_->[$i-1] = $_->[$i] - $_->[$i-1];
    }
    pop @$_;
  }
  shift @x;
  $n--;
}

my $xmin = $x[0];
my $xmax = $x[0];
my $ymin = $ally[0][0];
my $ymax = $ally[0][0];

foreach (@x) {
  if ($_ > $xmax) {
    $xmax = $_;
  } elsif ($_ < $xmin) {
    $xmin = $_;
  }
}

foreach (@ally) {
  foreach (@$_) {
    if ($_ > $ymax) {
      $ymax = $_;
    } elsif ($_ < $ymin) {
      $ymin = $_;
    }
  }
}

my $bound;


$xmin = $xxmin if (defined $xxmin);
$xmax = $xxmax if (defined $xxmax);
$ymin = $yymin if (defined $yymin);
$ymax = $yymax if (defined $yymax);

$bound = ($xmax - $xmin)*0.01;
$xmax += $bound if (! defined $xxmax);
$xmin -= $bound if (! defined $xxmin);

$bound = ($ymax - $ymin)*0.05;
$ymax += $bound if (! defined $yymax);
$ymin -= $bound if (! defined $yymin);

if ($panel) {
  my $nplot = scalar(@ally);
  my $nx = ceil(sqrt($nplot));
  my $ny = ceil($nplot/$nx);
  pgbeg(0,$dev,$nx,$ny);
} else {
pgbeg(0,$dev,1,1);
}

if (!$panel) {
if ($histo) {
    pgenv($ymin, $ymax, 0, 50, 0, 0);
} elsif ($time || $sectime) {
  pgpage();
  if ($invert) {
    pgswin($xmax, $xmin, $ymin, $ymax);
  } else {
    pgswin($xmin, $xmax, $ymin, $ymax);
  }
  pgtbox('BCNZHT', 0.0, 0, 'BCNTS', 0.0, 0);
  if (!defined $xlab) {
    $xlab = 'Time';
  }
} else {
  if ($invert) {
    pgenv($xmax, $xmin, $ymin, $ymax, 0, 0);
  } else {
    pgenv($xmin, $xmax, $ymin, $ymax, 0, 0);
  }
}

if (defined $title || defined $xlab || defined $ylab) {
  $title = '' if (! defined $title);
  $xlab = '' if (! defined $xlab);
  $ylab = '' if (! defined $ylab);

  pglab($xlab, $ylab, $title);
}
}

$i = 0;
foreach (@ally) {

  if ($panel) {
    pgsci(1);
    pgenv($xmin, $xmax, $ymin, $ymax, 0, 0);

    if (@xmark) {
      foreach (@xmark) {
	pgsci($colours[3]);
	pgmove($_, $ymin);
	pgdraw($_, $ymax);
      }
    }
  }

  pgsci(1);

  if ($colour) {
    pgsci($colours[$i%$ncol]);
  }

  if ($histo) {
    pghist($n, $_, 0, $ymax, 20, 1);
  } elsif ($line) {
    pgslw($width) if (defined $width);
    if ($dashed) {
      pgsls($styles[$i%$nstyle]);
    }
    pgline($n, \@x, $_);
  } else {
    pgsch($size);
    pgpt($n, \@x, $_, $symbol);
    pgsch(1);
  }
  $i++;
}

if ($smooth) {
  foreach (@ally) {
    my @smoothx;
    my @smoothy;
    $i=0;
    my $m=0;
    while ($i<$n) {
      my $sumx = 0;
      my $sumy = 0;
      my $j;
      for ($j=0; $j<$smooth && $i<$n; $j++) {
	$sumx+= $x[$i];
	$sumy+= $_->[$i];
	$i++;
      }
      push @smoothx, $sumx/$j;
      push @smoothy, $sumy/$j;
      $m++;
    }
    pgsci($colours[3]);
    pgline($m, \@smoothx, \@smoothy);
  }
}

if (@mark) {
  foreach (@mark) {
    pgsci($colours[3]);
    pgmove($xmin, $_);
    pgdraw($xmax, $_);
  }
}


if (@label) {
  if (scalar(@label>@ycol)) {
    warn "More labels passed than columns to plot!\n";
    # Do something about it...
  }
  $i = 0;
  pgsch(1.2);
  foreach (@label) {
    if ($colour) {
      pgsci($colours[$i%$ncol]);
    }
    pgmtxt('T', -2-$i*1.02, 0.98, 1, $_);

    $i++;
  }
}

pgend();