Announcement

Collapse
No announcement yet.

Cubic Spline

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

  • Cubic Spline

    Here is a cubic spline function I wrote. This turned out to be more complex than I thought and required a fast tridiagonal matrix solution. The following produces a 10 point interpolation of fixed input data. This was done for verification purposes so I could trace through the code and debug properly.

    Can't seem to understand how to retain formatting

  • #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]U_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)
    
    ' When interpolated to 10 interpolation points, this data should produce the following spline result
    ' 0.90000 0.00000
    ' 1.03333 0.59125
    ' 1.16667 1.11407
    ' 1.30000 1.50000
    ' 1.43333 1.70412
    ' 1.56667 1.77553
    ' 1.70000 1.78685
    ' 1.83333 1.81069
    ' 1.96667 1.91647
    ' 2.10000 2.10000
    #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
      Re: Formatting code -
      Posted Source Code Losing Indentation/Formatting - PowerBASIC Peer Support Community
      (in FAQ section)

      Cheers,
      Dale

      Comment


      • #4
        Sample output
        Click image for larger version

Name:	Cubic Spline.png
Views:	269
Size:	21.4 KB
ID:	812088

        Comment


        • #5
          OK, went and got the formatting worked out in the source code, posted at https://forum.powerbasic.com/forum/u...6-cubic-spline

          Some notes
          • There is a test code run from PBMAIN
          • Currently creates a 100 point spline for the 14 or so data points
          • Seems to run very fast given the numerical tridiagonal GS solution. The code copies the data to the clipboard.
          • Taking the interpolation our of the data structure is done as follows:
          Code:
          ' 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

          Comment


          • #6
            Apologies, I did not cut and paste the entire code. The correct version is in a second source posting.

            Comment


            • #7
              How to plot these interpolated data graphically?

              Comment


              • #8
                Ok, how to write the program to read in the x,y coordinates and plot them out in graphics?
                The x,y data points coordinates are :

                0 0
                .191919191919192 .179817579631133
                .383838383838384 .432868973254548
                .575757575757576 .832387994862525
                .767676767676768 1.45160845844734
                .95959595959596 2.36376417800129
                1.15151515151515 3.63367096244116
                1.34343434343434 5.2342282119494
                1.53535353535354 7.0996289267785
                1.72727272727273 9.16463423210082
                1.91919191919192 11.3640052530887
                2.11111111111111 13.6325031149145
                2.3030303030303 15.9048889427507
                2.49494949494949 18.1159238617695
                2.68686868686869 20.2003689971433
                2.87878787878788 22.0929854740445
                3.07070707070707 23.7286447848093
                3.26262626262626 25.0702315034977
                3.45454545454545 26.1211543371911
                3.64646464646465 26.8869723557255
                3.83838383838384 27.3732446289368
                4.03030303030303 27.5854380957347
                4.22222222222222 27.5243421419463
                4.41414141414141 27.1804484003616
                4.60606060606061 26.5421043141059
                4.7979797979798 25.5976573263042
                4.98989898989899 24.3354548800817
                5.18181818181818 22.7742770454954
                5.37373737373737 21.0580475087966
                5.56565656565656 19.3711461378896
                5.75757575757576 17.8981037435186
                5.94949494949495 16.8234511364277
                6.14141414141414 16.3105151448164
                6.33333333333333 16.3735607988952
                6.52525252525252 16.9425408429691
                6.71717171717172 17.9451706707874
                6.90909090909091 19.3091656760993
                7.1010101010101 20.962241252654
                7.29292929292929 22.8321127942007
                7.48484848484848 24.8464956944887
                7.67676767676768 26.9331053472672
                7.86868686868687 29.0196571462853
                8.06060606060606 31.0355057620906
                8.25252525252525 32.9363718561446
                8.44444444444444 34.7162803649456
                8.63636363636364 36.37587178856
                8.82828282828283 37.9157866270542
                9.02020202020202 39.3366653804949
                9.21212121212121 40.6391485489484
                9.4040404040404 41.8238766324812
                9.5959595959596 42.8914901311598
                9.78787878787879 43.8426295450507
                9.97979797979798 44.6779353742205
                10.1717171717172 45.3980481187354
                10.3636363636364 46.0036082786622
                10.5555555555556 46.4952563540672
                10.7474747474747 46.8736328450168
                10.9393939393939 47.1393782515777
                11.1313131313131 47.3019592991724
                11.3232323232323 47.453581992848
                11.5151515151515 47.7385009249251
                11.7070707070707 48.3022350526334
                11.8989898989899 49.2903033332027
                12.0909090909091 50.8397710534004
                12.2828282828283 52.8909149148208
                12.4747474747475 55.1734758192623
                12.6666666666667 57.4058296650598
                12.8585858585858 59.3063523505481
                13.050505050505 60.5951951292494
                13.2424242424242 61.1512496440995
                13.4343434343434 61.149031452269
                13.6262626262626 60.7965074693367
                13.8181818181818 60.3016446108815
                14.010101010101 59.8723595064829
                14.2020202020202 59.6760479429241
                14.3939393939394 59.7456956593398
                14.5858585858586 60.0844471096318
                14.7777777777778 60.6954467477018
                14.969696969697 61.5818390274516
                15.1616161616162 62.7467684027832
                15.3535353535354 64.1933793275981
                15.5454545454545 65.9248162557983
                15.7373737373737 67.9442236412855
                15.9292929292929 70.2547459379616
                16.1212121212121 72.8518016686624
                16.3131313131313 75.6278526904268
                16.5050505050505 78.4024799439694
                16.6969696969697 80.9937482450886
                16.8888888888889 83.2197224095828
                17.0808080808081 84.9015097753544
                17.2727272727273 85.9650123625394
                17.4646464646465 86.4649402169096
                17.6565656565656 86.4639127532224
                17.8484848484848 86.0245493862349
                18.040404040404 85.2094695307043
                18.2323232323232 84.0812926013878
                18.4242424242424 82.7026380130425
                18.6161616161616 81.1361251804258
                18.8080808080808 79.4443735182946
                19 77.6900024414062

                Comment


                • #9
                  Modifying the PB's sample prog called "Sine.bas" and adding to read in the data file in post #8
                  but still not been able to plot the graphical plot ?


                  '====================================================================
                  '

                  ' Scale a GRAPHIC control to pixel-based coordinates
                  ' and draw a plot of given points.
                  '
                  '====================================================================

                  #COMPILER PBWIN 10
                  #COMPILE EXE
                  #RESOURCE MANIFEST, 1, "XPTheme.xml"
                  #DIM ALL


                  %IDC_GRAPHIC1 = 110


                  ' Draw a plot of points
                  '
                  SUB DrawXYplot (BYVAL hDlg AS DWORD, BYVAL ID AS LONG)

                  LOCAL x, y, mp, wh , hzRadDeg, Pi AS DOUBLE ' AS LONG

                  ' read in the data file
                  DIM xcoord(1 TO 1000) AS LOCAL DOUBLE
                  DIM ycoord(1 TO 1000) AS LOCAL DOUBLE
                  LOCAL zf , jc , jj AS LONG
                  LOCAL chkst AS STRING
                  zf = FREEFILE
                  jc = 0
                  OPEN "Interpolated points.txt" FOR INPUT AS #zf
                  WHILE NOT EOF(zf)
                  INCR jc
                  INPUT #zf, xcoord(jc) , ycoord(jc)
                  WEND

                  CLOSE #zf

                  ? "number of data points" + STR$( jc-1 )
                  FOR jj = 1 TO jc-1
                  chkst = chkst+ STR$(xcoord(jj)) + " " +STR$( ycoord(jj)) + $CRLF
                  NEXT jj
                  ? chkst


                  GRAPHIC SCALE PIXELS ' Convert to pixels
                  GRAPHIC GET CANVAS TO x, y ' Get writable size
                  mp = y / 2 ' Midpoint
                  wh = mp - 2 ' Wave height
                  Pi = 4 * ATN(1) ' Calculate Pi
                  hzRadDeg = Pi / 180 ' Degrees to Radians

                  GRAPHIC ATTACH hDlg, ID, REDRAW ' Use faster, buffered draw
                  GRAPHIC COLOR %RGB_BLUE, %RGB_WHITE ' Blue line, white background
                  GRAPHIC CLEAR ' Paint the white background
                  GRAPHIC LINE (0, mp) - (x, mp), %GRAY ' Draw midline

                  ' y = SIN(hzRadDeg) * mp ' Calculate start point
                  ' GRAPHIC SET POS (0, mp) ' Set start position

                  ' FOR x = 0 TO 360 ' Run a full circle loop
                  ' y = SIN(hzRadDeg * x) * wh ' y-pos for each degree
                  ' GRAPHIC LINE STEP - (x, mp + y) ' Draw line
                  ' NEXT
                  ' GRAPHIC REDRAW ' Show it all

                  x = xcoord(1) ' the start point
                  y = ycoord(1)
                  GRAPHIC SET POS ( x,y)
                  FOR jj = 2 TO jc-1
                  x = xcoord(jj)
                  y = ycoord(jj)
                  GRAPHIC LINE STEP - (x ,y) , %RED
                  NEXT jj
                  GRAPHIC REDRAW

                  ' Finally, tell what is drawn in dialog caption.
                  DIALOG SET TEXT hDlg, "Plot of points"

                  END SUB



                  ' Main dialog callback procedure
                  '
                  CALLBACK FUNCTION DlgProc () AS LONG

                  SELECT CASE CB.MSG

                  CASE %WM_INITDIALOG ' <- Sent right before the dialog is shown

                  CASE %WM_COMMAND ' <- A control is calling
                  SELECT CASE CB.CTL ' <- Look at control's id

                  CASE %IDCANCEL
                  IF CB.CTLMSG = %BN_CLICKED THEN ' Exit on Esc
                  DIALOG END CB.HNDL
                  END IF

                  END SELECT

                  END SELECT

                  END FUNCTION



                  ' Program entry point
                  '
                  FUNCTION PBMAIN () AS LONG

                  LOCAL hDlg AS DWORD, w, h AS LONG

                  DIALOG NEW 0, "Plot of interpolated points",,, 200, 200, %WS_CAPTION OR %WS_SYSMENU, _
                  0 TO hDlg

                  CONTROL ADD GRAPHIC, hDlg, %IDC_GRAPHIC1, "", 0, 0, 200, 200

                  ' Set up a pixel-based coordinate system in the Graphic control
                  CONTROL GET CLIENT hDlg, %IDC_GRAPHIC1 TO w, h ' Get client size
                  DIALOG UNITS hDlg, w, h TO PIXELS w, h ' Convert to pixels
                  GRAPHIC ATTACH hDlg, %IDC_GRAPHIC1
                  GRAPHIC SCALE (0, 0) - (w, h) ' Scale to pixel coordinate system

                  DrawXYplot hDlg, %IDC_GRAPHIC1

                  DIALOG SHOW MODAL hDlg, CALL DlgProc

                  END FUNCTION






                  Comment


                  • #10
                    any suggestion to plot these points more accurately?

                    Comment


                    • #11
                      Originally posted by Mannish Bhandari View Post
                      Ok, how to write the program to read in the x,y coordinates and plot them out in graphics?
                      THat's a bit off topic for this thread.
                      There are plenty of examples of graphic displays of data if you search the forum.

                      The last example I posted is here: https://forum.powerbasic.com/forum/u...728#post795728

                      Here's another simple example. I'll leave it to you to draw in the axis labels etc.

                      '
                      Code:
                      #COMPILE EXE
                      #DIM ALL
                      FUNCTION PBMAIN() AS LONG
                          LOCAL x AS LONG 'loopcount
                          LOCAL lngDatapoint AS LONG
                          LOCAL strFile,strT AS STRING
                          LOCAL minY,maxY AS EXT
                          minY =  1.2*10^4932
                          maxY =  -1.2*10^4932
                          DIM arrX(1 TO 100) AS EXT ' Use separate arrays for ease of sort
                          DIM arrY(1 TO 100) AS EXT ' If initial data is not sorted
                      
                          DISPLAY OPENFILE , , , "Select data file", EXE.PATH$, CHR$("All Files", 0, "*.*", 0),"Mannishdata.tab","" , 0 TO strFile
                          IF NOT ISFILE(strFIle) THEN
                              ? "File Not Found - Exiting"
                              EXIT FUNCTION
                          END IF
                          OPEN strFile FOR INPUT AS #1
                          WHILE NOT EOF(#1)
                              INCR lngDatapoint
                              IF UBOUND(arrX()) < lngDatapoint THEN REDIM PRESERVE arrX(lngDataPoint + 100):REDIM PRESERVE arrY(lngDataPoint + 100)
                              LINE INPUT #1, strT
                              arrx(lngDatapoint) = VAL(PARSE$(strT,$TAB,1))
                              arrY(lngDatapoint) = VAL(PARSE$(strT,$TAB,2))
                              IF arrY(lngDatapoint) > MaxY THEN MaxY = arrY(lngDatapoint)
                              IF arrY(lngDatapoint) < MinY THEN MinY = arrY(lngDatapoint)
                      
                          WEND
                          CLOSE #1
                      
                          REDIM PRESERVE arrX(1 TO lngDatapoint)
                          REDIM PRESERVE arrY(1 TO lngDatapoint)
                          ARRAY SORT ArrX(), TAGARRAY ArrY()
                          Plotit ArrX(),ArrY(), MinY,MaxY, strFile
                      END FUNCTION
                      
                      FUNCTION Plotit(ArrX() AS EXT, arrY() AS EXT ,minY AS EXT,MaxY AS EXT,strF AS STRING) AS LONG
                          LOCAL hWin AS DWORD
                          LOCAL  minX,maxX, x AS EXT
                          LOCAL SMinX,SMinY,SMaxX,SMaxY AS EXT
                          'Get graph bounds
                          minX = ArrX(LBOUND(arrx()))
                          maxX = ArrX(UBOUND(arrx()))
                          'put 10% margins round data
                          SMinX = minx - (maxX-minX) *.1
                          sMinY = miny - (maxy-miny) *.1
                          sMaxX = maxx + (maxX-minX) *.1
                          sMaxY = maxy + (maxy-miny) *.1
                      
                          GRAPHIC WINDOW NEW  strF , 50, 50, 800, 400  TO hWin
                          GRAPHIC SCALE (Sminx,Smaxy) - (SmaxX,SMinY)
                          'Draw the axes
                          GRAPHIC COLOR %BLACK,-1
                          GRAPHIC WIDTH 2
                          GRAPHIC LINE (MinX,MinY)-(Minx,MAxY)
                          GRAPHIC LINE (Minx,MinY)-(MaxX,MinY)
                      
                          'Plot the data
                          GRAPHIC COLOR %BLUE,-1
                      
                          GRAPHIC SET POS (arrx(1),arry(1))
                          FOR x = 2 TO UBOUND(arrx())
                              GRAPHIC LINE - (arrX(x),arrY(x))
                          NEXT
                          GRAPHIC WAITKEY$
                          GRAPHIC DETACH
                          GRAPHIC WINDOW END hWin
                      END FUNCTION
                      '
                      Click image for larger version  Name:	MannishGraph.jpg Views:	0 Size:	31.6 KB ID:	812324

                      Comment


                      • #12
                        Thanks Stuart, that's what I need

                        Comment


                        • #13
                          Improved version of Stuart's prog


                          ' prog to read in a data file and plots out its curve


                          #COMPILE EXE
                          #DIM ALL



                          FUNCTION PBMAIN() AS LONG
                          LOCAL lngDatapoint AS LONG
                          LOCAL strFile AS STRING
                          LOCAL minY,maxY AS EXT
                          LOCAL ff AS LONG

                          'Inititally make these min and max extremely small and large
                          ' so that we can determine the real min and max later on
                          minY = 1.2*10^4932
                          maxY = -1.2*10^4932



                          ' Use separate arrays for ease of sort
                          ' If initial data is not sorted
                          DIM arrX() AS EXT
                          DIM arrY() AS EXT


                          ' use open file dialog to allow data file selection
                          DISPLAY OPENFILE , , , "Select data file", _
                          EXE.PATH$, CHR$("All Files", 0, "*.*", 0),_
                          "Interpolated points.txt","" , 0 TO strFile
                          IF NOT ISFILE(strFIle) THEN
                          ? "File Not Found - Exiting"
                          EXIT FUNCTION
                          END IF

                          ' Now open this selected data file
                          ff = FREEFILE
                          OPEN strFile FOR INPUT AS #ff
                          WHILE NOT EOF(#ff)
                          INCR lngDatapoint
                          REDIM PRESERVE arrX(1 TO lngDatapoint)
                          REDIM PRESERVE arrY(1 TO lngDatapoint)
                          INPUT #ff, arrX(lngDatapoint) , arrY(lngDatapoint)

                          IF arrY(lngDatapoint) > MaxY THEN
                          MaxY = arrY(lngDatapoint)
                          END IF
                          IF arrY(lngDatapoint) < MinY THEN
                          MinY = arrY(lngDatapoint)
                          END IF
                          WEND
                          CLOSE #ff

                          ' sort x coordinates array in its ascending order
                          ' tagging along its y coordinates array
                          ARRAY SORT ArrX(), TAGARRAY ArrY()

                          Plotit ArrX(),ArrY(), MinY,MaxY, PATHNAME$(NAMEX, strFile )
                          END FUNCTION



                          ' Plot the graph with straight lines
                          ' if there are sufficient many points then the graph can smoothen
                          ' into curves
                          FUNCTION Plotit(ArrX() AS EXT, arrY() AS EXT ,minY AS EXT,_
                          MaxY AS EXT, strF AS STRING) AS LONG
                          LOCAL hWin AS DWORD
                          LOCAL minX,maxX AS EXT
                          LOCAL xc , yc AS EXT
                          LOCAL jj , deltaY AS LONG
                          LOCAL SMinX,SMinY,SMaxX,SMaxY AS EXT


                          LOCAL hFont1, hFont2, hFont3 AS LONG
                          FONT NEW "Lucida Console", 10, 0, 0, 0, 0 TO hFont1 ' set up fonts to be used later
                          FONT NEW "Lucida Console", 10, 0, 0, 0, 900 TO hFont2 ' right angle font (needs to be true type)
                          FONT NEW "Lucida Console", 2, 0, 0, 0, 0 TO hFont3 ' small font

                          'Get graph bounds
                          minX = ArrX(LBOUND(arrx()))
                          maxX = ArrX(UBOUND(arrx()))
                          'put 10% margins round data
                          SMinX = minx - (maxX-minX) *.1
                          sMinY = miny - (maxy-miny) *.1
                          sMaxX = maxx + (maxX-minX) *.1
                          sMaxY = maxy + (maxy-miny) *.1

                          GRAPHIC WINDOW NEW strF , 50, 50, 800, 500 TO hWin
                          GRAPHIC SCALE (Sminx,Smaxy) - (SmaxX,SMinY)

                          ' move the axes upwards with deltaY
                          ' so that can plot the graph better
                          deltaY = 2
                          'Draw the axes
                          GRAPHIC COLOR %BLACK,-1
                          GRAPHIC WIDTH 2
                          GRAPHIC LINE (MinX,MinY+deltaY)-(Minx,MAxY+deltaY)
                          GRAPHIC LINE (Minx,MinY+deltaY)-(MaxX,MinY+deltaY)

                          ' display the max Y and Min Y
                          GRAPHIC SET FONT hFont1
                          LOCAL MaxYSt, MinYSt AS STRING
                          MaxYSt = TRIM$(LEFT$(STR$(MaxY)+ " ",5))
                          GPRINTX(MaxYst , Minx-1,MAxY,0)
                          MinYSt = TRIM$(LEFT$(STR$(MinY)+ " ",5))
                          GPRINTX(MinYSt ,MinX,MinY,0)

                          ' display the max X
                          LOCAL MaxXSt AS STRING
                          MaxXSt = TRIM$(LEFT$(STR$(MaxX)+ " ",5))
                          GPRINTX(MaxXst ,MaxX,MinY+deltaY,0)

                          ' display the X , Y axes labels
                          LOCAL xLab , yLab AS STRING
                          xLab = "X axis"
                          yLab = "Y axis"
                          GPRINTX(xLab ,MaxX - (MaxX/2),MinY - 2, 1)
                          GRAPHIC SET FONT hFont2
                          GPRINTX(yLab , Minx-1,MAxY - (MaxY/2)+deltaY, 1)


                          'Plot the data curve
                          ' sets the starting point first for the first line
                          GRAPHIC SET POS (arrx(1),arry(1)+deltaY)
                          FOR jj = 2 TO UBOUND(arrx())
                          GRAPHIC COLOR %CYAN,-1
                          GRAPHIC LINE - (arrX(jj),arrY(jj)+deltaY)
                          NEXT jj


                          ' draws the points with smaller font
                          GRAPHIC SET FONT hFont3
                          FOR jj = 2 TO UBOUND(arrx())
                          ' unfortunately the letter x is too large
                          ' GRAPHIC COLOR %RED,-1
                          ' xc = arrX(jj) - GRAPHIC(CELL.SIZE.X) / 2 ' adjust point positions to suit "X" center
                          ' yc = arrY(jj)+deltaY + GRAPHIC(CELL.SIZE.Y) / 2
                          ' GRAPHIC SET POS (xc, yc)
                          ' GRAPHIC PRINT "x"
                          GRAPHIC SET PIXEL (arrX(jj),arrY(jj)+deltaY), %RGB_DARKRED
                          NEXT jj

                          GRAPHIC WAITKEY$
                          GRAPHIC DETACH
                          GRAPHIC WINDOW END hWin
                          END FUNCTION


                          ' display the text
                          SUB GPRINTX(TX AS STRING,XS AS EXT,YS AS EXT, BGcol AS LONG)
                          IF BGCol = 1 THEN
                          'yellow color background
                          GRAPHIC COLOR %RGB_DARKGREEN, %RGB_LIGHTYELLOW
                          ELSE
                          ' default background
                          GRAPHIC COLOR %BLUE, -1
                          END IF
                          GRAPHIC SET POS(XS,YS)
                          GRAPHIC PRINT TX
                          END SUB

                          Comment


                          • #14
                            Click image for larger version  Name:	graph plot.png Views:	0 Size:	10.2 KB ID:	812353


                            only problem with the small size of the points drawn using the following statement, how to use ellipse to draw these points?

                            GRAPHIC SET PIXEL (arrX(jj),arrY(jj)+deltaY), %RGB_DARKRED

                            Comment


                            • #15
                              Improved version with marks on axes


                              ' prog to read in a data file and plots out its curve
                              ' adapted from Stuart McLachlan

                              #COMPILE EXE
                              #DIM ALL



                              FUNCTION PBMAIN() AS LONG
                              LOCAL lngDatapoint AS LONG
                              LOCAL strFile AS STRING
                              LOCAL minY,maxY AS EXT
                              LOCAL ff AS LONG

                              'Inititally make these min and max extremely small and large
                              ' so that we can determine the real min and max later on
                              minY = 1.2*10^4932
                              maxY = -1.2*10^4932



                              ' Use separate arrays for ease of sort
                              ' If initial data is not sorted
                              DIM arrX() AS EXT
                              DIM arrY() AS EXT


                              ' use open file dialog to allow data file selection
                              DISPLAY OPENFILE , , , "Select data file", _
                              EXE.PATH$, CHR$("All Files", 0, "*.*", 0),_
                              "Interpolated points.txt","" , 0 TO strFile
                              IF NOT ISFILE(strFIle) THEN
                              ? "File Not Found - Exiting"
                              EXIT FUNCTION
                              END IF

                              ' Now open this selected data file
                              ff = FREEFILE
                              OPEN strFile FOR INPUT AS #ff
                              WHILE NOT EOF(#ff)
                              INCR lngDatapoint
                              REDIM PRESERVE arrX(1 TO lngDatapoint)
                              REDIM PRESERVE arrY(1 TO lngDatapoint)
                              INPUT #ff, arrX(lngDatapoint) , arrY(lngDatapoint)

                              IF arrY(lngDatapoint) > MaxY THEN
                              MaxY = arrY(lngDatapoint)
                              END IF
                              IF arrY(lngDatapoint) < MinY THEN
                              MinY = arrY(lngDatapoint)
                              END IF
                              WEND
                              CLOSE #ff

                              ' sort x coordinates array in its ascending order
                              ' tagging along its y coordinates array
                              ARRAY SORT ArrX(), TAGARRAY ArrY()

                              Plotit ArrX(),ArrY(), MinY,MaxY, PATHNAME$(NAMEX, strFile )
                              END FUNCTION



                              ' Plot the graph with straight lines
                              ' if there are sufficient many points then the graph can smoothen
                              ' into curves
                              FUNCTION Plotit(ArrX() AS EXT, arrY() AS EXT ,minY AS EXT,_
                              MaxY AS EXT, strF AS STRING) AS LONG
                              LOCAL hWin AS DWORD
                              LOCAL minX,maxX AS EXT
                              LOCAL xc , yc AS EXT
                              LOCAL jj , deltaY AS LONG
                              LOCAL SMinX,SMinY,SMaxX,SMaxY AS EXT


                              LOCAL hFont1, hFont2, hFont3 AS LONG
                              FONT NEW "Lucida Console", 10, 0, 0, 0, 0 TO hFont1 ' set up fonts to be used later
                              FONT NEW "Lucida Console", 10, 0, 0, 0, 900 TO hFont2 ' vertical font (needs to be true type)
                              FONT NEW "Lucida Console", 2, 0, 0, 0, 0 TO hFont3 ' small font

                              'Get graph bounds
                              minX = ArrX(LBOUND(arrx()))
                              maxX = ArrX(UBOUND(arrx()))
                              'put 10% margins round data
                              SMinX = minx - (maxX-minX) *.1
                              sMinY = miny - (maxy-miny) *.1
                              sMaxX = maxx + (maxX-minX) *.1
                              sMaxY = maxy + (maxy-miny) *.1

                              GRAPHIC WINDOW NEW strF , 50, 50, 800, 500 TO hWin
                              GRAPHIC SCALE (Sminx,Smaxy) - (SmaxX,SMinY)

                              ' move the axes upwards with deltaY
                              ' so that can plot the graph better
                              deltaY = 2





                              ' display the X , Y axes labels
                              LOCAL xLab , yLab AS STRING
                              xLab = "X axis"
                              yLab = "Y axis"
                              GRAPHIC SET FONT hFont1
                              GPRINTX(xLab ,MaxX - (MaxX/2),MinY - 3.9, 1)
                              ' using the vertical fonts for the Y axis label
                              GRAPHIC SET FONT hFont2
                              GPRINTX(yLab , Minx-1.5,MAxY - (MaxY/2)+deltaY, 1)


                              ' for labeling the axes and scale marks
                              ' adapted from Dave Biggs
                              ' https://forum.powerbasic.com/forum/u...029#post601029
                              LOCAL Axmax, Aymax, Axmin, Aymin, xmark, ymark , r AS EXT
                              Axmax = MaxX
                              Aymax = MaxY
                              Axmin = MinX
                              Aymin = MinY
                              Axmin = IIF(Axmin > 10, FIX (Axmin/10)*10, FIX(Axmin)) ' round down minimum values &
                              Axmax = IIF(Axmax > 10, CEIL(Axmax/10)*10, CEIL(Axmax)) ' round up maximums for scale ranges
                              Aymin = IIF(Aymin > 10, FIX (Aymin/10)*10, FIX(Aymin))
                              Aymax = IIF(Aymax > 10, CEIL(Aymax/10)*10, CEIL(Aymax))

                              ' calculate mark spacing based on data range / values
                              xmark = SWITCH(Axmax-Axmin =< 10, 0.5, Axmax-Axmin =< 100, _
                              5, Axmax-Axmin =< 500, 50, Axmax-Axmin =< 1000, 100)
                              ymark = SWITCH(Aymax-Aymin =< 10, 0.5, Aymax-Aymin =< 100, _
                              5, Aymax-Aymin =< 500, 50, Aymax-Aymin =< 1000, 100)


                              'Draw the axes
                              GRAPHIC COLOR %BLACK,-1
                              GRAPHIC WIDTH 2
                              GRAPHIC LINE (Axmin,Aymin+deltaY)-(Axmin,Aymax+deltaY)
                              GRAPHIC LINE (Axmin,Aymin+deltaY)-(Axmax,Aymin+deltaY)


                              ' Marking the X axis
                              ' Note that we need to use trial and error method to adjust these labels
                              GRAPHIC SET FONT hFont1
                              GRAPHIC COLOR %BLUE, -1
                              FOR r = Axmin TO Axmax STEP xmark
                              ' draw the marks
                              GRAPHIC LINE (r, Aymin+ deltaY) - (r, Aymin-GRAPHIC(CELL.SIZE.Y)/2 + deltaY)
                              ' draw the mark labels
                              GRAPHIC SET POS STEP (-GRAPHIC(CELL.SIZE.X)*1.5, -GRAPHIC(CELL.SIZE.X)*5)
                              IF Axmax < 10 THEN
                              GRAPHIC PRINT FORMAT$(r, "0") ' print scale values
                              ELSE
                              GRAPHIC PRINT r
                              END IF
                              NEXT r


                              ' Marking the Y axis
                              ' Note that we need to use trial and error method to adjust these labels
                              FOR r = Aymin TO Aymax STEP ymark
                              ' draw the marks
                              GRAPHIC LINE (Axmin, r+ deltaY) - (Axmin-GRAPHIC(CELL.SIZE.X)/2, r+ deltaY)
                              ' draw the mark labels
                              GRAPHIC SET POS STEP (-GRAPHIC(CELL.SIZE.X)*5 +0.25 , GRAPHIC(CELL.SIZE.Y)/2)
                              IF Aymax =< 10 THEN
                              GRAPHIC PRINT FORMAT$(r, "0") ' print scale values
                              ELSE
                              GRAPHIC PRINT r
                              END IF
                              ' draw the marks
                              GRAPHIC LINE (Axmin-0.2, r+ deltaY) - (Axmin-GRAPHIC(CELL.SIZE.X)/2, r+ deltaY)
                              NEXT r


                              'Plot the data curve
                              ' sets the starting point first for the first line
                              GRAPHIC SET POS (arrx(1),arry(1)+deltaY)
                              FOR jj = 2 TO UBOUND(arrx())
                              GRAPHIC COLOR %CYAN,-1
                              GRAPHIC LINE - (arrX(jj),arrY(jj)+deltaY)
                              NEXT jj


                              ' draws the points with smaller font
                              GRAPHIC SET FONT hFont3
                              FOR jj = 2 TO UBOUND(arrx())
                              ' unfortunately the letter x is too large
                              ' GRAPHIC COLOR %RED,-1
                              ' xc = arrX(jj) - GRAPHIC(CELL.SIZE.X) / 2 ' adjust point positions to suit "X" center
                              ' yc = arrY(jj)+deltaY + GRAPHIC(CELL.SIZE.Y) / 2
                              ' GRAPHIC SET POS (xc, yc)
                              ' GRAPHIC PRINT "x"
                              GRAPHIC SET PIXEL (arrX(jj),arrY(jj)+deltaY), %RGB_DARKRED
                              NEXT jj

                              GRAPHIC WAITKEY$
                              GRAPHIC DETACH
                              GRAPHIC WINDOW END hWin
                              END FUNCTION


                              ' display the text
                              SUB GPRINTX(TX AS STRING,XS AS EXT,YS AS EXT, BGcol AS LONG)
                              IF BGCol = 1 THEN
                              'yellow color background
                              GRAPHIC COLOR %RGB_DARKGREEN, %RGB_LIGHTYELLOW
                              ELSE
                              ' default background
                              GRAPHIC COLOR %BLUE, -1
                              END IF
                              GRAPHIC SET POS(XS,YS)
                              GRAPHIC PRINT TX
                              END SUB

                              Comment


                              • #16
                                Click image for larger version

Name:	new plot.png
Views:	89
Size:	15.0 KB
ID:	812359

                                still the dots are too small , will need ellipses to show data points

                                Comment


                                • #17
                                  still the dots are too small , will need ellipses to show data points
                                  Code:
                                  'instead of -
                                  graphic set pixel (arrX(jj),arrY(jj)+deltaY), %rgb_darkred
                                  'try -
                                  graphic box (arrX(jj) - 1,arrY(jj)+deltaY - 1) - (arrX(jj) + 1,arrY(jj)+deltaY + 1), _
                                     50, %rgb_darkred, %rgb_darkred, 0
                                  For "Y" I may have the +1 and -1 reversed.

                                  If shape is a problem see Help about right and bottom, and/or change the "corner".
                                  Dale

                                  Comment


                                  • #18
                                    Thanks Dale that fix it, I have used

                                    ' Draw circular boxes with very tiny dimensions
                                    GRAPHIC BOX (arrX(jj) - 0.05,arrY(jj)+deltaY - 0.05) - _
                                    (arrX(jj) + 0.05,arrY(jj)+deltaY + 0.05), _
                                    100, %RGB_DARKRED, %RGB_DARKRED, 0

                                    Comment


                                    • #19
                                      Hello Mannish,
                                      Welcome to our forum, I noticed that you did not use Code tags to post your codes. The code tag button is located at the menu bar at #
                                      and when you click on this # button, you will get somethin like below
                                      now you just insert you code in between these tags, and everything would be ok, as shown below

                                      Code:
                                      ' Draw circular boxes with very tiny dimensions
                                              GRAPHIC BOX (arrX(jj) - 0.05,arrY(jj)+deltaY - 0.05) - _
                                                 (arrX(jj) + 0.05,arrY(jj)+deltaY + 0.05), _
                                                     100, %RGB_DARKRED, %RGB_DARKRED, 0
                                      you should also read this link for code indentations
                                      https://forum.powerbasic.com/forum/a...ing#post789625

                                      Comment


                                      • #20
                                        Mannish, you are starting where I did many years back. I needed a generic plotting program I could hook to my other tools. After twenty years, it has evolved considerably. The spline is the most recent addition. Here is a quick shot.

                                        You can create a DLL and call you routine from other tools which makes this nice.

                                        Some utilities

                                        here is a simple fast math function tool.
                                        Click image for larger version

Name:	Math Parser.png
Views:	61
Size:	23.8 KB
ID:	812380

                                        here I cut a portion of the plot and paste

                                        Click image for larger version

Name:	Copy Paste.png
Views:	48
Size:	41.3 KB
ID:	812381

                                        This is where I do data statistics on selections. The spline is the latest, but includes detrenders, smoothers, least squares fit, integration, differentiation etc.

                                        Click image for larger version

Name:	Data Statistics.png
Views:	49
Size:	11.3 KB
ID:	812382

                                        here is a discrete Fourier transform. Note the log-log plot. This was particularly challenging to work especially with zoom.

                                        Click image for larger version

Name:	Discrete Fourier Transform.png
Views:	49
Size:	50.3 KB
ID:	812383

                                        Comment

                                        Working...
                                        X