vector/interp.go

49 lines
1.4 KiB
Go
Raw Normal View History

2014-04-23 14:29:38 -07:00
package vector
import (
"errors"
"fmt"
)
// explicitly implements Neville's algorithm for three points at very
2014-01-12 23:05:25 -08:00
// particular parameter values: ts[]{0, 1/3, 2/3}, extrapolate to 1.
func Extrapolate(P []Point2d) (*Point2d, error) {
if len(P) != 3 {
return nil, errors.New("only works for len(P) == 3")
}
// P_{01} = -2 * P_0 + 3 * P_1
p01 := P[0].ToVector().Scale(-2.0).ToPoint().Add(P[1].ToVector().Scale(3.0))
// P_{12} = -1 * P_1 + 2 * P_2
p12 := P[1].ToVector().Scale(-1.0).ToPoint().Add(P[2].ToVector().Scale(2.0))
// P_{012} = 1/2 * P_{01} + 3/2 * P_{12}
p012 := p01.ToVector().Scale(-1.0 / 2.0).ToPoint().Add(p12.ToVector().Scale(3.0 / 2.0))
return &p012, nil
}
2014-01-12 23:05:25 -08:00
// Implements Neville's Algorithm for a slice of points (P) at parameter
// values(ts), to parameter value (t). Returns the point and an error
func Neville(P []Point2d, ts []float64, t float64) (*Point2d, error) {
if len(P) != len(ts) {
return nil, errors.New(
fmt.Sprintf(
"Incompatable slice lengths len(P): %d != len(ts): %d",
len(P),
len(ts),
),
)
}
2014-01-12 23:05:25 -08:00
Q := []Point2d(P)
order := len(P) - 1
for i := 0; i < order; i++ {
for j := 0; j < order-i; j++ {
var tLeft float32 = float32(ts[i+j+1] - t)
var tRight float32 = float32(t - ts[j])
Q[j] = Q[j].ToVector().Scale(tLeft).ToPoint().Add(P[j+1].ToVector().Scale(tRight))
Q[j] = Q[j].ToVector().Scale(float32(1.0 / (ts[i+j+1] - ts[j]))).ToPoint()
}
}
return &Q[0], nil
}