#!/usr/bin/perl -w

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

use strict;
use Sci::Science;

my $n = 10;
my $pi = Sci::pi();
my $xx = 0.5;

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

my $diff = DifferentiateBy3rdSpline(\@x, \@y, $xx);
my $z = cos($xx);

print "Differentiate by 3rd Spline\n";
print "  dy/dx=$diff at x=$xx (z=$z)\n";

exit;

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

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

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

	for(my $i = 1 ; $i < $nData ; $i++) {
		if($x[$i-1] <= $x and $x <= $x[$i]) {
			my $hi  = $ph->[$i];
			my $hi2 = $hi * $hi;
			my $dxp2 = ($pX->[$i] - $x) * ($pX->[$i] - $x);
			my $dxm2 = ($x - $pX->[$i-1]) * ($x - $pX->[$i-1]);
			my $diff = ($psp->[$i-1] * ($hi2 - 3.0 * $dxp2)
					  + $psp->[$i] * (3.0 * $dxm2 - $hi2)
					  + 6.0 * ($pY->[$i] - $pY->[$i-1]) ) / (6.0 * $hi);
			return $diff;
		}
	}
	return undef;
}

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);
}