' spline_interp.pb — Spline interpolation (linear, cubic)
#COMPILE EXE
#DIM ALL

FUNCTION PBMAIN() AS LONG
    LOCAL xData(20) AS DOUBLE
    LOCAL yData(20) AS DOUBLE
    LOCAL numPoints AS LONG
    LOCAL i AS LONG
    LOCAL j AS LONG
    LOCAL interpX AS DOUBLE
    LOCAL interpY AS DOUBLE
    LOCAL totalInterp AS LONG
    LOCAL leftIdx AS LONG
    LOCAL rightIdx AS LONG
    LOCAL tVal AS DOUBLE
    LOCAL y0 AS DOUBLE
    LOCAL y1 AS DOUBLE
    LOCAL y2 AS DOUBLE
    LOCAL y3 AS DOUBLE
    LOCAL a0 AS DOUBLE
    LOCAL a1 AS DOUBLE
    LOCAL a2 AS DOUBLE
    LOCAL a3 AS DOUBLE
    LOCAL mu AS DOUBLE
    LOCAL mu2 AS DOUBLE
    LOCAL interpolatedY AS DOUBLE
    LOCAL iter AS LONG
    LOCAL sampleIdx AS LONG
    LOCAL h AS DOUBLE
    LOCAL diffVal AS DOUBLE
    LOCAL maxError AS DOUBLE

    numPoints = 6

    xData(0) = 0.0
    yData(0) = 0.0
    xData(1) = 1.0
    yData(1) = 2.0
    xData(2) = 2.0
    yData(2) = 5.0
    xData(3) = 3.0
    yData(3) = 9.0
    xData(4) = 4.0
    yData(4) = 14.0
    xData(5) = 5.0
    yData(5) = 20.0

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

    PRINT "Linear interpolation:"
    FOR sampleIdx = 0 TO 19
        interpX = sampleIdx * 0.25

        leftIdx = 0
        FOR i = 0 TO numPoints - 2
            IF interpX >= xData(i) AND interpX <= xData(i + 1) THEN
                leftIdx = i
            END IF
        NEXT i
        rightIdx = leftIdx + 1

        tVal = (interpX - xData(leftIdx)) / (xData(rightIdx) - xData(leftIdx))
        interpolatedY = yData(leftIdx) + tVal * (yData(rightIdx) - yData(leftIdx))

        IF sampleIdx / 3 * 3 = sampleIdx THEN
            PRINT "  x=" + STR$(interpX) + " y_linear=" + STR$(interpolatedY)
        END IF
        totalInterp = totalInterp + 1
    NEXT sampleIdx

    PRINT "Cubic (Catmull-Rom) interpolation:"
    FOR sampleIdx = 0 TO 19
        interpX = sampleIdx * 0.25

        leftIdx = 0
        FOR i = 0 TO numPoints - 2
            IF interpX >= xData(i) AND interpX <= xData(i + 1) THEN
                leftIdx = i
            END IF
        NEXT i

        IF leftIdx = 0 THEN leftIdx = 1
        IF leftIdx >= numPoints - 2 THEN leftIdx = numPoints - 3

        y0 = yData(leftIdx - 1)
        y1 = yData(leftIdx)
        y2 = yData(leftIdx + 1)
        y3 = yData(leftIdx + 2)

        mu = (interpX - xData(leftIdx)) / (xData(leftIdx + 1) - xData(leftIdx))
        mu2 = mu * mu

        a0 = 0 - 0.5 * y0 + 1.5 * y1 - 1.5 * y2 + 0.5 * y3
        a1 = y0 - 2.5 * y1 + 2.0 * y2 - 0.5 * y3
        a2 = -0.5 * y0 + 0.5 * y2
        a3 = y1

        interpolatedY = a0 * mu * mu2 + a1 * mu2 + a2 * mu + a3

        IF sampleIdx / 3 * 3 = sampleIdx THEN
            PRINT "  x=" + STR$(interpX) + " y_cubic=" + STR$(interpolatedY)
        END IF
    NEXT sampleIdx

    PRINT "Total interpolations: " + STR$(totalInterp)
    PRINT "=== End Spline ==="
    FUNCTION = 0
END FUNCTION
