#!/usr/bin/perl -w

use lib "d:/Programs/Perl/lib";

use strict;
use Sci::Science;

my $n = 20;
my $pi = Sci::pi();

my @x;
my @y;
for(my $i = 0 ; $i <= $n ; $i++) {
	$x[$i] = $pi * $i / $n;
	$y[$i] = sin($x[$i]);
}

my $S = IntegrateBy3rdSpline(\@x, \@y);

print "Differentiate by 3rd Spline\n";
print "  S=$S\n";

exit;

sub IntegrateBy3rdSpline
{
	my ($pX, $pY) = @_;

	my $nData = @$pX;
	if($nData < 4) {
		return undef;
	}
	for(my $i = 0 ; $i < $nData-1 ; $i++) {
		if($pX->[$i+1] <= $pX->[$i]) {
			return undef;
		}
	}

	my ($ph, $psp) = subspl($pX, $pY);

	my $w = 0.0;
	for(my $i = 1 ; $i < $nData ; $i++) {
		my $hi = $ph->[$i];
		my $hi2 = $hi * $hi;
		$w += $hi * ( $pY->[$i] + $pY->[$i-1] 
					- $hi2 * ($psp->[$i-1] + $psp->[$i]) * 8.3333333333333e-2);
	}
	return $w * 0.5;
}

sub subspl
{
	my ($pX, $pY) = @_;

	my $nData = @$pX;
	my (@h, @sp, @du, @dc);
	$du[$nData-1] = 0.0;
	$dc[$nData-1] = 0.0;
	$sp[$nData-1] = 0.0;
	for(my $i = 1 ; $i < $nData ; $i++) {
		$h[$i] = $x[$i] - $x[$i-1];
	}
	for(my $i = 1 ; $i < $nData-1 ; $i++) {
		my $hip1 = $x[$i+1] - $x[$i];
		$sp[$i] = 6.0 * (($y[$i+1] - $y[$i]) / ($h[$i] * $hip1)
			    - ($y[$i] - $y[$i-1]) / ($h[$i] * $h[$i]));
		$du[$i] = $h[$i+1] / $h[$i];
		$dc[$i] = 2.0 * (1.0 + $du[$i]);
	}
	$dc[$nData-1] += 1.0;
	$dc[$nData-2] += $h[$nData-1] / $h[$nData-2];
	$du[1] /= $dc[1];
	$sp[1] /= $dc[1];
	for(my $i = 2 ; $i < $nData ; $i++) {
#print "$i: $dc[$i], $du[$i-1]\n";
		my $g = 1.0 / ($dc[$i] - $du[$i-1]);
		$du[$i] *= $g;
		$sp[$i] = ($sp[$i] - $sp[$i-1]) * $g;
	}
	for(my $i = $nData-2 ; $i >= 1 ; $i--) {
#print "$i: $sp[$i+1], $du[$i]\n";
		$sp[$i] -= $sp[$i+1] * $du[$i];
	}
	$sp[0] = $sp[1];
	$sp[$nData-1] = $sp[$nData-2];
	return (\@h, \@sp);
}