#!/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;
	if(($i != 0) and ($i != $n)) {
		$x[$i] += 0.005 * $i;
	}
	$y[$i] = sin($x[$i]);
#print "$i: $x[$i], $y[$i]\n";
}

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

print "Integral sin(x) from x=0 to pi\n";
print "  n=", $n+1, "\n";
print "  S=$S\n";

exit;

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

	my $nData = @$pX;
	if($nData < 3 or $nData % 2 == 0) {
		return undef;
	}

	my $w = 0.0;
	for(my $i = 1 ; $i < $nData ; $i += 2) {
		my $h   = 0.5 * ($x[$i+1] - $x[$i-1]);
		my $hpd = 1.0 / ($x[$i] - $x[$i-1]);
		my $hmd = 1.0 / ($x[$i+1] - $x[$i]);
		my $d   = $x[$i] - $x[$i-1] - $h;
		$w += ( (1.0 + 2.0*$d*$hpd) * $y[$i-1]
			  +  2.0*$h*($hpd+$hmd) * $y[$i]
			  + (1.0 - 2.0*$d*$hmd) * $y[$i+1] ) * $h;
	}
	my $S = $w / 3.0;
	return $S;
}