#!/usr/bin/perl
#			ignore white lines and #commented lines
#			formulas derived in Geometry.m
use strict;
use warnings;

my $SMALL=1.e-12; # threshold to decide 3 points are equally spaced
my $TINY=1.e-15; # threshold to eliminate too close points
my $NEGLIGIBLE=1.e-200; # threshold to eliminate too close points

sub print_help {
    print STDERR <<EOF;
usage:  peaks filenames [-b] [-t threshold]
	       compute all peaks in a x-sorted 2-column datafile x f(x)
	       by default use a 3-point Lorentzian/parabolic fit
	       except with option -b: yield the brutal max points in the file
	       -t threshold sets a minimum fractional peaky size for detection
	       -at athres sets a minimum absolute peak size for detection
	           by default both thresholds are = 0 
output: 3 columns with:         x_peak  extrapolated_f_peak(x) error_on_x_peak
                                              v. 2.7 by Nick Manini, 31-05-2011
						(bug fix provided by C. Negri)
EOF
return;
}

my $brutal=0;
my @inputfiles;
my $thresh=0;
my $athresh=0;
my @result;

while(my $p=shift(@ARGV)){
    if($p=~/^[^\-]/){#filename unless starting with -
        if(-e $p){
            push(@inputfiles,$p);
        }else{
            print STDERR "File $p does not exist!!!\n";
        }
        next;
    }
    if($p eq "-h"){# [-h] print help
        print_help;
        exit;
    }
    if($p eq "-b"){# forces brutal max
        $brutal=1;
        next;
    }
    if($p eq "-t"){# threshold
        $thresh=shift (@ARGV)+0.;
	($thresh>=0 && $thresh<1) || die "error! peaks -t requires 0 <= threshold < 1\n";
        next;
    }
    if($p eq "-at"){# absolute threshold
        $athresh=shift (@ARGV)+0.;
	$athresh>=0 || die "error! peaks -t requires 0 <= threshold < 1\n";
        next;
    }
    print_help;
    print STDERR "invalid option $p!\n ";
    exit;
}


if($#inputfiles<0){
    push(@inputfiles,"-");
}


my @xx;
my @yy;

my $xmax=0;
my $ymax=0;
my $width=0;
my $type="";
my $sus;
my $miny;

while (my $filename = shift(@inputfiles)){
    open(INF, $filename) || die "Cannot open file $filename\n";
    while (<INF>){
	if($_ !~/^\#/ && $_ !~/^\s*$/){ # ignore commented and empty lines
	    s/^\s+//;
	    my @fi =  split(/\s+/);
	    my $xnew=shift(@fi)+0;
	    my $ynew=shift(@fi)+0;
#	    print "debug: ",$#xx," ",$xx[$#xx]," ",$xnew,"\n";
	    if($#xx==-1 || $xnew-$xx[$#xx]>$TINY){# ignore near/unsorted points
		push @xx, $xnew;
		push @yy, $ynew;

		if($#xx<2){	# initialization
		    if($#xx==0){$miny=$ynew;}
		}else{# enough data: proceed
		    if($#result<0 && $ynew<$miny){
			$miny=$ynew; # the lowest y before here or peak
		    }
		    $sus=abs($yy[0])+abs($yy[1])+abs($yy[2]);
		    if($yy[2] -$yy[1] <0 && $yy[1] - $yy[0] >= 0 # local max 
		       && $sus>$NEGLIGIBLE){
#			print "#debug ",$xx[0]," ",$yy[0]," ",$xx[1]," ",$yy[1]," ",$xx[2]," ",$yy[2]," ";
#	  determination of the max coordinates $xmax $ymax and x-width $width
			if($brutal){
			    $xmax=$xx[1];
			    $ymax=$yy[1];
			    $width=0.5*($xx[2]-$xx[0]);
			}else{
#			    equally spaced case, very robust formula:
			    if((($xx[1]-$xx[0]- $xx[2]+$xx[1])/($xx[1]-$xx[0]))**2 < $SMALL){
				my $d=0.5*($xx[2]-$xx[0]);
				my $dena=$yy[1]*($yy[0]+$yy[2])-2*$yy[0]*$yy[2];
				my $numb=$d*$d * (
				    8.*$yy[1]*$yy[0]*$yy[2]*($yy[0] + $yy[2])
				    -$yy[1]**2*($yy[0] - $yy[2])**2
				    - (4.*$yy[0]*$yy[2])**2
				    );
#				print "debug ",$dena," ",$numb,"\n";
				if($numb*$dena*$dena>0){ # Lorentzian fit ok
				    $type="Lorentzian";
				    $xmax=$xx[1]+0.5*$d *$yy[1]*($yy[2]-$yy[0])/$dena;
				    my $a=2.*$d*$d*$yy[0]*$yy[1]*$yy[2]/$dena;
				    my $b=$numb/(4.*$dena*$dena);
				    $ymax=$a/$b;
				    $width=sqrt($b);
				}else{#			must use parabolic fit:
				    $type="parabola";
				    my $den=(2*$yy[1] -$yy[0] -$yy[2]);
				    $xmax=$xx[1] + $d*0.5*($yy[2]-$yy[0])/$den;
				    $ymax=
					($yy[0]**2 + (-4.*$yy[1] + $yy[2])**2
					 -2.*$yy[0]*(4.*$yy[1] + $yy[2]))/
					 (8.*$den);
				    $width=$d/(2*$den) *sqrt(abs($yy[0]**2 - 8*$yy[0]*$yy[1] + 16*$yy[1]**2 - 2*$yy[0]*$yy[2] - 8*$yy[1]*$yy[2] + $yy[2]**2));
				}
			    }else{ # generic spacing: much noisier formulae 
				my $den=4.*($xx[0]*$yy[0]*($yy[1] - $yy[2]) 
					    +$xx[2]*($yy[0] - $yy[1])*$yy[2] 
					    +$xx[1]*$yy[1]*(-$yy[0] + $yy[2]))**2;
				if($den>0){
				    $b=-(($xx[0]-$xx[1])**4*$yy[0]**2*$yy[1]**2
					 -2*($xx[0] - $xx[1])**2*$yy[0]*$yy[1]*
					 (($xx[0] - $xx[2])**2*$yy[0] 
					  +($xx[1] - $xx[2])**2*$yy[1])*$yy[2]
					 + (($xx[0] - $xx[2])**2*$yy[0] 
					    -($xx[1] - $xx[2])**2*$yy[1])**2
					 *$yy[2]**2)/$den;
				    if($b>0){#		Lorentzian fit ok:
					$type="Lorentzian";
					$xmax=($xx[0]**2*$yy[0]*($yy[1] - $yy[2]) + $xx[2]**2*($yy[0] - $yy[1])*$yy[2] + $xx[1]**2*$yy[1]*(-$yy[0] + $yy[2])) / (2*($xx[0]*$yy[0]*($yy[1] - $yy[2]) + $xx[2]*($yy[0] - $yy[1])*$yy[2] + $xx[1]*$yy[1]*(-$yy[0] + $yy[2])));
					$ymax= (-4*($xx[0] - $xx[1])*($xx[0] - $xx[2])*($xx[1] - $xx[2])*$yy[0]*$yy[1]*$yy[2]*($xx[0]*$yy[0]*($yy[1] - $yy[2]) + $xx[2]*($yy[0] - $yy[1])*$yy[2] + $xx[1]*$yy[1]*(-$yy[0] + $yy[2])))/(($xx[0] - $xx[1])**4*$yy[0]**2*$yy[1]**2 - 2*($xx[0] - $xx[1])**2*$yy[0]*$yy[1]*(($xx[0] - $xx[2])**2*$yy[0] + ($xx[1] - $xx[2])**2*$yy[1])*$yy[2] + (($xx[0] - $xx[2])**2*$yy[0] - ($xx[1] - $xx[2])**2*$yy[1])**2*$yy[2]**2);
					$width=sqrt($b);
				    }else{#	must use parabolic fit:
					$type="parabola";
					$xmax=($xx[2]**2*($yy[0] - $yy[1]) + $xx[0]**2*($yy[1] - $yy[2]) + $xx[1]**2*($yy[2] - $yy[0]))/(2*($xx[2]*($yy[0] - $yy[1]) + $xx[0]*($yy[1] - $yy[2]) + $xx[1]*($yy[2] - $yy[0])));
					$ymax=((($xx[1] - $xx[2])**2*$yy[0] - ($xx[0] - $xx[2])**2*$yy[1])**2 -2*($xx[0] - $xx[1])**2*(($xx[1] - $xx[2])**2*$yy[0] + ($xx[0] - $xx[2])**2*$yy[1])*$yy[2] +($xx[0] - $xx[1])**4*$yy[2]**2)/
					    (4.*($xx[0] - $xx[1])*($xx[0] - $xx[2])*($xx[1] - $xx[2])*($xx[2]*($yy[0] - $yy[1]) + $xx[0]*($yy[1] - $yy[2]) + $xx[1]*(-$yy[0] + $yy[2])));
					$den=2.*
					    ($xx[2]*($yy[0] -$yy[1]) + $xx[1]*($yy[2]-$yy[0]) +$xx[0]*($yy[1] -$yy[2]))/(($xx[0] -$xx[1])*($xx[0] -$xx[2])*($xx[1] -$xx[2]));
					if($den>0){
					    $width=sqrt(abs(((($xx[1]-$xx[2])**2*$yy[0]-($xx[0]-$xx[2])**2*$yy[1])**2 -2*($xx[0] - $xx[1])**2*(($xx[1] - $xx[2])**2*$yy[0] +($xx[0] - $xx[2])**2*$yy[1])*$yy[2] + ($xx[0] - $xx[1])**4*$yy[2]**2)/(($xx[0] -$xx[1])**2*($xx[0] -$xx[2])**2*($xx[1] - $xx[2])**2)))/$den;
					}else{
					    $width=0;
					}
				    }	# parabolic else
				}	# if($den>0): here a else is needed
			    }	# generic much noisier formulae 
			}	# if($brutal) ... else{
#		end determination of the max coordinates
			if($#result<0 || $ymax>$result[1]){	# insert/update
			    if(($ymax - $miny) >= $thresh*$sus
			       && ($ymax - $miny) >= $athresh){
				@result=($xmax,$ymax,$width,$type,$sus);
			    }
			}
		    } # max found
#		    print "#check ",$#result," ",$result[1]," ",$ynew," ",$miny,"\n";
		    if($#result>0
		       && ($result[1] - $ynew) >= $thresh*$result[4]
		       && ($result[1] - $ynew) >= $athresh){
			print $result[0],"\t",$result[1],"\t",$result[2],"\t",$result[3],"\n";
			@result=();
			$miny=$ynew;
		    }
		    shift @xx; shift @yy; # drop old info
		} # ($#xx>=2)
	    } # ($xnew-$xx[$#xx-1]>$TINY)
	} # ($_ !~/^\#/ && $_ !~/^\s*$/)
    } # while: $filename = shift(@ARGV))
    close(INF);
}