#!/usr/bin/perl
#			ignore white lines and #commented lines
use strict;
use warnings;

my $TINY=1.e-15; # threshold to eliminate too close points

sub print_help {
    print STDERR <<EOF;
usage:  add_kink_points filenames
	       compute and insert all downward kink points from
	       an x-sorted 2-column datafile x f(x)
	       by crossing two linear extrapolations from before and after
output: 2 columns with:                 x f(x) (including extra points)
                                              v. 1.1 by Nick Manini, 23-07-2017
EOF
return;
}

my @inputfiles;

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;
    }
    print_help;
    print STDERR "invalid option $p!\n ";
    exit;
}

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

my @xx;
my @yy;

my $xkink=0;
my $ykink=0;

my($ma,$qa,$mb,$qb);
my $after;

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;
#	    print "debug: ",$#xx," ",$xx[$#xx]," ",$xnew,"\n";
	    if($#xx==-1 || $xnew-$xx[$#xx]>$TINY){# ignore near/unsorted points
		push @xx, $xnew;
		push @yy, shift(@fi)+0; 
		if($#xx>=4){ # enough points: can proceed
		    if($yy[4]>=$yy[3] && $yy[3]>$yy[2] &&
		       $yy[2]<=$yy[1] && $yy[1]<=$yy[0] ){ # min condition at 2
#			print "debug ",$xx[1]," ",$yy[1]," ",$xx[2]," ",$yy[2]," ",$xx[3]," ",$yy[3],"\n";
			($ma,$qa)=straight_line_thru_2_points($xx[0],$yy[0],$xx[1],$yy[1]);
			($mb,$qb)=straight_line_thru_2_points($xx[3],$yy[3],$xx[4],$yy[4]);

			my $ytrya=$ma*$xx[2]+$qa;
			my $ytryb=$mb*$xx[2]+$qb;
			my $deltaa=abs($ytrya-$yy[2]);
			my $deltab=abs($ytryb-$yy[2]);

			print "# debug ",$deltaa," ",$deltab,"\n";
			if($deltaa<0.1*$deltab || $deltab<0.1*$deltaa){ # must be sharp!
			    if($deltaa<$deltab){
				$after=1;
				($ma,$qa)=straight_line_thru_2_points($xx[2],$yy[2],$xx[1],$yy[1]);
			    }else{
				$after=0;
				($mb,$qb)=straight_line_thru_2_points($xx[2],$yy[2],$xx[3],$yy[3]);
			    }
			    ($xkink,$ykink)=intersect_of_2_straight_lines($ma,$qa,$mb,$qb);

			    if($after){
				print $xx[2],"\t",$yy[2],"\n";
				print $xkink,"\t",$ykink," added\n";
			    }else{
				print $xkink,"\t",$ykink," added\n";
				print $xx[2],"\t",$yy[2],"\n";
			    }
			}else{ # a smooth minimum: just print out
			    print $xx[2],"\t",$yy[2],"\n";
			}
		    }else{ # non min, just print out:
			print $xx[2],"\t",$yy[2],"\n";
		    }
		    shift @xx; shift @yy; # drop old info
		}else{ # still filling the buffer:
		    my $ind=$#xx-2;
		    if($ind>=0){
			print $xx[$ind],"\t",$yy[$ind],"\n";
		    }
		}
	    } # ($xnew-$xx[$#xx-1]>$TINY)
	} # ($_ !~/^\#/ && $_ !~/^\s*$/)
    } # while: $filename = shift(@ARGV))
    close(INF);
}
print $xx[2],"\t",$yy[2],"\n";
print $xx[3],"\t",$yy[3],"\n";


sub straight_line_thru_2_points{ # returns the angular coeff m and intercept q
# input: x1,y1,x2,y2
    my $HUGE=1.e308;
    my ($m,$q);
    if($_[0] == $_[2]){
	$m=$HUGE;
	$q=-$HUGE;
    }else{
	$m=($_[3] - $_[1])/($_[2] - $_[0]);
	$q=0.5*(($_[3] + $_[1]) - $m*($_[2] + $_[0]));
    }
    my @out=($m,$q);
    return @out;
}

sub intersect_of_2_straight_lines{ # returns x,y of the common point
# input: m1,q1,m2,q2
    my $HUGE=1.e308;
    my ($x,$y);
    if($_[0] == $_[2]){
	if($_[1] == $_[3]){
	    die "intersect_of_2_straight_lines: indeterminate problem!\n";
	}else{
	    die "intersect_of_2_straight_lines: impossible problem!\n";
	}
    }else{
	$x=($_[1] - $_[3])/($_[2] - $_[0]);
	$y=0.5*(($_[0] + $_[2])*$x + ($_[1] + $_[3]));
    }
    my @out=($x,$y);
    return @out;
}