#!/usr/bin/perl -w

# Aufruf: ./monotone_cubic_interpolation.pl output_file

##############################################################################
# USER INPUT
##############################################################################

# output format
format OUT =
@>>>>>>>>>>>>>>>>>>>@>>>>>>>>>>>>>>>>>>>
$x, $y
.

# number of points for output curve
$points = 100;

# if output curve has to be monotone for monotone input data, set monotone = 1
$monotone = 1;

##############################################################################
open(IN, "UserPreproOpt.inp");
open(OUT,">$ARGV[0]");
open(OUT2,">$ARGV[0].history");

print OUT "*DEFINE_CURVE\n";
print OUT "\$Dehnrate für 1 1/s\n";
print OUT "\$     LCID      SIDR      SCLA      SCLO      OFFA      OFFO    DATTYP\n";
print OUT "         1         0       1.0  1.000000\n";
print OUT "\$                 a1                  o1\n";

# read input file
$k=0;
while (<IN>){
    $values{$k} = [ split ];
    $k++;
}


# number of sampling points
$n = $k;

$range = $values{$n-1}[0] - $values{0}[0]; 

# initialize tangents at every data point
$m{0} = ($values{1}[1]-$values{0}[1]) / ($values{1}[0]-$values{0}[0]);

for ($k=1;$k<$n-1;$k++){

    $m{$k} =  ($values{$k}[1]-$values{$k-1}[1]) / (2*($values{$k}[0]-$values{$k-1}[0])) + ($values{$k+1}[1]-$values{$k}[1]) / (2*($values{$k+1}[0]-$values{$k}[0]));

}

$m{$n-1} = ($values{$n-1}[1]-$values{$n-2}[1]) / ($values{$n-1}[0]-$values{$n-2}[0]);

if ($monotone == 1){
# monotony conditions
    for ($k=0;$k<$n-1;$k++){


	$d = ($values{$k+1}[1] - $values{$k}[1])/($values{$k+1}[0] - $values{$k}[0]);

	if ($d == 0){
	    $m{$k}   = 0;
	    $m{$k+1} = 0;
	}else{
	    $a = $m{$k}/$d;
	    $b = $m{$k+1}/$d;

	    if ($a+$b-2>0 && 2*$a+$b-3>0 && $a+2*$b-3>0 && $a-1/3*(2*$a+$b-3)*(2*$a+$b-3)/($a+$b-2)<0){
		
		$m{$k}   = 3/sqrt($a*$a+$b*$b)*$a*$d;
		$m{$k+1} = 3/sqrt($a*$a+$b*$b)*$b*$d;

	    }
	}

    }
}

# Cubic Hermite Splines
for ($k=0;$k<$n-1;$k++){
    $x0 = $values{$k}[0];
    $x1 = $values{$k+1}[0];
    $y0 = $values{$k}[1];
    $y1 = $values{$k+1}[1];

    $h = $x1 - $x0;
     
    $points_k = $points*$h/$range;

    for ($i=0; $i<$points_k; $i++ ){
	
      $x = $x0 + $i*$h/$points_k;
       
      if ($x != $x1){
	$t = ($x -$x0)/$h; 
        
        # Hermit basis functions
	$h00 = 2*$t*$t*$t -3*$t*$t +1;
	$h10 = $t*$t*$t -2*$t*$t + $t;
	$h01 = -2*$t*$t*$t + 3*$t*$t;
	$h11 = $t*$t*$t -$t*$t;

	# Interpolation
	$y = $h00*$y0 + $h10*$h*$m{$k} + $h01*$y1 + $h11*$h*$m{$k+1}; 

	write OUT;
	print OUT2 "$x $y\n";
      }
    }

}


close(IN);
close(OUT);
close(OUT2);