#!/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;
$xcolumn = 1;
$ycolumn = 2;
$conlin = 1;   # linear spacing for contours is default
$num_contours = 8;
$twod  = 0;
$fits  = 0;

# Read the command line arguments
while ($_=$ARGV[0]) {
  shift;
  if (/^-/) {
    /^-xnbins/   && ($xnbins =$ARGV[0],shift,next);
    /^-ynbins/   && ($ynbins =$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);
    /^-xbinsize/ && ($xbinsize=$ARGV[0],shift,next);
    /^-ybinsize/ && ($ybinsize=$ARGV[0],shift,next);
    /^-xcolumn/  && ($xcolumn=$ARGV[0],shift,next);
    /^-ycolumn/  && ($ycolumn=$ARGV[0],shift,next);
    /^-conlin/   && ($conlin=1,$num_contours=$ARGV[0],shift,next);
    /^-conlog/   && ($conlin=0,$num_contours=$ARGV[0],shift,next);
    /^-2d/       && ($twod=1,next);
    /^-fits/     && ($fits=1,next);
    /^-v/       && (print( STDERR "histofy version 1.0\n" ),exit);
    /^-h/        && (&usage);
    print STDERR "unknown argument: $_\n";
  }
}

undef( @x );
undef( @y );
# Skip non-numeric data in front, which is probably comments or qdp commands
while ( <> ) {
    if ($_ !~ /^[\d- ]/) {print STDERR "skipping comment: $_"} else {last}
}

# First load the data
print STDERR "loading x data from column number $xcolumn\n";
if ($twod == 1 ) {print STDERR "loading y data from column number $ycolumn\n"}
$i=0;
$x[$i] = (split(' ', $_))[$xcolumn-1];
$xmax_tmp = $xmin_tmp = $x[$i];
if ($twod == 1) {
    $y[$i] = (split(' ', $_))[$ycolumn-1];
    $ymax_tmp = $ymin_tmp = $y[$i];
}
$i++;
while (<>) {
    $x[$i] = (split(' ', $_))[$xcolumn-1];
    $xmax_tmp = $x[$i] > $xmax_tmp ? $x[$i] : $xmax_tmp;
    $xmin_tmp = $x[$i] < $xmin_tmp ? $x[$i] : $xmin_tmp;
    if ($twod == 1) {
	$y[$i] = (split(' ', $_))[$ycolumn-1];
	$ymax_tmp = $y[$i] > $ymax_tmp ? $y[$i] : $ymax_tmp;
	$ymin_tmp = $y[$i] < $ymin_tmp ? $y[$i] : $ymin_tmp;
    }
    $i++;
}
printf STDERR "Read %d points\n", $#x+1;

#
# Set the limits of the distribution and binsize for x (y later)
#

# 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($xnbins) && !defined($xbinsize)) {$xnbins = $DEFAULT_NBINS}
if (defined($xnbins) && defined($xbinsize)) {
    die "don't specify number of bins and binsize at the same time\n";
} else {
    if (!defined($xbinsize)) {$xbinsize = ($xmax - $xmin) / $xnbins}
    if (!defined($xnbins)) {
	$xnbins = int(($xmax - $xmin) / $xbinsize);
	$xmax = $xmin + $xnbins*$xbinsize;
	printf( STDERR "Redefining xmax to %5.3f\n", $xmax );
    }
}
printf( STDERR "$xnbins bins of %5.1f width in x\n", $xbinsize );

# reset the x axis array
for ( $i=0; $i<=$xnbins; $i++ ) {
    $xbincen[$i] = $xmin + $xbinsize/2.0 + $i*$xbinsize;
}

#
# Set the limits and binsize for y (x above)
#

if ($twod == 1) {
    # If not specified, find the limits of the distribution
    if (!defined($ymax)) {$ymax = $ymax_tmp}
    if (!defined($ymin)) {$ymin = $ymin_tmp}
    if ($ymax < $ymin) {die "xmax must be greater than xmin\n"}
    printf( STDERR "ymin = %5.3f, ymax = %5.3f.  Total height =%5.1f\n",
	   $ymin, $ymax, $ymax - $ymin );

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

    # Reset the y axis array
    for ( $i=0; $i<=$ynbins; $i++ ) {
	$ybincen[$i] = $ymin + $ybinsize/2.0 + $i*$ybinsize;
    }
}

# Reset the histogram array
if ($twod == 0) {
    undef( @hist );
    for ( $i=0; $i<=$xnbins; $i++ ) {$hist[$i]=0}
} else {
    undef( %hist );
    for ( $i=0; $i<=$xnbins; $i++ ) {
	for ( $j=0; $j<$ynbins; $j++ ) {$hist{$i,$j}=0}
    }
}

# Do the histogramming
for ( $i=0; $i<=$#x; $i++ ) {
  $xbin = int( ($x[$i] - $xmin) / $xbinsize );
  next if ( $xbin < 0 || $xbin > $xnbins );

  if ($twod == 1) {
    $ybin = int( ($y[$i] - $ymin) / $ybinsize );
    next if ( $ybin < 0 || $ybin > $ynbins );
  }

  # Increment a 1d array for 1d histograms and a 2d array for 2d histograms
  ($twod == 0) ? $hist[$xbin]++ : $hist{$xbin,$ybin}++;
}

# output the histogram
if ($twod == 0) {
  for ( $i=0; $i<$xnbins; $i++ ) {
    printf( "%5.3f %5d\n", $xbincen[$i], $hist[$i] );
  }
} else {
  if ($fits == 1) {
    print STDERR "Writing output in fits format\n";
    &make_fits_header( $xnbins, $ynbins, "f", "Joshua Spodek" );
    &write_to_fits( $xnbins, $ynbins, "f" );
    exit;
  } else {
    #printf( "%d %d\n", $xnbins, $ynbins );
    #for ( $i=0; $i<=$xnbins; $i++ ) {printf( "%5.2f ", $xbincen[$i] )}
    #print "\n";
    #for ( $i=0; $i<=$ynbins; $i++ ) {printf( "%5.2f ", $ybincen[$i] )}
    #print "\n";

    # Find the bins with the most and least number of counts
    $max_counts=0;
    $min_counts=99999;
    for ( $i=0; $i<$xnbins; $i++ ) {
      for ( $j=0; $j<$ynbins; $j++ ) {
	if ($hist{$i,$j} > $max_counts) {$max_counts = $hist{$i,$j}}
	if ($hist{$i,$j} < $min_counts) {$min_counts = $hist{$i,$j}}
      }
    }

    printf( "col of 1..999\n" );
    printf( "xaxis lin %d %5.3f\n", $xmin, $xbinsize );
    printf( "yaxis lin %d %5.3f\n", $ymin, $ybinsize );
    printf( "! %5.3f %5.3f %d xbins\n", $xmin, $xmax, $xnbins );
    printf( "! %5.3f %5.3f %d ybins\n", $ymin, $ymax, $ynbins );
    for ( $i=0; $i<$xnbins; $i++ ) {
      for ( $j=0; $j<$ynbins; $j++ ) {
	printf( "%5d", $hist{$i,$j} );
      }
      print "\n";
    }

    printf( "con lev" );
    # Write the contour levels according to user preference or defaults
    if ($conlin == 1) {
      # linear contour spacing
      $interval = int( ($max_counts - $min_counts)/$num_contours );
      for ( $i=$min_counts+$interval; $i<$max_counts; $i+=$interval ) {
	printf " %d", $i;	    
      }
    } else {
      # logarithmic contour spacing
      if ($min_counts == 0) {
	$interval = $max_counts**(1.0/$num_contours);
	for ( $i=1; $i<$max_counts; $i*=$interval ) {
	  printf " %5.3f", $i;
	}
      } else {
	$interval = ($max_counts/$min_counts)**(1.0/$num_contours);
	for ( $i=$min_counts; $i<$max_counts; $i*=$interval ) {
	  printf " %5.3f", $i;
	}
      }
    }
    print "\n";
  }
}

sub make_fits_header {
    local( $nx, $ny, $pixtype, $creator ) = @_;
    local( $foo, $str, $bitpix, @cards );

    if ( $pixtype=="f" ) { $bitpix = -32 }

    # fits file dimensions are in $nx and $ny.
    undef(@cards);
    $creator="Fits file made by 2dhistofy";
    $foo=pack("A8A2A10A60","SIMPLE",  "= ","T","/ $creator");push(@cards,$foo);
    $foo=pack("A8A2A10A60","BITPIX",  "= ",$bitpix,"/");     push(@cards,$foo);
    $foo=pack("A8A2A10A60","NAXIS",   "= ","2","/");         push(@cards,$foo);
    $foo=pack("A8A2A10A60","NAXIS1",  "= ",$nx,"/");         push(@cards,$foo);
    $foo=pack("A8A2A10A60","NAXIS2",  "= ",$ny,"/");         push(@cards,$foo);
    $foo=pack("A80","END");                                  push(@cards,$foo);
    $foo=pack("A80"," ");
    $i=0;
    foreach $str ( @cards ) { $i++; print STDOUT $str; }
    while ($i % 36 != 0) {$i++;print STDOUT $foo;}
}

sub write_to_fits {
  local( $naxis1, $naxis2, $pixtype ) = @_;
  local( $linenumber, $i, $pixel_index, $bitpix, $pixperblock );

  if ( $pixtype=="f" ) { $bitpix = -32 }

  for ( $i=0; $i<$naxis2; $i++ ) {
    for ( $j=0; $j<$naxis1; $j++ ) {
      print STDOUT pack( "$pixtype", $hist{$j,$i} );
      printf( STDERR "%5d", $hist{$j,$i} );
    }
    printf( STDERR "\n" );
  }

  # fits files should have an integer multiple of 2880 bytes.
  # pad the end with zeroes if necessary
  $histsize = $naxis1 * $naxis2;

  $pixperblock = &abs( int( 2880 * 8 / $bitpix ) );

  for ( $i=1; ($histsize+$i-1)%$pixperblock != 0; $i++ ) {
    print STDOUT pack( "$pixtype", 0.0 );
  }
}
 
sub abs {
    # return the absolute value of a number

    local( $x ) = @_;
    return( ($x>0) ? $x : -$x );
}

sub usage {
    print STDERR "\n\tusage: histofy [-xcolumn \#] [-xnbins \#] ";
    print STDERR "[-xmin \#] [-xmax \#] [-xbinsize \#]\n\n";
    print STDERR "\t\t-xcolumn sets which column to read, default is first.\n";
    print STDERR "All the options can also be set for y in a 2d histogram\n";
    print STDERR "\n";
    exit;
}