#!/usr/bin/perl  
require 5.002;
use Getopt::Long;
$opt_line=".*";
$opt_c=1;
$opt_w=-1;
$opt_nbin=40;
$opt_not='^$'; #  '
$opt_s=0;
$opt_norm=0;
GetOptions(qw(help v line:s cum log d s e a norm:i not:s dx:f c:i w:i xmax:f xmin:f n nbin:i));

if ( $opt_help ) {
    print "$0 makes simple histograms \n";
    print "options:  -line line-match-pattern\n";
    print "          -not pattern   exclude pattern\n";
    print "          -dx bin-width calculates nbin \n";
    print "          -c   value column in data-file to do histo on \n";
    print "          -w   weight column in data-file to do histo on \n";
    print "          -s   calculates spread instead of error of mean \n";
    print "          -a   auto binning \n";
    print "          -log         \n";
    print "          -xmax xmax    otherwise dynamic    \n";
    print "          -xmin xmin    otherwise dynamic    \n";
    print "          -nbin number-of-bins   \n";
    print "\n";
    exit;
}

$xmin=1e30;
$xmax=-1e30;
$n=0;
while (<>) {
    if ( /$opt_line/ && $opt_d) {
    	print STDERR $_;
    }
    if ( /$opt_line/ && ! /$opt_not/ ) {
	(@line)=split;
	if ($#line >= $opt_c-1) {
	    $val[$n]=$line[$opt_c-1];
         
            if ( $opt_w > 0 ) { 
		if ( $opt_w <= $#line+1 ) { 
		    $weight[$n]=$line[$opt_w-1];
		}else {
		    $weight[$n]=0;
		}    
            } else {
		$weight[$n]=1.0;
	    }
            print " val:$val[$n] weight:$weight[$n] -- weight_col:$opt_w colums:$#line+1\n" if $opt_v ;
	    $xmax=$val[$n] if $xmax<$val[$n] ;
	    $xmin=$val[$n] if $xmin>$val[$n] ;
	    $n++;
	}
    }
}

if( $opt_a ) {
    @val2=sort {$a <=> $b} @val;
    $nb=0;
    for ($i=1;$i<$n;$i++){
	if ( $val2[$i] != $val2[$i-1] ) {$nb++;$val2[$nb]=$val2[$i]};
    }
    $nb++;
    # try lin 
    $go=1;
    for ($i=2;$i<$nb;$i++){
       if( abs(($val2[$i-1]-$val2[$i-2])-($val2[$i]-$val2[$i-1]))>0.01) {
	   $go=0;
       }
    }
    if ( $go ) {
	$opt_nbin=$nb;
	$dx=($val2[$nb-1]-$val2[0])/($nb-1);
	$xmin=$val2[0]-$dx/2;
	$xmax=$val2[$nb-1]+$dx/2;
	print STDERR "auto-lin $val2[0] $xmin $xmax  $dx $opt_nbin\n";
    } 
    # try log
    $go=1;
    for ($i=2;$i<$nb;$i++){
	if($val2[$i-2] == 0.0 || $val2[$i] == 0.0 ||
	   abs(($val2[$i-1]/$val2[$i-2])/($val2[$i]/$val2[$i-1])-1.0)>0.01) {
	   $go=0;
       }
    }
    if ( $go ) {
	$opt_nbin=$nb;
	$dx=($val2[$nb-1]/$val2[0])**(1.0/($nb-1));
	$xmin=$val2[0]/sqrt($dx);
	$xmax=$val2[$nb-1]*sqrt($dx);
	print STDERR "auto-log $xmin $xmax  $dx $opt_nbin\n";
	#print "auto-log $xmin $xmax $dx $opt_nbin\n";
    } 
}

if (defined($opt_xmax)) {$xmax=$opt_xmax;}
if (defined($opt_xmin)) {$xmin=$opt_xmin;}
if (defined($opt_dx)) {$opt_nbin=($xmax-$xmin)/$opt_dx;}
for ($ix=1; $ix<=$opt_nbin; $ix++){
    $nbin[$ix]=0.;
    $bin[$ix]=0.;
    $bin2[$ix]=0.;
}

for ($i=0; $i<$n; $i++) {
    if ( $opt_log ) {
	if($val[$i]>$xmin){
	    $ix=int(log($val[$i]/$xmin)/log($xmax/$xmin) * ($opt_nbin)+1);
	} else {
	    $ix=-1;
	}
    }else{
	$ix=int(($val[$i]-$xmin)/($xmax-$xmin) * ($opt_nbin)+1);
    }
    $nbin[$ix] += 1 if $ix > 0 ;
    $bin[$ix] += $weight[$i] if $ix > 0 ;
    $bin2[$ix] += $weight[$i]**2 if $ix > 0 ;
    print "bin:xmin:$xmin val:$val[$i] xmax:$xmax w:$weight[$i] bin:$bin[$ix] ix:$ix\n" if $opt_v;
}

$s=0.0;
$w=0.0;
for ($ix=1; $ix<=$opt_nbin; $ix++){
    if( $opt_log ) {	
	$x=$xmin*($xmax/$xmin)**(($ix-0.5)/($opt_nbin));
	$dx=$xmin*($xmax/$xmin)**(($ix)/($opt_nbin))-$xmin*($xmax/$xmin)**(($ix-1.0)/($opt_nbin));
    } else {
	$x=$xmin+(1.0*$ix-0.5)*($xmax-$xmin)/($opt_nbin);
	$dx=($xmax-$xmin)/($opt_nbin);
    }
    $ddx=$dx;
    if ( $opt_norm == 0 ) {$dx=1;} else { $dx=$dx*$opt_norm;}
    $dy=0.0;
    if ( ! $opt_n ) {
	$y=$bin[$ix]/$dx;#/$n;
    }
    if( $opt_n && $nbin[$ix] != 0 ) {
	$y=$bin[$ix]/$nbin[$ix]/$dx,"   ";
	if ($opt_s) {
	    $dy=sqrt($bin2[$ix]/$nbin[$ix]-($bin[$ix]/$nbin[$ix])*($bin[$ix]/$nbin[$ix]))/$dx;
        } else {
	    $dy=sqrt($bin2[$ix]/$nbin[$ix]-($bin[$ix]/$nbin[$ix])*($bin[$ix]/$nbin[$ix]))/sqrt($nbin[$ix])/$dx;
        }
    }
    print " $x  $y  $dy  \n" unless $opt_n && $nbin[$ix] ==0 ;
    $s+=$x*$y;
    $w+=$y;
}
if ( $n != 0 ) {
    print STDERR "total:",$n,"mean:",$s/$w,"\n";
} else {
    print STDERR "total:",$n,"\n";
}
print "---------- \n" if $opt_v;
print "$xmin $xmax \n" if $opt_v;