Announcement

Collapse
No announcement yet.

Cubic Spline

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

  • Mannish Bhandari
    replied
    Thanks Dean, your cubic spline prog is very good, I will try to add the graphics plotting by Stuart and Dave
    and list out here. I don't really need the sopshicated plotting by PV2D

    Leave a comment:


  • Dean Schrage
    replied
    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:	134
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:	102
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:	106
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:	104
Size:	50.3 KB
ID:	812383

    Leave a comment:


  • Tim Lakinir
    replied
    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

    Leave a comment:


  • Mannish Bhandari
    replied
    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

    Leave a comment:


  • Dale Yarker
    replied
    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".

    Leave a comment:


  • Mannish Bhandari
    replied
    Click image for larger version

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

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

    Leave a comment:


  • Mannish Bhandari
    replied
    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

    Leave a comment:


  • Mannish Bhandari
    replied
    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

    Leave a comment:


  • Mannish Bhandari
    replied
    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

    Leave a comment:


  • Mannish Bhandari
    replied
    Thanks Stuart, that's what I need

    Leave a comment:


  • Stuart McLachlan
    replied
    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

    Leave a comment:


  • Mannish Bhandari
    replied
    any suggestion to plot these points more accurately?

    Leave a comment:


  • Mannish Bhandari
    replied
    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






    Leave a comment:


  • Mannish Bhandari
    replied
    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

    Leave a comment:


  • Mannish Bhandari
    replied
    How to plot these interpolated data graphically?

    Leave a comment:


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

    Leave a comment:


  • Dean Schrage
    replied
    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

    Leave a comment:


  • Dean Schrage
    replied
    Sample output
    Click image for larger version

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

    Leave a comment:


  • Dale Yarker
    replied
    Re: Formatting code -
    Posted Source Code Losing Indentation/Formatting - PowerBASIC Peer Support Community
    (in FAQ section)

    Cheers,

    Leave a comment:


  • Dean Schrage
    replied
    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 = 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

    Leave a comment:

Working...
X