use FindBin;
use local::lib "$FindBin::Bin/tmplib";

use warnings;
use strict;
use feature ':5.10';
my @accel = ([-0.7437,0.1118,-0.5367],
             [-0.5471,0.0062,-0.6338],
             [-0.6437,0.1216,-0.5255],
             [-0.4437,0.3216,-0.3255],
             # add three more observations, to show it works with more than 4 points (with any 3N+1)
             [-0.2718,0.3142,-0.1618],
             [-0.4567,0.1234,+0.1618],
             [-0.9876,0.5432,+0.3162],
);  # note: I changed from an array whose only element is an arrayref of arrayrefs to an array who has n arrayrefs, to be a simple 2D array

use Math::Spline;



die "requires 3N+1 points, supplied ".scalar(@accel)." points" unless 1 == scalar(@accel) % 3;
my $first = 0;
my @interp_accel = ();
my $NSTEPS = 100;
my $dt = 1/$NSTEPS;
my $grp = 0;
while($first+3 < @accel) {
    my (@idx) = map {$first+$_} 0..3;
    my $t = @interp_accel ? $dt : 0;
    while($t < 1+$dt) {
        my @m = (   (1-$t)**3,
                    (1-$t)**2 * ($t)**1 * 3,
                    (1-$t)**1 * ($t)**2 * 3,
                                ($t)**3
                );
        my $x = 0; $x += $m[$_] * $accel[$idx[$_]][0] for 0..3;
        my $y = 0; $y += $m[$_] * $accel[$idx[$_]][1] for 0..3;
        my $z = 0; $z += $m[$_] * $accel[$idx[$_]][2] for 0..3;

        printf "interpolate t=%.3f => [%+-7.4f,%+-7.4f,%+-7.4f]\n", $grp+$t, $x, $y, $z if $t<=2*$dt or $t>=$dt*($NSTEPS-2);  # debug print
        print "...\n" if $t==3*$dt;
        push @interp_accel, [$x,$y,$z];

        $t += $dt;
    }
    ++$grp;
    $first += 3;
}