Announcement

Collapse

Forum Guidelines

This forum is for finished source code that is working properly. If you have questions about this or any other source code, please post it in one of the Discussion Forums, not here.
See more
See less

Cubic Spline

Collapse
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • Cubic Spline

    Code:
    FUNCTION CSP_FIND_INDEX(                                             _
                                   N_CSP      AS LONG                    , _
                            BYVAL pU_CSP      AS UDT_CSP        POINTER  , _
                            BYVAL pU_CSP_PARM AS UDT_CSP_PARM   POINTER  , _
                                   t_INTP     AS DOUBLE                     )  AS LONG
    
       '  Routine finds the ith index that an interpolated data point would lie between
       '  Note this ranges 00 to n - 01 - 01, thus returns -1 on fail which is tested
    
       DIM i                                                 AS LONG
    
    
       IF       N_CSP = 00                                   OR _
                pU_CSP = 00                                  OR _
                pU_CSP_PARM = 00                             THEN
    
          FUNCTION = -1
    
       ELSEIF  @pU_CSP_PARM.N_DATA = 00                      OR _     '  Test all allocations and pointers
               @pU_CSP_PARM.pt_DATA = 00                     OR _
               @pU_CSP_PARM.py_DATA = 00                     OR _
               @pU_CSP_PARM.N_INTP = 00                      OR _
               @pU_CSP_PARM.pt_INTP = 00                     OR _
               @pU_CSP_PARM.py_INTP = 00                     OR _
               @pU_CSP_PARM.N_LDU = 00                       OR _
               @pU_CSP_PARM.pU_LDU = 00                      THEN
    
    
       ELSEIF  NOT @pU_CSP_PARM.N_DATA = @pU_CSP_PARM.N_LDU  THEN     '  LDU allocation must match the input data size
    
          FUNCTION = -1
    
       ELSE
    
          '  Trap the boundaries
    
          IF       t_INTP <= @[email protected]_DATA[00]                THEN
    
             FUNCTION = 00
    
          ELSEIF   t_INTP >= @[email protected]_DATA[N_CSP - 01 - 01]   THEN
    
             FUNCTION = N_CSP - 01 - 01
    
          ELSE
    
             FOR i = 00 TO N_CSP - 01 - 01
    
                IF t_INTP >= @[email protected]_DATA[i]                 AND _
                   t_INTP <  @[email protected]_DATA[i + 01]            THEN
    
                   FUNCTION = i
                   EXIT FUNCTION
    
                END IF
    
             NEXT i
    
             '  Fail condition is encountered
             FUNCTION = -1
    
          END IF
    
       END IF
    
    END FUNCTION
    
    
    FUNCTION CSP_SOLVE     (                                             _
                                   N_CSP      AS LONG                    , _
                            BYVAL pU_CSP      AS UDT_CSP        POINTER  , _
                            BYVAL pU_CSP_PARM AS UDT_CSP_PARM   POINTER     )  AS LONG
    
       '  Routine computes a cubic spline solution to a set of input data by pointer
    
       DIM i                                                 AS LONG
       DIM i_INTP                                            AS LONG
    
    
       IF       N_CSP = 00                                   OR _
                pU_CSP = 00                                  OR _
                pU_CSP_PARM = 00                             THEN
    
          FUNCTION = %F
    
       ELSEIF  @pU_CSP_PARM.N_DATA = 00                      OR _     '  Test all allocations and pointers
               @pU_CSP_PARM.pt_DATA = 00                     OR _
               @pU_CSP_PARM.py_DATA = 00                     OR _
               @pU_CSP_PARM.N_INTP = 00                      OR _
               @pU_CSP_PARM.pt_INTP = 00                     OR _
               @pU_CSP_PARM.py_INTP = 00                     OR _
               @pU_CSP_PARM.N_LDU = 00                       OR _
               @pU_CSP_PARM.pU_LDU = 00                      THEN
    
    
       ELSEIF  NOT @pU_CSP_PARM.N_DATA = @pU_CSP_PARM.N_LDU  THEN     '  LDU allocation must match the input data size
    
          FUNCTION = %F
    
       ELSE
    
          '  Compute b and h coefficients
    
          FOR i = 00 TO N_CSP - 01 - 01
    
             @pU_CSP[i].h = @[email protected]_DATA[i + 01] - @[email protected]_DATA[i]
    
             IF NOT @pU_CSP[i].h = 0.0#                      THEN
    
                @pU_CSP[i].b = 1.0# / @pU_CSP[i].h * (@[email protected]_DATA[i + 01] - @[email protected]_DATA[i])
    
             END IF
    
          NEXT i
    
    
          '  Compute v and u coefficients
    
          FOR i = 01 TO N_CSP - 01 - 01
    
             @pU_CSP[i].v = 2.0# * (@pU_CSP[i - 01].h + @pU_CSP[i].h)
             @pU_CSP[i].u = 6.0# * (@pU_CSP[i].b - @pU_CSP[i - 01].b)
    
          NEXT i
    
    
          '  Load the LDU coefficients to prepare for the tridiagonal solution
    
          FOR i = 00 TO @pU_CSP_PARM.N_LDU - 01
    
             '  Load the LDU vectors
             '  Note the IIF operator for the boundary condition (i = 00 OR i = N_CSP - 01) which processes S''(x) = 0 by setting L, U = 0 and D = 1
             @[email protected]_LDU[i].L = IIF(i = 00 OR i = N_CSP - 01, 0.0#, IIF(i = 01             , %F, @pU_CSP[i - 01].h))
             @[email protected]_LDU[i].D = IIF(i = 00 OR i = N_CSP - 01, 1.0#, @pU_CSP[i].v                                   )
             @[email protected]_LDU[i].U = IIF(i = 00 OR i = N_CSP - 01, 0.0#, IIF(i = N_CSP - 01 - 01, %F, @pU_CSP[i].h)     )
    
             '  Load the b solution source vector which in this case is the [u] vector
             @[email protected]_LDU[i].b = @pU_CSP[i].u
    
          NEXT i
    
    
          '  Solve the cubic spline [z] vector for S.i''(x) using the tridiagonal Gauss Seidel LDU solver
    
          GS_LDU_SOLVE   @pU_CSP_PARM.N_LDU, @pU_CSP_PARM.pU_LDU, VARPTR(@pU_CSP_PARM.U_LDU_PARM)
          #IF %F
          MSGBOX GS_LDU_GET_S(@pU_CSP_PARM.N_LDU, @pU_CSP_PARM.pU_LDU), %MB_ICONINFORMATION OR %MB_TASKMODAL OR %MB_TOPMOST, $S_MSGBOX_HDR
          #ENDIF
    
          '  Range over the span of the interpolated t points
    
          FOR i_INTP = 00 TO @pU_CSP_PARM.N_INTP - 01
    
             '  Create the uniform interpolation points
             @[email protected]_INTP[i_INTP] = @[email protected]_DATA[00] + (@[email protected]_DATA[@pU_CSP_PARM.N_DATA - 01] - @[email protected]_DATA[00]) * i_INTP / (@pU_CSP_PARM.N_INTP - 01)
    
             '  Find the data index that each point lies in
             CSP_FIND_INDEX    N_CSP                            , _
                               pU_CSP                           , _
                               pU_CSP_PARM                      , _
                               @[email protected]_INTP[i_INTP]    TO i
    
             '  Error check the index for -1 from the find index function
             IF NOT i = -1                                   THEN
    
                '  Compute the cubic spline interpolation y = ao + a1 t + a2 t^2 + a3 t^3
                '  Note the LDU solved [x]kp1 term replaces the [z] solution vector by variable name
                @[email protected]_INTP[i_INTP] = @[email protected]_LDU[i + 01].xkp1 / (6.0# * @pU_CSP[i].h) * (@[email protected]_INTP[i_INTP] - @[email protected]_DATA[i]) ^ 3                                                  + _
                                                @[email protected]_LDU[i].xkp1      / (6.0# * @pU_CSP[i].h) * (@[email protected]_DATA[i + 01] - @[email protected]_INTP[i_INTP]) ^ 3                                             + _
                                                (@[email protected]_DATA[i + 01] / @pU_CSP[i].h - @[email protected]_LDU[i + 01].xkp1 / 6.0# * @pU_CSP[i].h) * (@[email protected]_INTP[i_INTP] - @[email protected]_DATA[i])       + _
                                                (@[email protected]_DATA[i]      / @pU_CSP[i].h - @[email protected]_LDU[i].xkp1      / 6.0# * @pU_CSP[i].h) * (@[email protected]_DATA[i + 01] - @[email protected]_INTP[i_INTP])
    
             END IF
    
          NEXT i_INTP
    
          FUNCTION = %T
    
       END IF
    
    END FUNCTION
    
    
    FUNCTION CSP_SET_CLIP  (                                             _
                                   N_CSP      AS LONG                    , _
                            BYVAL pU_CSP      AS UDT_CSP        POINTER  , _
                            BYVAL pU_CSP_PARM AS UDT_CSP_PARM   POINTER     )  AS LONG
    
       '  Routine sets cubic spline solution to the clipboard for easy plotting
    
       DIM i                                                 AS LONG
       DIM i_INTP                                            AS LONG
    
       DIM S_CLIP                                            AS STRING
    
    
       IF       N_CSP = 00                                   OR _
                pU_CSP = 00                                  OR _
                pU_CSP_PARM = 00                             THEN
    
          FUNCTION = %F
    
       ELSEIF  @pU_CSP_PARM.N_DATA = 00                      OR _     '  Test all allocations and pointers
               @pU_CSP_PARM.pt_DATA = 00                     OR _
               @pU_CSP_PARM.py_DATA = 00                     OR _
               @pU_CSP_PARM.N_INTP = 00                      OR _
               @pU_CSP_PARM.pt_INTP = 00                     OR _
               @pU_CSP_PARM.py_INTP = 00                     OR _
               @pU_CSP_PARM.N_LDU = 00                       OR _
               @pU_CSP_PARM.pU_LDU = 00                      THEN
    
    
       ELSEIF  NOT @pU_CSP_PARM.N_DATA = @pU_CSP_PARM.N_LDU  THEN     '  LDU allocation must match the input data size
    
          FUNCTION = %F
    
       ELSE
    
          '  Range over the span of the interpolated t points
    
          FOR i_INTP = 00 TO @pU_CSP_PARM.N_INTP - 01
    
             S_CLIP  += STR$(@[email protected]_INTP[i_INTP]) + " " + STR$(@[email protected]_INTP[i_INTP]) + $CRLF
    
          NEXT i_INTP
    
          CLIPBOARD SET TEXT S_CLIP
    
          FUNCTION = %T
    
       END IF
    
    END FUNCTION
    
    
    FUNCTION CSP_GET_S     (                                             _
                                   N_CSP      AS LONG                    , _
                            BYVAL pU_CSP      AS UDT_CSP        POINTER  , _
                            BYVAL pU_CSP_PARM AS UDT_CSP_PARM   POINTER     )  AS STRING
    
       DIM i                                                 AS LONG
       DIM i_INTP                                            AS LONG
       DIM S_CSP                                             AS STRING
    
    
       IF       N_CSP = 00                                   OR _
                pU_CSP = 00                                  OR _
                pU_CSP_PARM = 00                             THEN
    
          FUNCTION = ""
    
       ELSEIF  @pU_CSP_PARM.N_DATA = 00                      OR _     '  Test all allocations and pointers
               @pU_CSP_PARM.pt_DATA = 00                     OR _
               @pU_CSP_PARM.py_DATA = 00                     OR _
               @pU_CSP_PARM.N_INTP = 00                      OR _
               @pU_CSP_PARM.pt_INTP = 00                     OR _
               @pU_CSP_PARM.py_INTP = 00                     OR _
               @pU_CSP_PARM.N_LDU = 00                       OR _
               @pU_CSP_PARM.pU_LDU = 00                      THEN
    
    
       ELSEIF  NOT @pU_CSP_PARM.N_DATA = @pU_CSP_PARM.N_LDU  THEN     '  LDU allocation must match the input data size
    
          FUNCTION = ""
    
       ELSE
    
          #IF 00
          FOR i = 00 TO N_CSP - 01
    
             S_CSP   += TRIM$("b[" + STR$(i+1) + "]") + " " + STR$(@pU_CSP[i].b) + " " + _
                        TRIM$("v[" + STR$(i+1) + "]") + " " + STR$(@pU_CSP[i].v) + " " + _
                        TRIM$("u[" + STR$(i+1) + "]") + " " + STR$(@pU_CSP[i].u) + $CRLF
    
          NEXT i
          #ENDIF
    
          #IF 00
          FOR i = 00 TO @pU_CSP_PARM.N_LDU - 01
    
             S_CSP   += TRIM$("L[" + STR$(i+1) + "]") + " " + STR$(@[email protected]_LDU[i].L) + " " + _
                        TRIM$("D[" + STR$(i+1) + "]") + " " + STR$(@[email protected]_LDU[i].D) + " " + _
                        TRIM$("U[" + STR$(i+1) + "]") + " " + STR$(@[email protected]_LDU[i].U) + $CRLF
    
          NEXT i
          #ENDIF
    
          #IF -1
          FOR i_INTP = 00 TO @pU_CSP_PARM.N_INTP - 01
    
             S_CSP   += "[i INTP, t_INTP, y_INTP] " + STR$(i_INTP) + " " + STR$(@[email protected]_INTP[i_INTP]) + " " + STR$(@[email protected]_INTP[i_INTP]) + $CRLF
    
          NEXT i_INTP
          #ENDIF
    
          FUNCTION = S_CSP
    
       END IF
    
    END FUNCTION
    
    
    FUNCTION CSP_TEST      (                                          )  AS LONG
    
       '  Routine computes an cubic spline test function
    
       DIM N_CSP                                             AS LONG
       DIM U_CSP()                                           AS UDT_CSP
       DIM U_CSP_PARM                                        AS UDT_CSP_PARM
    
       DIM N_DATA                                            AS LONG
       DIM t_DATA()                                          AS DOUBLE
       DIM y_DATA()                                          AS DOUBLE
    
       DIM N_INTP                                            AS LONG
       DIM t_INTP()                                          AS DOUBLE
       DIM y_INTP()                                          AS DOUBLE
    
       DIM N_LDU                                             AS LONG
       DIM U_LDU()                                           AS UDT_LDU
    
    
       '  Create a test input data set
       '  Note the assumption is that the load data is sorted
       '     0.9   0
       '     1.3   1.5
       '     1.9   1.85
       '     2.1   2.1
       N_DATA = 04
       REDIM t_DATA(00 TO N_DATA - 01) : ARRAY ASSIGN t_DATA() = 0.9# , 1.3# , 1.9# , 2.1#
       REDIM y_DATA(00 TO N_DATA - 01) : ARRAY ASSIGN y_DATA() = 0.0# , 1.5# , 1.85#, 2.1#
       U_CSP_PARM.N_DATA    = N_DATA
       U_CSP_PARM.pt_DATA   = VARPTR(t_DATA(00))
       U_CSP_PARM.py_DATA   = VARPTR(y_DATA(00))
    
       '  Allocate the LDU solution data structure
       N_LDU = N_DATA
       REDIM U_LDU(00 TO N_LDU - 01)
       U_CSP_PARM.N_LDU     = N_DATA
       U_CSP_PARM.pU_LDU    = VARPTR(U_LDU(00))
    
       '  Allocate the cubic spline coefficients to the number of data points
       N_CSP = N_DATA
       REDIM U_CSP(00 TO N_CSP - 01)
    
       '  Set the LDU solve parameters
       U_CSP_PARM.U_LDU_PARM.N_ITR_MAX = 10
       U_CSP_PARM.U_LDU_PARM.e_ITR_MIN = 1.0E-6#
    
       '  Create a test input interpolation set
       '  Note this is larger than the input data set to create smoothness
       N_INTP = 10
       REDIM t_INTP(00 TO N_INTP - 01)
       REDIM y_INTP(00 TO N_INTP - 01)
       U_CSP_PARM.N_INTP    = N_INTP
       U_CSP_PARM.pt_INTP   = VARPTR(t_INTP(00))
       U_CSP_PARM.py_INTP   = VARPTR(y_INTP(00))
    
    
       '  Solve for cubic spline
       CSP_SOLVE      N_CSP, VARPTR(U_CSP(00)), VARPTR(U_CSP_PARM)
       CSP_SET_CLIP   N_CSP, VARPTR(U_CSP(00)), VARPTR(U_CSP_PARM)
    
       '  Show the interpolation
       #IF -1
       MSGBOX      CSP_GET_S(N_CSP, VARPTR(U_CSP(00)), VARPTR(U_CSP_PARM)), %MB_ICONINFORMATION OR %MB_TASKMODAL OR %MB_TOPMOST, $S_MSGBOX_HDR
       #ENDIF
    
    END FUNCTION

  • #2
    Code:
    #COMPILE EXE "Cubic Spline.exe"
    
    #DIM ALL
    
       $S_MSGBOX_HDR = "Cubic Spline"
       %T = -1
       %F = 00
    
       '  Supports tridiagonal matrix solution
    
       TYPE UDT_LDU
    
          '  Lower, diagonal and upper vectors
    
          L                                                     AS DOUBLE
          D                                                     AS DOUBLE
          U                                                     AS DOUBLE
    
          '  Solution vector x<k> and x<k+1>
          xk                                                    AS DOUBLE
          xkp1                                                  AS DOUBLE
    
          '  Right hand side source vector
          b                                                     AS DOUBLE
    
       END TYPE
    
       '  Supports tridiagonal matrix solver parameters
    
       TYPE UDT_LDU_PARM
    
          N_ITR_MAX                                             AS LONG
          e_ITR_MIN                                             AS DOUBLE
    
       END TYPE
    
    
       '  Supports cubic spline interpolation
    
       TYPE UDT_CSP
    
          '  Spline coefficients
          h                                                     AS DOUBLE
          b                                                     AS DOUBLE
          v                                                     AS DOUBLE
    
          '  Source vector
          u                                                     AS DOUBLE
    
          '  Solution state
          z                                                     AS DOUBLE
    
       END TYPE
    
       '  Supports cubic spline input parameters
    
       TYPE UDT_CSP_PARM
    
          '  References the input data
    
          N_DATA                                                AS LONG
         pt_DATA                                                AS DOUBLE POINTER
         py_DATA                                                AS DOUBLE POINTER
    
          '  References the interpolated data which forms the spline
          N_INTP                                                AS LONG
         pt_INTP                                                AS DOUBLE POINTER
         py_INTP                                                AS DOUBLE POINTER
    
          '  Supports tridiagonal coefficient solution
    
          N_LDU                                                 AS LONG
         pU_LDU                                                 AS UDT_LDU POINTER
          U_LDU_PARM                                            AS UDT_LDU_PARM
    
       END TYPE
    
    
    FUNCTION PBMAIN () AS LONG
    
       '  Call the test function
       '  Note this sets the interpolated data to the clipboard
       CSP_TEST
    
    END FUNCTION
    
    
    FUNCTION GS_LDU_xkp1   (                                       _
                                   N_LDU      AS LONG              , _
                            BYVAL pU_LDU      AS UDT_LDU POINTER      )  AS LONG
    
       '  Routine computes a Gauss Seidel single step iteration
    
       DIM i                                                 AS LONG
    
       DIM Lxk                                               AS DOUBLE
       DIM Uxk                                               AS DOUBLE
    
       IF       N_LDU = 00                                   OR _
                pU_LDU = 00                                  THEN
    
          FUNCTION = %F
    
       ELSE
    
          FOR i = 00 TO N_LDU - 01
    
             '  Form the matrix (L+U)[x.k] noting boundaries for lower and upper
             Lxk = IIF(i = 00        , %F, @pU_LDU[i].L * @pU_LDU[i - 01].xkp1)   '  Uses most advance computed [x], Gauss Seidel method
             Uxk = IIF(i = N_LDU - 01, %F, @pU_LDU[i].U * @pU_LDU[i + 01].xk)
    
             '  Solve [x.kp1]
             IF NOT @pU_LDU[i].D = 0.0#                      THEN
    
                @pU_LDU[i].xkp1 = 1.0# / @pU_LDU[i].D * (@pU_LDU[i].b - (Lxk + Uxk))
    
             END IF
    
          NEXT i
    
          FUNCTION = %T
    
       END IF
    
    END FUNCTION
    
    
    FUNCTION GS_LDU_BFR    (                                       _
                                   N_LDU      AS LONG              , _
                            BYVAL pU_LDU      AS UDT_LDU POINTER      )  AS LONG
    
       '  Routine computes a Gauss Seidel solution vector buffering step
    
       DIM i                                                 AS LONG
    
       IF       N_LDU = 00                                   OR _
                pU_LDU = 00                                  THEN
    
          FUNCTION = %F
    
       ELSE
    
          FOR i = 00 TO N_LDU - 01
    
             @pU_LDU[i].xk = @pU_LDU[i].xkp1
    
          NEXT i
    
          FUNCTION = %T
    
       END IF
    
    END FUNCTION
    
    
    FUNCTION GS_LDU_e_ITR  (                                       _
                                   N_LDU      AS LONG              , _
                            BYVAL pU_LDU      AS UDT_LDU POINTER      )  AS DOUBLE
    
       '  Routine computes a Gauss Seidel solution vector error norm
    
       DIM i                                                 AS LONG
    
       DIM e_ITR                                             AS DOUBLE
    
    
       IF       N_LDU = 00                                   OR _
                pU_LDU = 00                                  THEN
    
          FUNCTION = %F
    
       ELSE
    
          FOR i = 00 TO N_LDU - 01
    
             e_ITR += ABS(@pU_LDU[i].xkp1 - @pU_LDU[i].xk)
    
          NEXT i
    
          FUNCTION = e_ITR
    
       END IF
    
    END FUNCTION
    
    
    FUNCTION GS_LDU_SOLVE  (                                             _
                                   N_LDU      AS LONG                    , _
                            BYVAL pU_LDU      AS UDT_LDU        POINTER  , _
                            BYVAL pU_LDU_PARM AS UDT_LDU_PARM   POINTER     )  AS LONG
    
       '  Routine computes a Gauss Seidel solution to the GS LDU parameters
    
       DIM k_ITR                                             AS LONG
       DIM e_ITR                                             AS DOUBLE
    
    
       IF       N_LDU = 00                                   OR _
                pU_LDU = 00                                  OR _
                pU_LDU_PARM = 00                             THEN
    
          FUNCTION = %F
    
       ELSE
    
          DO
    
             INCR k_ITR
             GS_LDU_xkp1    N_LDU, pU_LDU              '  Comute single step
             GS_LDU_e_ITR   N_LDU, pU_LDU TO e_ITR     '  Compute error norm
             GS_LDU_BFR     N_LDU, pU_LDU              '  Back buffer [x.k] = [x.kp1]
    
             #IF %F
             MSGBOX "k_ITR, e_ITR " + STR$(k_itr) + " " + STR$(e_itr) + $CRLF + GS_LDU_GET_S(N_LDU, pU_LDU), %MB_ICONINFORMATION OR %MB_TASKMODAL OR %MB_TOPMOST, $S_MSGBOX_HDR
             #ENDIF
    
          LOOP WHILE (k_ITR < @pU_LDU_PARM.N_ITR_MAX) AND (e_ITR > @pU_LDU_PARM.e_ITR_MIN)
    
          FUNCTION = %T
    
       END IF
    
    END FUNCTION
    
    
    FUNCTION GS_LDU_GET_S  (                                       _
                                   N_LDU      AS LONG              , _
                            BYVAL pU_LDU      AS UDT_LDU POINTER      )  AS STRING
    
       DIM i                                                 AS LONG
       DIM S_LDU                                             AS STRING
    
    
       IF       N_LDU = 00                                   OR _
                pU_LDU = 00                                  THEN
    
          FUNCTION = ""
    
       ELSE
    
          FOR i = 00 TO N_LDU - 01
    
             S_LDU   += TRIM$("x.kp1[" + STR$(i+1) + "]") + " " + STR$(@pU_LDU[i].xkp1) + $CRLF
    
          NEXT i
    
          FUNCTION = S_LDU
    
       END IF
    
    END FUNCTION
    
    
    FUNCTION GS_LDU_TEST   (                                          )  AS LONG
    
       '  Routine computes an Gauss Seidel test function
    
       DIM N_LDU                                             AS LONG
       DIM U_LDU()                                           AS UDT_LDU
       DIM U_LDU_PARM                                        AS UDT_LDU_PARM
    
    
       N_LDU = 3
       REDIM U_LDU(00 TO N_LDU - 01)
    
       U_LDU(00).L = %F
       U_LDU(01).L = 1.0#
       U_LDU(02).L = 1.0#
    
       U_LDU(00).D = 8.0#
       U_LDU(01).D = -7.0#
       U_LDU(02).D = 9.0#
    
       U_LDU(00).U = 1.0#
       U_LDU(01).U = 2.0#
       U_LDU(02).U = %F
    
       U_LDU(00).b = 8.0#
       U_LDU(01).b = -4.0#
       U_LDU(02).b = 12.0#
    
       U_LDU(00).xk = 0.0#
       U_LDU(01).xk = 0.0#
       U_LDU(02).xk = 0.0#
    
       U_LDU_PARM.N_ITR_MAX = 10
       U_LDU_PARM.e_ITR_MIN = 1.0E-6#
    
      'GS_LDU_xkp1    N_LDU, VARPTR(U_LDU(00))
       GS_LDU_SOLVE   N_LDU, VARPTR(U_LDU(00)), VARPTR(U_LDU_PARM)
    
       MSGBOX         GS_LDU_GET_S(N_LDU, VARPTR(U_LDU(00))), %MB_ICONINFORMATION OR %MB_TASKMODAL OR %MB_TOPMOST, $S_MSGBOX_HDR
    
    END FUNCTION
    
    
    FUNCTION CSP_FIND_INDEX(                                             _
                                   N_CSP      AS LONG                    , _
                            BYVAL pU_CSP      AS UDT_CSP        POINTER  , _
                            BYVAL pU_CSP_PARM AS UDT_CSP_PARM   POINTER  , _
                                   t_INTP     AS DOUBLE                     )  AS LONG
    
       '  Routine finds the ith index that an interpolated data point would lie between
       '  Note this ranges 00 to n - 01 - 01, thus returns -1 on fail which is tested
    
       DIM i                                                 AS LONG
    
    
       IF       N_CSP = 00                                   OR _
                pU_CSP = 00                                  OR _
                pU_CSP_PARM = 00                             THEN
    
          FUNCTION = -1
    
       ELSEIF  @pU_CSP_PARM.N_DATA = 00                      OR _     '  Test all allocations and pointers
               @pU_CSP_PARM.pt_DATA = 00                     OR _
               @pU_CSP_PARM.py_DATA = 00                     OR _
               @pU_CSP_PARM.N_INTP = 00                      OR _
               @pU_CSP_PARM.pt_INTP = 00                     OR _
               @pU_CSP_PARM.py_INTP = 00                     OR _
               @pU_CSP_PARM.N_LDU = 00                       OR _
               @pU_CSP_PARM.pU_LDU = 00                      THEN
    
    
       ELSEIF  NOT @pU_CSP_PARM.N_DATA = @pU_CSP_PARM.N_LDU  THEN     '  LDU allocation must match the input data size
    
          FUNCTION = -1
    
       ELSE
    
          '  Trap the boundaries
    
          IF       t_INTP <= @[email protected]_DATA[00]                THEN
    
             FUNCTION = 00
    
          ELSEIF   t_INTP >= @[email protected]_DATA[N_CSP - 01 - 01]   THEN
    
             FUNCTION = N_CSP - 01 - 01
    
          ELSE
    
             FOR i = 00 TO N_CSP - 01 - 01
    
                IF t_INTP >= @[email protected]_DATA[i]                 AND _
                   t_INTP <  @[email protected]_DATA[i + 01]            THEN
    
                   FUNCTION = i
                   EXIT FUNCTION
    
                END IF
    
             NEXT i
    
             '  Fail condition is encountered
             FUNCTION = -1
    
          END IF
    
       END IF
    
    END FUNCTION
    
    
    FUNCTION CSP_SOLVE     (                                             _
                                   N_CSP      AS LONG                    , _
                            BYVAL pU_CSP      AS UDT_CSP        POINTER  , _
                            BYVAL pU_CSP_PARM AS UDT_CSP_PARM   POINTER     )  AS LONG
    
       '  Routine computes a cubic spline solution to a set of input data by pointer
    
       DIM i                                                 AS LONG
       DIM i_INTP                                            AS LONG
    
    
       IF       N_CSP = 00                                   OR _
                pU_CSP = 00                                  OR _
                pU_CSP_PARM = 00                             THEN
    
          FUNCTION = %F
    
       ELSEIF  @pU_CSP_PARM.N_DATA = 00                      OR _     '  Test all allocations and pointers
               @pU_CSP_PARM.pt_DATA = 00                     OR _
               @pU_CSP_PARM.py_DATA = 00                     OR _
               @pU_CSP_PARM.N_INTP = 00                      OR _
               @pU_CSP_PARM.pt_INTP = 00                     OR _
               @pU_CSP_PARM.py_INTP = 00                     OR _
               @pU_CSP_PARM.N_LDU = 00                       OR _
               @pU_CSP_PARM.pU_LDU = 00                      THEN
    
    
       ELSEIF  NOT @pU_CSP_PARM.N_DATA = @pU_CSP_PARM.N_LDU  THEN     '  LDU allocation must match the input data size
    
          FUNCTION = %F
    
       ELSE
    
          '  Compute b and h coefficients
    
          FOR i = 00 TO N_CSP - 01 - 01
    
             @pU_CSP[i].h = @[email protected]_DATA[i + 01] - @[email protected]_DATA[i]
    
             IF NOT @pU_CSP[i].h = 0.0#                      THEN
    
                @pU_CSP[i].b = 1.0# / @pU_CSP[i].h * (@[email protected]_DATA[i + 01] - @[email protected]_DATA[i])
    
             END IF
    
          NEXT i
    
    
          '  Compute v and u coefficients
    
          FOR i = 01 TO N_CSP - 01 - 01
    
             @pU_CSP[i].v = 2.0# * (@pU_CSP[i - 01].h + @pU_CSP[i].h)
             @pU_CSP[i].u = 6.0# * (@pU_CSP[i].b - @pU_CSP[i - 01].b)
    
          NEXT i
    
    
          '  Load the LDU coefficients to prepare for the tridiagonal solution
    
          FOR i = 00 TO @pU_CSP_PARM.N_LDU - 01
    
             '  Load the LDU vectors
             '  Note the IIF operator for the boundary condition (i = 00 OR i = N_CSP - 01) which processes S''(x) = 0 by setting L, U = 0 and D = 1
             @[email protected]_LDU[i].L = IIF(i = 00 OR i = N_CSP - 01, 0.0#, IIF(i = 01             , %F, @pU_CSP[i - 01].h))
             @[email protected]_LDU[i].D = IIF(i = 00 OR i = N_CSP - 01, 1.0#, @pU_CSP[i].v                                   )
             @[email protected]_LDU[i].U = IIF(i = 00 OR i = N_CSP - 01, 0.0#, IIF(i = N_CSP - 01 - 01, %F, @pU_CSP[i].h)     )
    
             '  Load the b solution source vector which in this case is the [u] vector
             @[email protected]_LDU[i].b = @pU_CSP[i].u
    
          NEXT i
    
    
          '  Solve the cubic spline [z] vector for S.i''(x) using the tridiagonal Gauss Seidel LDU solver
    
          GS_LDU_SOLVE   @pU_CSP_PARM.N_LDU, @pU_CSP_PARM.pU_LDU, VARPTR(@pU_CSP_PARM.U_LDU_PARM)
          #IF %F
          MSGBOX GS_LDU_GET_S(@pU_CSP_PARM.N_LDU, @pU_CSP_PARM.pU_LDU), %MB_ICONINFORMATION OR %MB_TASKMODAL OR %MB_TOPMOST, $S_MSGBOX_HDR
          #ENDIF
    
          '  Range over the span of the interpolated t points
    
          FOR i_INTP = 00 TO @pU_CSP_PARM.N_INTP - 01
    
             '  Create the uniform interpolation points
             @[email protected]_INTP[i_INTP] = @[email protected]_DATA[00] + (@[email protected]_DATA[@pU_CSP_PARM.N_DATA - 01] - @[email protected]_DATA[00]) * i_INTP / (@pU_CSP_PARM.N_INTP - 01)
    
             '  Find the data index that each point lies in
             CSP_FIND_INDEX    N_CSP                            , _
                               pU_CSP                           , _
                               pU_CSP_PARM                      , _
                               @[email protected]_INTP[i_INTP]    TO i
    
             '  Error check the index for -1 from the find index function
             IF NOT i = -1                                   THEN
    
                '  Compute the cubic spline interpolation y = ao + a1 t + a2 t^2 + a3 t^3
                '  Note the LDU solved [x]kp1 term replaces the [z] solution vector by variable name
                @[email protected]_INTP[i_INTP] = @[email protected]_LDU[i + 01].xkp1 / (6.0# * @pU_CSP[i].h) * (@[email protected]_INTP[i_INTP] - @[email protected]_DATA[i]) ^ 3                                                  + _
                                                @[email protected]_LDU[i].xkp1      / (6.0# * @pU_CSP[i].h) * (@[email protected]_DATA[i + 01] - @[email protected]_INTP[i_INTP]) ^ 3                                             + _
                                                (@[email protected]_DATA[i + 01] / @pU_CSP[i].h - @[email protected]_LDU[i + 01].xkp1 / 6.0# * @pU_CSP[i].h) * (@[email protected]_INTP[i_INTP] - @[email protected]_DATA[i])       + _
                                                (@[email protected]_DATA[i]      / @pU_CSP[i].h - @[email protected]_LDU[i].xkp1      / 6.0# * @pU_CSP[i].h) * (@[email protected]_DATA[i + 01] - @[email protected]_INTP[i_INTP])
    
             END IF
    
          NEXT i_INTP
    
          FUNCTION = %T
    
       END IF
    
    END FUNCTION
    
    
    FUNCTION CSP_SET_CLIP  (                                             _
                                   N_CSP      AS LONG                    , _
                            BYVAL pU_CSP      AS UDT_CSP        POINTER  , _
                            BYVAL pU_CSP_PARM AS UDT_CSP_PARM   POINTER     )  AS LONG
    
       '  Routine sets cubic spline solution to the clipboard for easy plotting
    
       DIM i                                                 AS LONG
       DIM i_INTP                                            AS LONG
    
       DIM S_CLIP                                            AS STRING
    
    
       IF       N_CSP = 00                                   OR _
                pU_CSP = 00                                  OR _
                pU_CSP_PARM = 00                             THEN
    
          FUNCTION = %F
    
       ELSEIF  @pU_CSP_PARM.N_DATA = 00                      OR _     '  Test all allocations and pointers
               @pU_CSP_PARM.pt_DATA = 00                     OR _
               @pU_CSP_PARM.py_DATA = 00                     OR _
               @pU_CSP_PARM.N_INTP = 00                      OR _
               @pU_CSP_PARM.pt_INTP = 00                     OR _
               @pU_CSP_PARM.py_INTP = 00                     OR _
               @pU_CSP_PARM.N_LDU = 00                       OR _
               @pU_CSP_PARM.pU_LDU = 00                      THEN
    
    
       ELSEIF  NOT @pU_CSP_PARM.N_DATA = @pU_CSP_PARM.N_LDU  THEN     '  LDU allocation must match the input data size
    
          FUNCTION = %F
    
       ELSE
    
          '  Range over the span of the interpolated t points
    
          FOR i_INTP = 00 TO @pU_CSP_PARM.N_INTP - 01
    
             S_CLIP  += STR$(@[email protected]_INTP[i_INTP]) + " " + STR$(@[email protected]_INTP[i_INTP]) + $CRLF
    
          NEXT i_INTP
    
          CLIPBOARD SET TEXT S_CLIP
    
          FUNCTION = %T
    
       END IF
    
    END FUNCTION
    
    
    FUNCTION CSP_GET_S     (                                             _
                                   N_CSP      AS LONG                    , _
                            BYVAL pU_CSP      AS UDT_CSP        POINTER  , _
                            BYVAL pU_CSP_PARM AS UDT_CSP_PARM   POINTER     )  AS STRING
    
       DIM i                                                 AS LONG
       DIM i_INTP                                            AS LONG
       DIM S_CSP                                             AS STRING
    
    
       IF       N_CSP = 00                                   OR _
                pU_CSP = 00                                  OR _
                pU_CSP_PARM = 00                             THEN
    
          FUNCTION = ""
    
       ELSEIF  @pU_CSP_PARM.N_DATA = 00                      OR _     '  Test all allocations and pointers
               @pU_CSP_PARM.pt_DATA = 00                     OR _
               @pU_CSP_PARM.py_DATA = 00                     OR _
               @pU_CSP_PARM.N_INTP = 00                      OR _
               @pU_CSP_PARM.pt_INTP = 00                     OR _
               @pU_CSP_PARM.py_INTP = 00                     OR _
               @pU_CSP_PARM.N_LDU = 00                       OR _
               @pU_CSP_PARM.pU_LDU = 00                      THEN
    
    
       ELSEIF  NOT @pU_CSP_PARM.N_DATA = @pU_CSP_PARM.N_LDU  THEN     '  LDU allocation must match the input data size
    
          FUNCTION = ""
    
       ELSE
    
          #IF 00
          FOR i = 00 TO N_CSP - 01
    
             S_CSP   += TRIM$("b[" + STR$(i+1) + "]") + " " + STR$(@pU_CSP[i].b) + " " + _
                        TRIM$("v[" + STR$(i+1) + "]") + " " + STR$(@pU_CSP[i].v) + " " + _
                        TRIM$("u[" + STR$(i+1) + "]") + " " + STR$(@pU_CSP[i].u) + $CRLF
    
          NEXT i
          #ENDIF
    
          #IF 00
          FOR i = 00 TO @pU_CSP_PARM.N_LDU - 01
    
             S_CSP   += TRIM$("L[" + STR$(i+1) + "]") + " " + STR$(@[email protected]_LDU[i].L) + " " + _
                        TRIM$("D[" + STR$(i+1) + "]") + " " + STR$(@[email protected]_LDU[i].D) + " " + _
                        TRIM$("U[" + STR$(i+1) + "]") + " " + STR$(@[email protected]_LDU[i].U) + $CRLF
    
          NEXT i
          #ENDIF
    
          #IF -1
          FOR i_INTP = 00 TO @pU_CSP_PARM.N_INTP - 01
    
             S_CSP   += "[i INTP, t_INTP, y_INTP] " + STR$(i_INTP) + " " + STR$(@[email protected]_INTP[i_INTP]) + " " + STR$(@[email protected]_INTP[i_INTP]) + $CRLF
    
          NEXT i_INTP
          #ENDIF
    
          FUNCTION = S_CSP
    
       END IF
    
    END FUNCTION
    
    
    FUNCTION CSP_TEST      (                                          )  AS LONG
    
       '  Routine computes an cubic spline test function
    
       DIM N_CSP                                             AS LONG
       DIM U_CSP()                                           AS UDT_CSP
       DIM U_CSP_PARM                                        AS UDT_CSP_PARM
    
       DIM N_DATA                                            AS LONG
       DIM t_DATA()                                          AS DOUBLE
       DIM y_DATA()                                          AS DOUBLE
    
       DIM N_INTP                                            AS LONG
       DIM t_INTP()                                          AS DOUBLE
       DIM y_INTP()                                          AS DOUBLE
    
       DIM N_LDU                                             AS LONG
       DIM U_LDU()                                           AS UDT_LDU
    
    
       '  Create a test input data set
       '  Note the assumption is that the load data is sorted
       '     0.9   0
       '     1.3   1.5
       '     1.9   1.85
       '     2.1   2.1
       N_DATA = 14
       REDIM t_DATA(00 TO N_DATA - 01) : ARRAY ASSIGN t_DATA() = 0.00, 1.00, 3.00, 4.00, 5.00, 6.00, 8.00, 11.00, 12.00, 13.00, 14.00, 16.00, 17.00, 19.00
       REDIM y_DATA(00 TO N_DATA - 01) : ARRAY ASSIGN y_DATA() = 0.00, 2.60, 23.16, 27.57, 24.26, 16.63, 30.41, 47.20, 50.03, 60.33, 59.89, 71.18, 84.27, 77.69
       U_CSP_PARM.N_DATA    = N_DATA
       U_CSP_PARM.pt_DATA   = VARPTR(t_DATA(00))
       U_CSP_PARM.py_DATA   = VARPTR(y_DATA(00))
    
       '  Allocate the LDU solution data structure
       N_LDU = N_DATA
       REDIM U_LDU(00 TO N_LDU - 01)
       U_CSP_PARM.N_LDU     = N_DATA
       U_CSP_PARM.pU_LDU    = VARPTR(U_LDU(00))
    
       '  Allocate the cubic spline coefficients to the number of data points
       N_CSP = N_DATA
       REDIM U_CSP(00 TO N_CSP - 01)
    
       '  Set the LDU solve parameters
       U_CSP_PARM.U_LDU_PARM.N_ITR_MAX = 5
       U_CSP_PARM.U_LDU_PARM.e_ITR_MIN = 1.0E-3#
    
       '  Create a test input interpolation set
       '  Note this is larger than the input data set to create smoothness
       N_INTP = 100
       REDIM t_INTP(00 TO N_INTP - 01)
       REDIM y_INTP(00 TO N_INTP - 01)
       U_CSP_PARM.N_INTP    = N_INTP
       U_CSP_PARM.pt_INTP   = VARPTR(t_INTP(00))
       U_CSP_PARM.py_INTP   = VARPTR(y_INTP(00))
    
    
       '  Solve for cubic spline
       CSP_SOLVE      N_CSP, VARPTR(U_CSP(00)), VARPTR(U_CSP_PARM)
       CSP_SET_CLIP   N_CSP, VARPTR(U_CSP(00)), VARPTR(U_CSP_PARM)
    
       #IF %T
       MSGBOX      CSP_GET_S(N_CSP, VARPTR(U_CSP(00)), VARPTR(U_CSP_PARM)), %MB_ICONINFORMATION OR %MB_TASKMODAL OR %MB_TOPMOST, $S_MSGBOX_HDR
       #ENDIF
    
    END FUNCTION

    Comment


    • #3
      Can you show how to plot out the interpolated data graphically?

      Comment

      Working...
      X