After running the program, please hit the "help" button for operating info.
Code: Select all
' curve fitting with polynomials in 2D
'
np_max=200 ' maximum # of points, may be modified
m_max=9 ' maximum degree of fitting polynomial
f$="ls_data" ' file used to store and restore point data
option base 1
dim xy(np_max,2),xlist$(np_max),ylist$(np_max)
dim coeff(m_max+1),coef$(m_max+1)
'
prog_init()
do slowdown
if b_p("new_x") then add_new_X()
k=list_selected("xp") ! if k>0 and k<=.np then edit_X(k)
k=list_selected("yp") ! if k>0 then edit_Y(k)
if b_p("del_p") then del_point()
if b_p("degree") then set_degree()
if b_p("calc") then calc()
if b_p("help") then help()
if b_p("load") then load_data(.f$)
if b_p("save") then save_data(.f$)
if b_p("stop") then stop
until forever
end
def prog_init()
graphics ! graphics clear ! draw color 0,0,0
get screen size sw,sh ! cx=sw/2
.gx1=20 ! .gy1=500 ! .gx2=748 ! .gy2=908
.np=0 ! help_page()
c_list("xp","X-coordinates",.xlist$,.np,20,20,190,420)
c_list("yp","Y-coordinates",.ylist$,.np,210,20,380,420)
button "new_x" text "Add X-value" at 20,440 size 170,30
button "del_p" text "Delete point" at 210,440 size 170,30
field "points" text " "&"0 Points" at 450,20 size 170,30 RO
field "points" font name "Courier-Bold"
button "calc" text "y = f(x)" at 450,440 size 290,30
button "degree" text "Degree" at 450,70 size 170,30
button "load" text "Load data" at 240,sh-50 size 120,30
button "save" text "Save data" at 400,sh-50 size 120,30
button "help" text "Help" at 20,sh-50 size 120,30
button "stop" text "Quit" at sw-140,sh-50 size 120,30
button "calc" hide
init_numpad(450,120,60,.7,.7,.7,1)
end def
def add_new_X()
if .np=.np_max then return
.np+=1 ! .xy(.np,2)=0 ! edit_X(.np)
field "points" text " " & .np & " Points"
list "xp" select .np ! list "xp" select -1
list "yp" select .np ! list "yp" select -1
end def
def edit_X(k)
.xy(k,1)=numpad(0,0)
tablesort (.np,2,.xy,1)
disp_lists() ! plot()
end def
def edit_Y(k)
.xy(k,2)=numpad(0,0)
for i=1 to .np ! .ylist$(i)=i & ": " &.xy(i,2) ! next i
c_list("yp","Y-coordinates",.ylist$,.np,210,20,380,420)
plot()
end def
def del_point()
k=numpad(1,.np)
if k<.np then
for i=k to .np-1
.xy(i,1)=.xy(i+1,1) ! .xy(i,2)=.xy(i+1,2)
next i
end if
.np-=1
field "points" text " " & .np & " Points"
disp_lists() ! plot()
end def
def disp_lists()
for i=1 to .np
.xlist$(i)=i & ": " &.xy(i,1) ! .ylist$(i)=i & ": " &.xy(i,2)
next i
c_list("xp","X-coordinates",.xlist$,.np,20,20,190,420)
c_list("yp","Y-coordinates",.ylist$,.np,210,20,380,420)
end def
def stat()
.xmin=1E+8 ! .xmax=-1E+8 ! .ymin=1E+8 ! .ymax=-1E+8
for i=1 to .np
x=.xy(i,1) ! .xmin=min(.xmin,x) ! .xmax=max(.xmax,x)
y=.xy(i,2) ! .ymin=min(.ymin,y) ! .ymax=max(.ymax,y)
next i
end def
def set_degree()
dim x(.np),y(.np)
if .np<2 then return
.m=numpad(0,min(.np-1,.m_max))
button "degree" text "Degree = "&.m
for i=1 to .np ! x(i)=.xy(i,1) ! y(i)=.xy(i,2) ! next i
ls (.np,.m,x,y,.coeff)
button "calc" text "y = f(x)" ! button "calc" show
c_table() ! plot()
end def
def calc()
x=numpad(.xmin,.xmax)
y=poly(.m,.coeff,x)
button "calc" text "x="&x&" => "&"y="&y
end def
def load_data(f$)
i=0
do
i+=1 ! file f$ input .xy(i,1),.xy(i,2)
until not data_exist(f$)
.np=i ! file f$ reset
field "points" text " " & .np & " Points"
disp_lists() ! plot()
end def
def save_data(f$)
file f$ delete
for i=1 to .np ! file f$ print .xy(i,1),.xy(i,2) ! next i
end def
def plot()
stat() ! if .xmin=.xmax or .ymin=.ymax or .np<2 then return
if .xmin>0 then
sc_x=720/.xmax ! y_as=.gx1+4 ! goto lab1
end if
if .xmax<0 then
sc_x=720/-.xmin ! y_as=.gx2-4 ! goto lab1
end if
sc_x=720/(.xmax-.xmin) ! y_as=.gx1+4-.xmin*sc_x
lab1:
if .ymin>0 then
sc_y=400/.ymax ! x_as=.gy2-4 ! goto lab2
end if
if .ymax<0 then
sc_y=400/-.ymin ! x_as=.gy1+4 ! goto lab2
end if
sc_y=400/(.ymax-.ymin) ! x_as=.gy2-4+.ymin*sc_y
lab2:
refresh off
fill color 0.88,0.85,0.65 ! fill rect .gx1,.gy1 to .gx2,.gy2
draw rect .gx1,.gy1 to .gx2,.gy2
draw color .7,.7,.7 ! draw size 1 ' draw grid
for x=.gx1+4 to .gx2-4 step 20
draw line x,.gy1+4 to x,.gy2-4
next x
for y=.gy1+4 to .gy2-4 step 20
draw line .gx1+4,y to .gx2-4,y
next y
draw color 0,0,1 ! draw size 3
draw line y_as,.gy1+4 to y_as,.gy2-4 ' draw axes
draw line .gx1+4,x_as to .gx2-4,x_as
draw color 0,0,0
for i=1 to .np
x=xp(.xy(i,1)) ! y=yp(.xy(i,2))
draw circle x,y size 2
next i
if .m>0 then
draw color 1,0,0 ! draw size 2
draw to y_as+sc_x*.xmin,x_as-sc_y*poly(.m,.coeff,.xmin)
dx=(.xmax-.xmin)/48
for x=.xmin+dx to .xmax step dx
draw line to y_as+sc_x*x,x_as-sc_y*poly(.m,.coeff,x)
next x
end if
refresh
draw color 0,0,0
end def
def xp(x)=plot.y_as+plot.sc_x*x
def yp(y)=plot.x_as-plot.sc_y*y
def c_table()
for i=1 to .m+1 ! .coef$(i)="a("&(i-1)&") = "&.coeff(i) ! next i
c_list("cf","Coefficients",.coef$,.m+1,450,120,620,420)
end def
def poly (n,coef(),x)
res=coef(n+1)
for i=n to 1 step-1
res=res*x+coef(i)
next i
return res
end def
def n2a$(num,lang,dec)
b$=" "
fac=10^dec
num$=floor(fac*num)/fac
tot=lang-len(num$)
if tot<1 then tot=1
rem test ("", tot)
a$=substr$(b$,1,tot) & num$
n2a$=a$
end def
def help()
n$="help" ! page n$ set ! page n$ alpha 0 ! page n$ show
for i=0 to 1 step 0.01 ! page n$ alpha i ! pause 0.01 ! next i
do slowdown ! until b_p("h_stop")
for i=1 to 0 step -0.01 ! page n$ alpha i ! pause 0.01 ! next i
page n$ hide ! page "" set
end def
def help_page()
xs=40 ! ys=200 ! w=680 ! h=680
h$="Curve fitting with polynomials in 2 dimensions."&chr$(10)&"A set of x-y points can be fitted with a polynomial of a suitable degree."&chr$(10)&"In this program a maximum number of 200 point is set, but this can easily be changed in the dim function."&chr$(10)&"The maximum degree of the polynomial is set at 9, as a higher degrees seldom add to the quality of the fit."&chr$(10)&"x- and y-coordinates are kept in separate lists with their point number (for synchronizing purpose). Both lists are kept in sequence of the x-coordinates."&chr$(10)&"Points are added by touching the ""new x"" button. the corresponding y-value is defaulted to zero."&chr$(10)&"Existing x- and y-values may be edited bij touching them in the appropriate list."&chr$(10)&"Points may be deleted by touching the ""delete point"" button and entering the point number."&chr$(10)&"The collection of points can be saved to a file and loaded later on by touching the ""save"" and ""load"" buttons."&chr$(10)&"A graph will appear as soon as two points are entered with different x- and y-values. Each modification in the point data is right away reflected in the graph."&chr$(10)&"Also, when more than one point is present, a fit command may be issued by touching the ""degree"" button and entering the degree of the fitting polynomial."&chr$(10)&"The coefficients of the polynomial are then shown and the calculated curve is drawn in the graph."&chr$(10)&"An y-value may be calculated for a given x-value by using the ""y=f(x)"" button."
n$="help" ! page n$ set ! page n$ frame xs,ys,w,h
field "hh" text h$ at 0,0 size w,h ML RO
field "hh" font color 0,0,1 ! field "hh" font size 20
button "h_stop" text "Close" at w/2-30,h-40 size 60,25
page n$ hide
page "" set
end def
' {LS-utils}
' id$ = object name
' cont$ = array met elementen
' size = aantal elementen in de list
'
def c_list(id$,title$,cont$(),size,xt,yt,xb,yb)
size=max(1,size)
dim temp$(size)
for i=1 to size ! temp$(i)=cont$(i) ! next i
list id$ text temp$ at xt+2,yt+32 size xb-xt-4,yb-yt-34
draw size 3
draw rect xt,yt to xb,yb ! draw line xt,yt+30 to xb,yt+30
draw color 0,0,1 ! draw text title$ at xt+5,yt+5
draw color 0,0,0
end def
' numeric table sorting function (option base 1 assumed)
' r = number of rows
' c = number of columns
' in(,) = table to be sorted
' col = column to be used for sorting
'
def tablesort (r,c,in(,),col)
dim sortcol(r),index(r),out(r,c)
for i=1 to r ! sortcol(i)=in(i,col) ! next i
sort sortcol as index
for i=1 to r ! for j=1 to c
out(i,j)=in(index(i),j)
next j ! next i
for i=1 to r ! for j=1 to c ! in(i,j)=out(i,j) ! next j ! next i
return
end def
' numerical keypad object
'
' produce a simple keypad to quickly enter a number in an app
' upon entry, the keypad disappears
' initialize once, multiple use after
' left upper corner is placed at "xtop,ytop"
' "bs" is the button size (keypad becomes 4.3 times larger)
' size of number is accepted between "minval" and "maxval"
' if both "minval" and "maxval" are zero, then no restrictions
' max number of tokens in the number is 10 (minus and dot included)
' works for option base 0 and 1
'
def init_numpad(xtop,ytop,bs,R,G,B,alpha)
name$="numpad" ! cn=10
page name$ set
page name$ frame xtop,ytop,0,0
set buttons custom
if bs<20 then bs=20
sp=4 ! th=.5*bs+4 ! ww=4*bs+5*sp ! hh=th+4*bs+6*sp
fsize=.5*bs
draw font size fsize ! set buttons font size fsize
draw color 1,1,1 ! fill color .7,.7,.7
button "rec" title "" at 0,0 size ww,hh
button "res" title "" at 0,0 size ww,th+4
fill color R,G,B ! fill alpha alpha
button "0" title "0" at sp,th+3*bs+5*sp size bs,bs
for k=1 to 9
x=(k-1)%3 ! y=2-floor((k-1)/3)
button k title k at (x+1)*sp+x*bs,th+y*bs+(y+2)*sp size bs,bs
next k
button "-" title "-" at 2*sp+bs,th+3*bs+5*sp size bs,bs
button "." title "." at 3*sp+2*bs,th+3*bs+5*sp size bs,bs
button "Cl" title "C" at 4*sp+3*bs,th+2*sp size bs,bs
button "del" title "<-" at 4*sp+3*bs,th+bs+3*sp size bs,bs
button "ok" title "ok" at 4*sp+3*bs,th+2*bs+4*sp size bs,2*bs+sp
page name$ hide
page name$ frame xtop,ytop,ww,hh
set buttons default ! set buttons font size 20
draw font size 20 ! draw color 0,0,0
end def
def numpad(minval,maxval)
page "numpad" set ! page "numpad" show
a$="" ! pflag=0 ! sflag=0 ! ob=1-option_base()
nump1:
if b_p("ok") then
number=val(a$) ! a$="" ! button "res" text ""
if minval<>0 or maxval<>0 then
if number<minval or number>maxval then
button "res" text "range error"
pflag=0 ! a$="" ! pause 1
button "res" text ""
goto nump1
end if
end if
page "numpad" hide ! page "" set
return number
end if
if b_p("Cl") then
a$ = "" ! pflag=0 ! sflag=0 ! goto nump3
end if
if b_p("del") and len(a$) then
ll=len(a$) ! if substr$(a$,ll-ob,ll-ob)="." then pflag=0
a$ = left$(a$,ll-1) ! sflag=0 ! goto nump3
end if
if b_p("-") then
a$ = "-" ! pflag=0 ! sflag=0 ! goto nump3
end if
if b_p(".") and not pflag and not sflag then
a$ &= "." ! pflag=1 ! goto nump3
end if
for k=0 to 9
t$=k
if b_p(t$) and not sflag then
a$ &= t$ ! goto nump3
end if
next k
goto nump1
nump3:
if len(a$)>10 then ! sflag=1 ! goto nump1 ! end if
button "res" text a$
goto nump1
end def
' this is a m-degree polynomial least squares fit of a number of
' point in 2 dimensional space
' there are n points, with x- and y-coordinates in vectors x() and y()
' m is the degree of the polynomial (take 1 for straight line fit,
' m=2 for parabola, and so on)
' the coefficients of the best fit polynomial are returned in vector c()
' f.i. for m=2 : y = c(1) + c(2)*x + c(3)*x^2
'
def ls (n,m,x(),y(),c())
dim a1(20,20),a2(20,20),a3(20,20),rl(20)
m=m+1
for i=1 to n ! a1(i,1)=1
for j=2 to m ! a1(i,j)=a1(i,j-1)*x(i) ! next j
next i
mat_trans (n,m,a1,a2)
mat_mul (m,n,a2,a1,a3)
mat_vec(m,n,a2,y,rl)
mat_inv (m,a3,a1)
mat_vec(m,n,a1,rl,c)
end def
def vec_uit (n$,n,v(),x,y,l,d)
ys=y-30
if n$<>"" then
draw text n$ at x,y
ys=ys+30
end if
for i=1 to n
a$=n2a$(i,3,0) & " " & n2a$(v(i),l,d)
draw text a$ at x-12,ys+25*i
next i
end def
def mat_trans (n,m,mat(,),matt(,))
for i=1 to n
for j=1 to m
matt(j,i)=mat(i,j)
next j
next i
end def
def mat_mul (n,m,mata(,),matb(,),matc(,))
for i=1 to n
for j=1 to n
tot=0
for k=1 to m
tot=tot+mata(i,k)*matb(k,j)
next k
matc(i,j)=tot
next j
next i
end def
def mat_vec (n,m,mat(,),vin(),vuit())
for i=1 to n
tot=0
for j=1 to m
tot=tot+mat( i,j)*vin(j)
next j
vuit(i)=tot
next i
end def
' produce the inverse of square matrix a() giving matrix ainv()
' returns 1 on success, 0 if matrix is singular
'
def mat_inv (nvar,a(,),ainv(,))
dim w(nvar,2*nvar)
for i=1 to nvar
for j=1 to nvar ! w(i,j)=a(i,j) ! w(i,j+nvar)=0 ! next j
w(i,i+nvar)=1
next i
for piv=1 to nvar
fac=w(piv,piv) ! if fac=0 then return 0
for j=piv to piv+nvar ! w(piv,j)=w(piv,j)/fac ! next j
for i=1 to nvar
if i<>piv then
fac=w(i,piv)
for j=piv to piv+nvar ! w(i,j)=w(i,j)-fac*w(piv,j) ! next j
endif
next i
next piv
for i=1 to nvar
for j=1 to nvar ! ainv(i,j)=w(i,j+nvar) ! next j
next i
return 1
end def
def b_p(a$) = button_pressed(a$)
- IMG_1487.PNG (936.15 KiB) Viewed 3285 times
)