' data_splines.pb — Cubic spline interpolation
#COMPILE EXE
#DIM ALL

FUNCTION PBMAIN() AS LONG
    LOCAL xData(20) AS DOUBLE
    LOCAL yData(20) AS DOUBLE
    LOCAL coefficients(20) AS DOUBLE
    LOCAL numPoints AS LONG
    LOCAL i AS LONG
    LOCAL j AS LONG
    LOCAL interpX AS DOUBLE
    LOCAL interpY AS DOUBLE
    LOCAL h AS DOUBLE
    LOCAL aVal AS DOUBLE
    LOCAL bVal AS DOUBLE
    LOCAL cVal AS DOUBLE
    LOCAL dVal AS DOUBLE
    LOCAL tVal AS DOUBLE
    LOCAL diff AS DOUBLE
    LOCAL totalError AS DOUBLE
    LOCAL numInterp AS LONG
    LOCAL iter AS LONG
    LOCAL idx AS LONG
    LOCAL xVal AS DOUBLE
    LOCAL yVal AS DOUBLE

    numPoints = 5

    xData(0) = 0.0
    yData(0) = 0.0
    xData(1) = 1.0
    yData(1) = 1.0
    xData(2) = 2.0
    yData(2) = 4.0
    xData(3) = 3.0
    yData(3) = 9.0
    xData(4) = 4.0
    yData(4) = 16.0

    PRINT "=== Cubic Spline Interpolation ==="
    PRINT "Data points: " + STR$(numPoints)
    FOR i = 0 TO numPoints - 1
        PRINT "  x=" + STR$(xData(i)) + " y=" + STR$(yData(i))
    NEXT i

    FOR i = 0 TO numPoints - 2
        h = xData(i + 1) - xData(i)
        aVal = yData(i)
        bVal = 0
        diff = yData(i + 1) - yData(i)
        cVal = (3.0 * diff / h - 2.0 * bVal) / h
        dVal = (bVal - 2.0 * diff / h) / (h * h)
        coefficients(i) = cVal
    NEXT i

    PRINT "Interpolated values:"
    numInterp = 0
    FOR i = 0 TO numPoints - 2
        idx = i
        aVal = yData(idx)
        h = xData(idx + 1) - xData(idx)
        cVal = coefficients(idx)

        FOR j = 1 TO 4
            tVal = j / 4.0
            interpX = xData(idx) + tVal * h
            interpY = aVal + cVal * (interpX - xData(idx)) * (interpX - xData(idx + 1))
            PRINT "  x=" + STR$(interpX) + " y=" + STR$(interpY)
            numInterp = numInterp + 1
        NEXT j
    NEXT i

    PRINT "Interpolated points: " + STR$(numInterp)
    PRINT "=== End Splines ==="
    FUNCTION = 0
END FUNCTION
