interp

package
v3.0.0-...-8217f41 Latest Latest
Warning

This package is not in the latest version of its module.

 Go to latest Published: Jun 1, 2018 License: MIT Imports: 3 Imported by: 0

Documentation

Overview

Interp: Chapter 3, Interpolation.

Len3 and Len5 types

These types allow interpolation from a table of equidistant x values and corresponding y values. Since the x values are equidistant, only the first and last values are supplied as arguments to the constructors. The interior x values are implicit. All y values must be supplied however. They are passed as a slice, and the length of y is fixed. For Len3 it must be 3 and for Len5 it must be 5.

For these Len3 and Len5 functions, Meeus notes the importance of choosing the 3 or 5 rows of a larger table that will minimize the interpolating factor n. He does not provide algorithms for doing this however.

For an example of a selection function, see Len3ForInterpolateX. This was useful for computing Delta T.

Index

Examples

Constants

This section is empty.

Variables

View Source 
var (
	ErrorNot3            = errors.New("Argument y must be length 3")
	ErrorNot4            = errors.New("Argument y must be length 4")
	ErrorNot5            = errors.New("Argument y must be length 5")
	ErrorNoXRange        = errors.New("Argument x3 (or x5) cannot equal x1")
	ErrorNOutOfRange     = errors.New("Interpolating factor n must be in range -1 to 1")
	ErrorXOutOfRange     = errors.New("Argument x outside of range x1 to x3 (or x5)")
	ErrorNoExtremum      = errors.New("No extremum in table")
	ErrorExtremumOutside = errors.New("Extremum falls outside of table")
	ErrorZeroOutside     = errors.New("Zero falls outside of table")
	ErrorNoConverge      = errors.New("Failure to converge")
)

Error values returned by functions and methods in this package. Defined here to help testing for specific errors.

Functions

func Lagrange 

func Lagrange(x float64, table []struct{ X, Y float64 }) (y float64)

Lagrange performs interpolation with unequally-spaced abscissae.

Given a table of X and Y values, interpolate a new y value for argument x.

X values in the table do not have to be equally spaced; they do not even have to be in order. They must however, be distinct. table 中包含了 n 个点且xi 必须互异,x 为目标插值点

Example

exercise, p. 34. 

package main

import (
	"fmt"

	"github.com/mooncaker816/learnmeeus/v3/interp"
)

func main() {
	table := []struct{ X, Y float64 }{
		{29.43, .4913598528},
		{30.97, .5145891926},
		{27.69, .4646875083},
		{28.11, .4711658342},
		{31.58, .5236885653},
		{33.05, .5453707057},
	}
	// 10 significant digits in input, no more than 10 expected in output
	fmt.Printf("30: %.10f\n", interp.Lagrange(30, table))
	fmt.Printf("0:  %.10f\n", interp.Lagrange(0, table))
	fmt.Printf("90: %.10f\n", interp.Lagrange(90, table))
}

Output:

30: 0.5000000000
0:  0.0000512249
90: 0.9999648100

func LagrangePoly 

func LagrangePoly(table []struct{ X, Y float64 }) []float64

LagrangePoly uses the formula of Lagrange to produce an interpolating polynomial.

X values in the table do not have to be equally spaced; they do not even have to be in order. They must however, be distinct.

The returned polynomial will be of degree n-1 where n is the number of rows in the table. It can be evaluated for x using common.Horner.

Example 
package main

import (
	"fmt"

	"github.com/mooncaker816/learnmeeus/v3/interp"
)

func main() {
	// Example 3.g, p, 34.
	table := []struct{ X, Y float64 }{
		{1, -6},
		{3, 6},
		{4, 9},
		{6, 15},
	}
	p := interp.LagrangePoly(table)
	// output format contrived to fit expected result
	for _, c := range p {
		fmt.Printf("%.0f\n", c*5)
	}
}

Output:

-87
69
-13
1

func Len4Half 

func Len4Half(y []float64) (float64, error)

Len4Half interpolates a center value from a table of four rows.

Example 
// Example 3.f, p. 32.
half, err := interp.Len4Half([]float64{
	unit.FromSexa(0, 10, 18, 48.732),
	unit.FromSexa(0, 10, 23, 22.835),
	unit.FromSexa(0, 10, 27, 57.247),
	unit.FromSexa(0, 10, 32, 31.983),
})
if err != nil {
	fmt.Println(err)
	return
}
fmt.Printf("%.3d", sexa.FmtRA(unit.RAFromHour(half)))

Output:

10ʰ25ᵐ40ˢ.001

Types

type Len3 

type Len3 struct {
	// contains filtered or unexported fields
}

Len3 allows second difference interpolation. 等距三点插值结构

func Len3ForInterpolateX 

func Len3ForInterpolateX(x, x1, xn float64, y []float64) (*Len3, error)

Len3ForInterpolateX is a special purpose Len3 constructor.

Like NewLen3, it takes a table of x and y values, but it is not limited to tables of 3 rows. An X value is also passed that represents the interpolation target x value. Len3ForInterpolateX will locate the appropriate three rows of the table for interpolating for x, and initialize the Len3 object for those rows. 

x is the target for interpolation
x1 is the x value corresponding to the first y value of the table.
xn is the x value corresponding to the last y value of the table.
y is all y values in the table.  len(y) should be >= 3.

给定n个点，但是我们只需选取离目标点x最接近的三个点来做三点插值， 此时可用以下函数来自动选择最优三点，来构造三点插值 同样，前提是n个点等距，且与y一一对应

func NewLen3 

func NewLen3(x1, x3 float64, y []float64) (*Len3, error)

NewLen3 prepares a Len3 object from a table of three rows of x and y values.

X values must be equally spaced, so only the first and last are supplied. X1 must not equal x3. Y must be a slice of three y values. 根据上述定义，创建三点插值结构

func (*Len3) Extremum 

func (d *Len3) Extremum() (x, y float64, err error)

Extremum returns the x and y values at the extremum.

Results are restricted to the range of the table given to the constructor NewLen3.

Example 
// Example 3.b, p. 26.
d3, err := interp.NewLen3(12, 20, []float64{
	1.3814294,
	1.3812213,
	1.3812453,
})
if err != nil {
	fmt.Println(err)
	return
}
x, y, err := d3.Extremum()
if err != nil {
	fmt.Println(err)
	return
}
fmt.Printf("distance:  %.7f AU\n", y)
fmt.Printf("date:     %.4f\n", x)
i, frac := math.Modf(x)
fmt.Printf("1992 May %d, at %h TD",
	int(i), sexa.FmtTime(unit.TimeFromDay(frac)))

Output:

distance:  1.3812030 AU
date:     17.5864
1992 May 17, at 14ʰ TD

func (*Len3) InterpolateN 

func (d *Len3) InterpolateN(n float64) (y float64)

InterpolateN interpolates for a given interpolating factor n.

This is interpolation formula (3.3)

The interpolation factor n is x-x2 in units of the tabular x interval. (See Meeus p. 24.) 非严格插值计算，不用保证目标点插值因子绝对值小于等于1， 即不用保证离我们所选三点中点距离小于一个步长

Example 
package main

import (
	"fmt"

	"github.com/mooncaker816/learnmeeus/v3/interp"
	"github.com/soniakeys/unit"
)

func main() {
	// Example 3.a, p. 25.
	d3, err := interp.NewLen3(7, 9, []float64{
		.884226,
		.877366,
		.870531,
	})
	if err != nil {
		fmt.Println(err)
		return
	}
	h := unit.FromSexa(0, 4, 21, 0)
	fmt.Println(h, "hours")
	n := h / 24
	y := d3.InterpolateN(n)
	fmt.Printf("%.6f\n", y)
}

Output:

4.35 hours
0.876125

func (*Len3) InterpolateNStrict 

func (d *Len3) InterpolateNStrict(n float64) (y float64, err error)

InterpolateNStrict interpolates for a given interpolating factor n.

N is restricted to the range [-1..1] corresponding to the range x1 to x3 given to the constructor NewLen3. 严格插值计算，必须保证目标点插值因子绝对值小于等于1， 即必须保证离我们所选三点中点距离小于一个步长

func (*Len3) InterpolateX 

func (d *Len3) InterpolateX(x float64) (y float64)

InterpolateX interpolates for a given x value. 计算插值因子n，调用非严格插值计算

Example 
package main

import (
	"fmt"

	"github.com/mooncaker816/learnmeeus/v3/interp"
	"github.com/soniakeys/unit"
)

func main() {
	// Example 3.a, p. 25.
	d3, err := interp.NewLen3(7, 9, []float64{
		.884226,
		.877366,
		.870531,
	})
	if err != nil {
		fmt.Println(err)
		return
	}
	x := 8 + unit.NewTime(' ', 4, 21, 0).Day() // 8th day at 4:21
	y := d3.InterpolateX(x)
	fmt.Printf("%.6f\n", y)
}

Output:

0.876125

func (*Len3) InterpolateXStrict 

func (d *Len3) InterpolateXStrict(x float64) (y float64, err error)

InterpolateXStrict interpolates for a given x value, restricting x to the range x1 to x3 given to the constructor NewLen3. 计算插值因子n，调用严格插值计算

func (*Len3) Zero 

func (d *Len3) Zero(strong bool) (x float64, err error)

Len3Zero finds a zero of the quadratic function represented by the table.

That is, it returns an x value that yields y=0.

Argument strong switches between two strategies for the estimation step. when iterating to converge on the zero.

Strong=false specifies a quick and dirty estimate that works well for gentle curves, but can work poorly or fail on more dramatic curves.

Strong=true specifies a more sophisticated and thus somewhat more expensive estimate. However, if the curve has quick changes, This estimate will converge more reliably and in fewer steps, making it a better choice.

Results are restricted to the range of the table given to the constructor NewLen3. strong 为考虑修正量的迭代方式，更为精确

Example 
// Example 3.c, p. 26.
x1 := 26.
x3 := 28.
// the y unit doesn't matter.  working in degrees is fine
yTable := []float64{
	unit.FromSexa('-', 0, 28, 13.4),
	unit.FromSexa(' ', 0, 6, 46.3),
	unit.FromSexa(' ', 0, 38, 23.2),
}
d3, err := interp.NewLen3(x1, x3, yTable)
if err != nil {
	fmt.Println(err)
	return
}
x, err := d3.Zero(false)
if err != nil {
	fmt.Println(err)
	return
}
fmt.Printf("February %.5f\n", x)
i, frac := math.Modf(x)
fmt.Printf("February %d, at %m TD",
	int(i), sexa.FmtTime(unit.TimeFromDay(frac)))

Output:

February 26.79873
February 26, at 19ʰ10ᵐ TD
Example (Strong) 
package main

import (
	"fmt"

	"github.com/mooncaker816/learnmeeus/v3/interp"
)

func main() {
	// Example 3.d, p. 27.
	x1 := -1.
	x3 := 1.
	yTable := []float64{-2, 3, 2}
	d3, err := interp.NewLen3(x1, x3, yTable)
	if err != nil {
		fmt.Println(err)
		return
	}
	x, err := d3.Zero(true)
	if err != nil {
		fmt.Println(err)
		return
	}
	fmt.Printf("%.12f\n", x)
}

Output:

-0.720759220056

type Len5 

type Len5 struct {
	// contains filtered or unexported fields
}

Len5 allows fourth difference interpolation. 五点等距插值结构

func NewLen5 

func NewLen5(x1, x5 float64, y []float64) (*Len5, error)

NewLen5 prepares a Len5 object from a table of five rows of x and y values.

X values must be equally spaced, so only the first and last are supplied. X1 must not equal x5. Y must be a slice of five y values. 构造5点插值结构

func (*Len5) Extremum 

func (d *Len5) Extremum() (x, y float64, err error)

Extremum returns the x and y values at the extremum.

Results are restricted to the range of the table given to the constructor NewLen5. (Meeus actually recommends restricting the range to one unit of the tabular interval, but that seems a little harsh.) 5点插值极值

func (*Len5) InterpolateN 

func (d *Len5) InterpolateN(n float64) (y float64)

InterpolateN interpolates for a given interpolating factor n.

The interpolation factor n is x-x3 in units of the tabular x interval. (See Meeus p. 28.) Horner 为工具函数，求解多项式之和，interpCoeff为多项式系数

func (*Len5) InterpolateNStrict 

func (d *Len5) InterpolateNStrict(n float64) (y float64, err error)

InterpolateNStrict interpolates for a given interpolating factor n.

N is restricted to the range [-1..1]. This is only half the range given to the constructor NewLen5, but is the recommendation given on p. 31.

func (*Len5) InterpolateX 

func (d *Len5) InterpolateX(x float64) (y float64)

InterpolateX interpolates for a given x value.

Example 
// Example 3.e, p. 28.
x1 := 27.
x5 := 29.
// work in degrees
yTable := []float64{
	unit.FromSexa(' ', 0, 54, 36.125),
	unit.FromSexa(' ', 0, 54, 24.606),
	unit.FromSexa(' ', 0, 54, 15.486),
	unit.FromSexa(' ', 0, 54, 08.694),
	unit.FromSexa(' ', 0, 54, 04.133),
}
d5, err := interp.NewLen5(x1, x5, yTable)
if err != nil {
	fmt.Println(err)
	return
}
n := (3 + 20./60) / 24
x := 28 + n
y := d5.InterpolateX(x)
fmt.Printf("%.3d", sexa.FmtAngle(unit.AngleFromDeg(y)))

Output:

54′13″.369

func (*Len5) InterpolateXStrict 

func (d *Len5) InterpolateXStrict(x float64) (y float64, err error)

InterpolateXStrict interpolates for a given x value, restricting x to the range x1 to x5 given to the the constructor NewLen5.

func (*Len5) Zero 

func (d *Len5) Zero(strong bool) (x float64, err error)

Len5Zero finds a zero of the quartic function represented by the table.

That is, it returns an x value that yields y=0.

Argument strong switches between two strategies for the estimation step. when iterating to converge on the zero.

Strong=false specifies a quick and dirty estimate that works well for gentle curves, but can work poorly or fail on more dramatic curves.

Strong=true specifies a more sophisticated and thus somewhat more expensive estimate. However, if the curve has quick changes, This estimate will converge more reliably and in fewer steps, making it a better choice.

Results are restricted to the range of the table given to the constructor NewLen5. strong 为带修正模式

Example 
// Exercise, p. 30.
x1 := 25.
x5 := 29.
yTable := []float64{
	unit.FromSexa('-', 1, 11, 21.23),
	unit.FromSexa('-', 0, 28, 12.31),
	unit.FromSexa(' ', 0, 16, 07.02),
	unit.FromSexa(' ', 1, 01, 00.13),
	unit.FromSexa(' ', 1, 45, 46.33),
}
d5, err := interp.NewLen5(x1, x5, yTable)
if err != nil {
	fmt.Println(err)
	return
}
z, err := d5.Zero(false)
if err != nil {
	fmt.Println(err)
	return
}
fmt.Printf("1988 January %.6f\n", z)
zInt, zFrac := math.Modf(z)
fmt.Printf("1988 January %d at %m TD\n", int(zInt),
	sexa.FmtTime(unit.TimeFromDay(zFrac)))

// compare result to that from just three central values
d3, err := interp.NewLen3(26, 28, yTable[1:4])
if err != nil {
	fmt.Println(err)
	return
}
z3, err := d3.Zero(false)
if err != nil {
	fmt.Println(err)
	return
}
dz := z - z3
fmt.Printf("%.6f day\n", dz)
fmt.Printf("%.1f minute\n", dz*24*60)

Output:

1988 January 26.638587
1988 January 26 at 15ʰ20ᵐ TD
0.000753 day
1.1 minute

Source Files

View all Source files

Jump to

Keyboard shortcuts

? : This menu
/ : Search site
f or F : Jump to
y or Y : Canonical URL
go.dev uses cookies from Google to deliver and enhance the quality of its services and to analyze traffic. Learn more.