XProfan

WindowTitle mkstr$(" ",26)+upper$("AKIMA-Interpolation mittels Kubischer Splines [H. Akima, 1978]")
'  Portiert nach XProfan11 (CL) CopyLeft 2014-08 by P. Specht, Vienna (Austria).
'{ Keine wie auch immer geartete Garantie. Nutzung auf alleinige Gefahr des Anwenders!
'
' Übersetzungsbasis: Turbo Pascal Handbuch 1974; Für Detailklärungen wurde verwendet:
' F90/95: https://www.tacc.utexas.edu/documents/13601/162125/fortran_class.pdf
' F95: https://www.mrao.cam.ac.uk/~rachael/compphys/SelfStudyF95.pdf
'
' Pascal-DemoSource: https://jean-pierre.moreau.pagesperso-orange.fr/Pascal/akima_pas.txt
' Orig. 2D-Algorithmus: Hiroshi Akima,U.S. Department of Commerce, NTIA/ITS, F90-Version of 1995/08,
' veröffentlicht in "Hiroshi Akima, 'A Method of Bivariate Interpolation and Smooth Surface Fitting
' for Irregularly Distributed Data Points", ACM Transactions on Math. Software Vol.4 No.2, June 1978
' pp. 148-159. Copyright 1978, Association for Computing Machinery, Inc.
' Fortran-90: TOMS 22.yr,No3 (Sept. 1996) p.357ff,  https://calgo.acm.org/760.gz (Abfrg. 2014-08)
' Test Data: https://cran.r-project.org/web/packages/akima/akima.pdf
'
'********************************************************
'*   Program to demonstrate the Akima spline fitting    *
'*       of Function SIN(X) in double precision         *
'* ---------------------------------------------------- *
'*   Reference: Di base Scientific Subroutines, Vol. II   *
'*   By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1].  *
'*                                                      *
'*             Pascal Version by J.-P. Moreau, Paris.   *
'*                       (www.jpmoreau.fr)              *
'* ---------------------------------------------------- *
'* SAMPLE RUN:                                          *
'*                                                      *
'* Akima spline fitting of SIN(X):                      *
'*                                                      *
'*   X   SIN(X) HANDBOOK   AKIMA INTERPOLATION   ERROR  *
'* ---------------------------------------------------- *
'* 0.00    0.0000000          0.0000000       0.0000000 *
'* 0.05    0.0499792          0.0500402      -0.0000610 *
'* 0.10    0.0998334          0.0998435      -0.0000101 *
'* 0.15    0.1494381          0.1494310       0.0000072 *
'* 0.20    0.1986693          0.1986459       0.0000235 *
'* 0.25    0.2474040          0.2474157      -0.0000118 *
'* 0.30    0.2955202          0.2955218      -0.0000016 *
'* 0.35    0.3428978          0.3428916       0.0000062 *
'* 0.40    0.3894183          0.3894265      -0.0000081 *
'* 0.45    0.4349655          0.4349655       0.0000000 *
'* 0.50    0.4794255          0.4794204       0.0000051 *
'* 0.55    0.5226872          0.5226894      -0.0000021 *
'* 0.60    0.5646425          0.5646493      -0.0000068 *
'* 0.65    0.6051864          0.6051821       0.0000043 *
'* 0.70    0.6442177          0.6442141       0.0000035 *
'* 0.75    0.6816388          0.6816405      -0.0000017 *
'* ---------------------------------------------------- *
'*                                                      *
'*******************************************************}
'}

PROC Interpol_Akima

    '********************************************************
    '*       Akima Spline Fitting Subroutine                *
    '* -----------------------------------------------------*
    '* The input table is X(i), Y(i), where Y(i) is the     *
    '* dependant variable. The interpolation point is xx,   *
    '* which is assumed to be in the interval of the table  *
    '* with at least one table value to the left, and three *
    '* to the right. The interpolated returned value is yy. *
    '* n is returned as an error check (n=0 implies error). *
    '* It is also assumed that the X(i) are in ascending    *
    '* order!                                               *
    '********************************************************
    parameters i&
    n&=1
    ' {special case xx=0}

    if xx!=0: yy!=0:return:endif

        ' {Check to see if interpolation point is correct}

        if (xx!<X![1]) or (xx!>=X![iv&-3]): n&=0:return:endif

            X![0]=2*X![1]-X![2]
            ' {Calculate Akima coefficients, a and b}

            whileloop iv&-1:i&=&Loop

                '   {Shift i to i+2}
                XM![i&+2]=(Y![i&+1]-Y![i&])/(X![i&+1]-X![i&])

            endwhile

            XM![iv&+2]=2*XM![iv&+1]-XM![iv&]
            XM![iv&+3]=2*XM![iv&+2]-XM![iv&+1]
            XM![2]=2*XM![3]-XM![4]
            XM![1]=2*XM![2]-XM![3]

            whileloop iv&:i&=&Loop

                a!=ABS(XM![i&+3]-XM![i&+2])
                b!=ABS(XM![i&+1]-XM![i&])

                if (a!+b!)<>0

                    Z![i&]=(a!*XM![i&+1]+b!*XM![i&+2])/(a!+b!)

                else

                    Z![i&]=(XM![i&+2]+XM![i&+1])/2

                endif

            endwhile

            ' {Find relevant table interval}
            i&=0

            repeat

                inc i&

            until xx!<=X![i&]

            dec i&
            ' {Begin interpolation}
            b!=X![i&+1]-X![i&]
            a!=xx!-X![i&]
            yy!=Y![i&]+Z![i&]*a!+(3*XM![i&+2]-2*Z![i&]-Z![i&+1])*a!*a!/b!
            yy!=yy!+(Z![i&]+Z![i&+1]-2*XM![i&+2])*a!*a!*a!/(b!*b!)

        EndProc

        '{ Main Program AKIMA-SINTEST
        Windowstyle 24:font 0:set("decimals",8)
        CLS'Window 0,0-%maxx,%maxy
        var fmt$=" 0.0000000;-0.0000000"
        var SIZE& = 14
        declare X![SIZE&+1],Y![SIZE&+1],XM![SIZE&+4],Z![SIZE&+1]
        declare a!,b!,xx!,yy!, i&,iv&,n&
        iv&=14'{Number of points in table}
        CLS
        Font 2:Print "\n";tab(10);"  Hiroshi Akima's Interpolation mit kubischen Splines"
        Font 0:print tab(10);" "+mkstr$("-",53)
        print tab(10);" Basis sind 14 unregelmäßig verteilte Stützwerte einer "
        print tab(10);" Eigenbau-Sinusfunktion. Errechnet wird der absolute   "
        Print tab(10);" Fehler (als Differenz zur internen Sin(x)-Funktion).  "
        print tab(10);" Im Programmtext sind die Vergleichswerte zu finden!   \n"
        ' -----------------------------------------------------------------
        ' Input y=Sine(x) table
        ' -----------------------------------------------------------------
        ' Sine table values from  Handbook of mathematical functions
        ' by M. Abramowitz and I.A. Stegun, NBS, june 1964
        ' -----------------------------------------------------------------
        X![1]=0.000:Y![1]=0.00000000  :  X![2]=0.125:Y![2]=0.12467473
        X![3]=0.217:Y![3]=0.21530095  :  X![4]=0.299:Y![4]=0.29456472
        X![5]=0.376:Y![5]=0.36720285  :  X![6]=0.450:Y![6]=0.43496553
        X![7]=0.520:Y![7]=0.49688014  :  X![8]=0.589:Y![8]=0.55552980
        X![9]=0.656:Y![9]=0.60995199  :  X![10]=0.721:Y![10]=0.66013615
        X![11]=0.7853981634           : Y![11]=0.7071067812
        X![12]=0.849:Y![12]=0.75062005:  X![13]=0.911:Y![13]=0.79011709
        X![14]=0.972:Y![14]=0.82601466
        ' -----------------------------------------------------------------
        print "\n";mkstr$("-",77)
        font 2
        print tab(3+3);"X";tab(%pos+15);"SIN(X)";\
        tab(%pos+3);"AKIMA-Stützwertinterpolation";tab(%pos+2);"ABS.FEHLER"
        font 0
        print mkstr$("=",77)
        ' {main loop}
        xx!=0

        whileloop 16:i&=&Loop

            Interpol_Akima(i&)
            print tab(4);xx!;tab(%pos+6);SIN(xx!);tab(%pos+9);yy!;tab(%pos+11);format$(fmt$,SIN(xx!) - yy!)
            casenot i& mod 4:print
            xx! = xx! + 0.05

        endwhile

        ' {print footer}
        print mkstr$("-",77)
        beep
        waitinput
        END
        '} of file akima.prf