Thiele's interpolation formula: Difference between revisions

=={{header|FreeBASIC}}==

{{trans|Phix}}

<syntaxhighlight lang="vbnet">Const As Integer n1 = 32

Const As Integer n2 = (n1 * (n1 - 1) / 2)

Const As Double paso = 0.05

Const As Double INF = 1e308

Const As Double NaN = -(INF/INF)

Dim As Double xVal(n1), tSin(n1), tCos(n1), tTan(n1)

For i As Integer = 1 To n1

xVal(i) = (i-1) * paso

tSin(i) = Sin(xVal(i))

tCos(i) = Cos(xVal(i))

tTan(i) = Tan(xVal(i))

Next i

Dim As Integer rSin, rCos, rTan, rTrig

Dim Shared rhot(rTrig, n2) As Double

For i As Integer = 0 To rTrig

For j As Integer = 0 To n2

rhot(i, j) = NaN

Next j

Next i

Function rho(x() As Double, y() As Double, Byval rdx As Integer, _

Byval i As Integer, Byval n As Integer) As Double

If n < 0 Then Return 0

If n = 0 Then Return y(i+1)

Dim As Integer idx = (n1 - 1 - n) * (n1 - n) / 2 + i + 1

If rhot(rdx, idx) = NaN Then 'valor aún no calculado

rhot(rdx, idx) = (x(i+1) - x(i+1 + n)) _

/ (rho(x(), y(), rdx, i, n-1) - rho(x(), y(), rdx, i+1, n-1)) _

+ rho(x(), y(), rdx, i+1, n-2)

End If

Return rhot(rdx, idx)

End Function

Function thieleInterpolator(x() As Double, y() As Double, _

Byval rdx As Integer, Byval xin As Double, Byval n As Integer) As Double

If n > n1-1 Then Return 1

Return rho(x(), y(), rdx, 0, n) - rho(x(), y(), rdx, 0, n-2) _

+ (xin-x(n+1)) / thieleInterpolator(x(), y(), rdx, xin, n+1)

End Function

Print " PI : 3.141592653589793"

Print " 6*arcsin(0.5) : "; 6 * Asin(0.5)

Print " 3*arccos(0.5) : "; 3 * Acos(0.5)

Print " 4*arctan(1.0) : "; 4 * Atn(1.0)

Print "6*thiele(tSin,xVal,rSin,0.5,0) : "; 6 * thieleInterpolator(tSin(), xVal(), rSin, 0.5, 0)

Print "3*thiele(tCos,xVal,rCos,0.5,0) : "; 3 * thieleInterpolator(tCos(), xVal(), rCos, 0.5, 0)

Print "4*thiele(tTan,xVal,rTan,1.0,0) : "; 4 * thieleInterpolator(tTan(), xVal(), rTan, 1.0, 0)

Sleep</syntaxhighlight>