| home | send comment | send link | add bookmark |

Interpolation using Bessel's formula

by Stephen R. Schmitt

Introduction

Interpolation lets you estimate the value of a process between two known tabulated values. Suppose that values of a process y(t) are tabulated at tn with a constant interval of dt = tn+1 - tn. To find the value of y(tn + u·dt), where 

0 < u < 1,
linear interpolation gives 
y = y(tn) + u·[y(tn+1) - y(tn)]
Often, this is good enough. Sometimes greater accuracy is needed. Given the tabulated values of a process y(t) at six evenly spaced points t1 ... t6; we want to evaluate y(t) at t3 + u·dt, where dx is the interval between samples. In this situation, we can use nonlinear or polynomial interpolation to obtain a more precise value for y(t). A theorem attributed to Weierstrass, which states that every continuous function in an interval (a,b) can be represented in that interval to any desired accuracy by a polynomial, justifies the use of polynomial interpolation. Several methods work, including:

These methods all use a table of differences for equal-interval interpolation. The first two columns are tn and y(tn). The entries in the remaining columns are the differences between the two entries in the column immediately to the left. Bessel's finite difference interpolation uses a triangular table with entries as shown below: 

t1    y1   
           dA1.5
t2    y2          dB2
           dA2.5        dC2.5
t3    y3          dB3         dD3
           dA3.5        dC3.5       dE3.5
t4    y4          dB4         dD4
           dA4.5        dC4.5
t5    y5          dB5
           dA5.5  
t6    y6
For the above table, with six samples, Bessel's formula is: 
y =  0.5*( y3 + y4 ) + B1(u)*dA3.5 
  + B2(u)*(dB3 + dB4) + B3(u)*dC3.5 
  + B4(u)*(dD3 + dD4) + B5(u)*dE3.5
Where 
u = (t - t3) / (t4 - t3)
Bessel's coefficients defined as follows: 
        u - 1/2           (u + n/2 - 3/2)!
Bn(u) = ------- · ---------------------------------     n is odd
           n      (u + n/2 - 3/2 - n + 1)! (n - 1)!

         1        (u + n/2 - 1)!
Bn(u) = --- · ---------------------                     n is even
         2    (u + n/2 - 1 - n)! n!
The first few coefficients are: 
B1(u) = u - 0.5 
B2(u) = 0.5·u·(u - 1)/2 
B3(u) = (u - 0.5)·u·(u - 1)/6
B4(u) = 0.5·(u + 1)·u·(u - 1)·(u - 2)/24
B5(u) = (u - 0.5)·(u + 1)·u·(u - 1)·(u - 2)/120
B6(u) = 0.5·(u + 2)·(u + 1)·u·(u - 1)·(u - 2)·(u - 3)/720
Interpolation formula taken as far as B1 is identical to linear interpolation. The addition of higher terms fits a curve to more and more points on both sides of the desired value to give a more accurate value for y(t).

Below is a sample program showing an implementation of Bessel interpolation to calculate the position of the moon using positions tabulated daily.

Zeno source code

Zeno 1.2 is an interpreter for the Zeno programming language. It is an easy to learn and is suitable for educational purposes. 

type table : array[6] of real

program

    var ras, dec, dst : table

    % Moon : Apparent Geocentric Positions, True Equator and Equinox of Date
    % UT1  :       h  m   s       h  m   s           o  '   "           km          
    % ------------------------------------------------------------------------      
    % 2004 Apr 17 00:00:00.0     23 57 08.676     -  4 16 29.27     389306.996      
    % 2004 Apr 18 00:00:00.0      0 42 27.849     +  1 37 35.17     392652.502      
    % 2004 Apr 19 00:00:00.0      1 27 20.166     +  7 22 35.29     395849.635      
    % 2004 Apr 20 00:00:00.0      2 12 38.939     + 12 45 44.69     398812.496      
    % 2004 Apr 21 00:00:00.0      2 59 08.486     + 17 34 58.11     401416.265      
    % 2004 Apr 22 00:00:00.0      3 47 20.082     + 21 38 41.48     403502.112

    ras[1] := ras2real( 23, 57, 08.676 ) - 24    % due to discontinuity below
    ras[2] := ras2real(  0, 42, 27.849 )
    ras[3] := ras2real(  1, 27, 20.166 )
    % result between element 3 and 4
    ras[4] := ras2real(  2, 12, 38.939 )
    ras[5] := ras2real(  2, 59, 08.486 )
    ras[6] := ras2real(  3, 47, 20.082 )

    dec[1] := dec2real( - 4, 16, 29.27)
    dec[2] := dec2real( + 1, 37, 35.17)
    dec[3] := dec2real( + 7, 22, 35.29)
    % result between element 3 and 4
    dec[4] := dec2real( +12, 45, 44.69)
    dec[5] := dec2real( +17, 34, 58.11)
    dec[6] := dec2real( +21, 38, 41.48)

    dst[1] := 389306.996      
    dst[2] := 392652.502      
    dst[3] := 395849.635      
    % result between element 3 and 4
    dst[4] := 398812.496      
    dst[5] := 401416.265      
    dst[6] := 403502.112      

    var y : real
    var u : real := time2real( 18, 0, 0 )
    
    y := interpolate( ras, u )
    put "RA =", ha2str( y )...
    y := interpolate( dec, u )
    put ", Dec =", dec2str( y )...
    y := interpolate( dst, u )
    put ", Dist =", y:12:3, " km"

end program

% interpolate from six element table using Bessel's formula
%  u - fraction between elements 3 and 4 of table
%  y - table of six values assumed to be spaced evenly and to be precise
function interpolate( var y : table, u : real ) : real

    var rv : real

    % first level of divided differences
    var dA15 : real := y[2] - y[1]  
    var dA25 : real := y[3] - y[2]  
    var dA35 : real := y[4] - y[3] 
    var dA45 : real := y[5] - y[4]  
    var dA55 : real := y[6] - y[5]  

    % second level
    var dB2 : real := dA25 - dA15
    var dB3 : real := dA35 - dA25
    var dB4 : real := dA45 - dA35
    var dB5 : real := dA55 - dA45

    % third level
    var dC25 : real := dB3 - dB2
    var dC35 : real := dB4 - dB3
    var dC45 : real := dB5 - dB4

    % fourth level
    var dD3 : real := dC35 - dC25
    var dD4 : real := dC45 - dC35

    % fifth level
    var dE35 : real := dD4 - dD3

    rv := 0.5*(y[3] + y[4]) + B1(u)*dA35 + B2(u)*(dB3 + dB4) 
        + B3(u)*dC35 + B4(u)*(dD3 + dD4) + B5(u)*dE35
    
    return rv

end function

% first bessel coefficient
function B1( u : real ) : real

    return (u - 0.5)

end function

% second bessel coefficient
function B2( u : real ) : real

    return 0.5*u*(u - 1)/2

end function

% third bessel coefficient
function B3( u : real ) : real

    return (u - 0.5)*u*(u - 1)/6

end function

% fourth bessel coefficient
function B4( u : real ) : real

    return 0.5*(u + 1)*u*(u - 1)*(u - 2)/24

end function

% fifth bessel coefficient
function B5( u : real ) : real

    return (u - 0.5)*(u + 1)*u*(u - 1)*(u - 2)/120

end function

% convert (hr min sec) to fraction of day
function time2real( hr, min, sec : int ) : real

    return (hr + min/60 + sec/3600)/24 

end function

% convert (hr min sec) to right ascension
function ras2real( hr, min : int, sec : real ) : real

    return hr + min/60 + sec/3600

end function

% convert (deg min sec) to declination
function dec2real( deg, min : int, sec : real ) : real

    var rv : real
    
    if deg > 0 then
        rv := deg + min/60 + sec/3600
    else
        rv := -( -deg + min/60 + sec/3600 )
    end if

    return rv

end function

% converts hour angle into hour angle string
function ha2str( ra : real ) : string

    assert (0 <= ra) and (ra < 24)       % assure valid range
    var h : int := floor( ra )
    ra := ra - h
    var m : int := floor( 60*ra )
    ra := ra - m / 60
    var s : real := 3600*ra
    return intstr( h, 3 ) & "h " & intstr( m, 2 ) & "m " & frealstr( s, 5, 3 ) & "s"

end function

% converts declination angle in degrees into string
function dec2str( x : real ) : string

    assert (-90 <= x) and (x <= +90)     % assure valid range
    var dec : real := abs( x )           % work with absolute value
    var d : int := floor( dec )
    dec := dec - d
    var m : int := floor( 60*dec )
    dec := dec - m / 60
    var s : real := 3600*dec
    return intstr( sgn(x)*d, 3 ) & "° " & intstr( m, 2 ) & "' " & frealstr( s, 4, 2 ) & "\""

end function

% returns absolute value of real number
function abs( x : real ) : real

    if x < 0 then
        return -x
    end if

    return x

end function

% returns sign of real number as int
function sgn( x : real ) : int

    if x < 0 then
        return -1
    end if

    return +1

end function

Sample output

RA =  2h  1m 14.150s, Dec = 11° 27' 40.49", Dist =  398099.951 km

AbCd Classics - free on-line books

Copyright © 2005, Stephen R. Schmitt