#!/usr/bin/perl

# short program to make a histogram from a column of ascii data
# Written by Joshua Spodek
# Copyright (c) 1998 Joshua Spodek
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 2
# of the License, or (at your option) any later version.
# 
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
# 
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

$DEFAULT_NBINS = 50;
$column = 1;

# Read the command line arguments
while ($_=$ARGV[0]) {
  shift;
  if (/^-/) {
    /^-nbins/   && ($nbins =$ARGV[0],shift,next);
    /^-xmin/    && ($xmin =$ARGV[0],shift,next);
    /^-xmax/    && ($xmax =$ARGV[0],shift,next);
    /^-ymin/    && ($ymin =$ARGV[0],shift,next);
    /^-ymax/    && ($ymax =$ARGV[0],shift,next);
    /^-binsize/ && ($binsize=$ARGV[0],shift,next);
    /^-column/  && ($column=$ARGV[0],shift,next);
    /^-ycolumn/ && ($ycolumn=$ARGV[0],shift,next);
    /^-v/       && (print( STDERR "histofy version 1.0\n" ),exit);
    /^-h/       && (&usage);
    print STDERR "unknown argument: $_\n";
  }
}

# Set default ycolumn if user has given ymin or ymax but not ycolumn
if (!defined($ycolumn) && (defined($ymin) || defined($ymax))) {$ycolumn=2}

# First load the data
print STDERR "loading data from column number $column ... \n";
$i=0;
undef( @x );
# Skip non-numeric data in front, which is probably comments or qdp commands
while ( <> ) {
  if ($_ !~ /^[\d- ]/) {print STDERR "skipping comment: $_"} else {last}
}

$x[$i] = (split(' ', $_))[$column-1];
$xmax_tmp = $xmin_tmp = $x[$i];
$i++;
while (<>) {
  # Here is where the data is read in.  If outside the bounds of ymax and ymin
  # don't include the point.
  if(!defined($ycolumn)) {
    $x[$i] = (split(' ', $_))[$column-1];
  } else {
    @temp = split(' ', $_);
    $y = $temp[$ycolumn-1];
    next if (defined($ymin) && $y<$ymin);
    next if (defined($ymax) && $y>$ymax);
    #next if ($y ne "no" && $y ne "small" && $y ne "large");
    $x[$i] = $temp[$column-1];
  }
  $xmax_tmp = $x[$i] > $xmax_tmp ? $x[$i] : $xmax_tmp;
  $xmin_tmp = $x[$i] < $xmin_tmp ? $x[$i] : $xmin_tmp;
  $i++;
}
printf STDERR "Read %d points\n", $#x+1;

# If not specified, find limits of distribution
if (!defined($xmax)) {$xmax = $xmax_tmp}
if (!defined($xmin)) {$xmin = $xmin_tmp}
if ($xmax < $xmin) {die "xmax must be greater than xmin\n"}
printf( STDERR "xmin = %5.3f, xmax = %5.3f.  Total width=%5.1f\n", $xmin,
	$xmax, $xmax - $xmin );

# Figure out the binsize or number of bins
if (!defined($nbins) && !defined($binsize)) {$nbins = $DEFAULT_NBINS}
if (defined($nbins) && defined($binsize)) {
  die "don't specify number of bins and binsize at the same time\n";
} else {
  if (!defined($binsize)) {$binsize = ($xmax - $xmin) / $nbins}
  if (!defined($nbins)) {
      $nbins = int(($xmax - $xmin) / $binsize);
      $xmax = $xmin + $nbins*$binsize;
      printf( STDERR "redefining xmax to %5.3f\n", $xmax );
  }
}
printf( STDERR "%5.3f bins of %5.1f width\n", $nbins, $binsize );

# reset the histogram array
for ( $i=0; $i<=$nbins; $i++ ) {
  $bincen[$i] = $xmin + $binsize/2.0 + $i*$binsize;
  $hist[$i] = 0;
}

# do the histogramming
for ( $i=0; $i<=$#x; $i++ ) {
  $bin = int( ($x[$i] - $xmin) / $binsize );
  next if ( $bin < 0 || $bin > $nbins );
  $hist[$bin]++;
}

# output the histogram
for ( $i=0; $i<=$nbins; $i++ ) {
  printf( "%8.6g %6.4f\n", $bincen[$i], $hist[$i] );
}

sub usage {
  print STDERR "\n\tusage: histofy [-column \#] [-nbins \#] ";
  print STDERR "[-xmin \#] [-xmax \#] [-binsize \#]\n\n";
  print STDERR "\t\t-column sets which column to read, default is first.\n";
  print STDERR "\n";
  exit;
}