Code: Select all
' Library "matlib" version 4-01-2014
' vector and matrix calculations and some additional thingies.
' produced by "Henko"
' no copyrights and no guaranties
'
' IMPORTANT: all functions are based upon OPTION BASE 1
'
' vec_in (n,v,x,y)
' vec_out (n$,n,v,x,y,l,d)
' vec_zero (n,v)
' vec_one (n,v)
' vec_rnd (n,v,mini,maxi)
' vec_copy (n,from,into)
' vec_scal (n,v,s)
' vec_len (n,v)
' vec_norm (n,v)
' vec_plus (n,v1,v2)
' vec_min (n,v1,v2)
' inprod (n,v1,v2)
' mat_in (n,m,mat,x,y)
' mat_out (n$,n,m,mat,x,y,l,d)
' mat_zero (n,m,mat)
' mat_unit (n,mat)
' mat_det(n,mat)
' mat_rnd (n,m,mat,mini,maxi)
' mat_scal (n,m,mat,s)
' mat_copy (n,m,from,into)
' mat_trans (n,m,mat,matt)
' mat_plus (n,m,mata,matb)
' mat_min (n,m,mata,matb)
' mat_mul (n,m,mata,matb,matc)
' mat_inv (n,mata,ainv)
' mat_vec (n,m,mat,vin,vout)
' eigen_2 (mat,ev,lab)
' eigen_n (n,mat,ev)
' ls (n,m,x,y,c)
' rot2(deg,vin(),vout())
' rot_x(deg,vin(),vout())
' rot_y(deg,vin(),vout())
' rot_z(deg,vin(),vout())
' rot_xyz(dgx,dgy,dgz,vin(),vout())
' poly (n,coef,x)
' surf_n(n,c(,))
' surf_3(a(),b())
' mod (a,m)
' rest (a,b)
' def pi
' def rad
' input a n-sized vector v() at screen position x,y
'
def vec_in (n,v(),x,y)
dim b$(20)
for i=1 to n
b$(i)="a" & i ! field b$(i) text b$(i) at x,y+30*(i-1) size 80,25
next i
som=0
loop:
for i=1 to n
if field_changed(b$(i)) then
v(i)=field_text$(b$(i)) ! field b$(i) delete
draw text n2a$(v(i),6,2) at x,y+30*(i-1)
som=som+1
end if
next i
if som<n then goto loop
end def
' print a vector at position x,y with number length l and precision d
'
def vec_out (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
' produce a vector filled with zero's
'
def vec_zero (n,v())
for i=1 to n ! v(i)=0 ! next i
end def
' produce a vector filled with one's
'
def vec_unit (n,v())
for i=1 to n ! v(i)=1 ! next i
end def
' produce a vector filled with random numbers between mini and maxi
'
def vec_rnd (n,v(),mini,maxi)
for i=1 to n ! v(i)=mini+(maxi-mini)*rnd(1) ! next i
end def
' copy vector from() into vector into()
'
def vec_copy (n,from(),into())
for i=1 to n ! into(i)=from(i) ! next i
end def
' multiply a vector with scalar s
'
def vec_scal (n,v(),s)
for i=1 to n ! v(i)=s*v(i) ! next i
end def
' produce the length of a vector
'
def vec_len (n,v()) = sqrt(inprod(n,v,v))
' normalize a vector to unit length
'
def vec_norm (n,v())
vec_scal(n,v,1/vec_len(n,v))
end def
' add vector v2() to vector v1()
'
def vec_plus (n,v1(),v2())
for i=1 to n ! v1(i)=v1(i)+v2(i) ! next i
end def
' subtract vector v2() from vector v1()
'
def vec_min (n,v1(),v2())
for i=1 to n ! v1(i)=v1(i)-v2(i) ! next i
end def
' produce the scalar product of two vectors
'
def inprod (n,v1(),v2())
som=0
for i=1 to n ! som=som+v1(i)*v2(i) ! next i
inprod=som
end def
' input a nxm matrix at position x,y (left upper corner)
' n is the number of rows, m is the number of columns
'
def mat_in (n,m,mat(,),x,y)
dim b$(100)
b$(1)=0
for i=1 to n
for j=1 to m
k=m*(i-1)+j ! b$(k)=" a" & i & j ! tt$=b$(k)
field tt$ text tt$ at x+100*(j-1),y+30*(i-1) size 80,25
next j
next i
som=0
loop1:
for i=1 to n
for j=1 to m
k=m*(i-1)+j
if field_changed(b$(k)) then
mat(i,j)=field_text$(b$(k))
field b$(k) delete
draw text n2a$(mat(i,j),6,2) at x+100*(j-1),y+30*(i-1)
som=som+1
end if
next j
next i
if som<n*m then goto loop1
end def
' produce a nxm matrix with zero elements
'
def mat_zero (n,m,mat(,))
for i=1 to n
for j=1 to m ! mat(i,j)=0 ! next j
next i
end def
' produce a unit matrix with one's in the diagonal and zeros elsewhere
'
def mat_unit (n,mat(,))
for i=1 to n
for j=1 to n ! if i=j then mat(i,j)=1 else mat(i,j)=0 ! next j
next i
end def
' produce the determinant of a nxn matrix
' the input matrix mat(,) remains unchanged
' the value of the determinant is given back by the function
'
def mat_det(n,mat(,))
dim a(n,n)
mat_copy(n,n,mat,a)
for i=1 to n-1
if i=1 then det=a(i,i) else det=det*a(i,i)
for j=i+1 to n
fac=a(j,i)/a(i,i)
for k=i+1 to n ! a(j,k)=a(j,k)-fac*a(i,k) ! next k
next j
next i
mat_det=a(n,n)*det
end def
' produce a matrix with random elements between mini and maxi
'
def mat_rnd (n,m,mat(,),mini,maxi)
for i=1 to n
for j=1 to m ! mat(i,j)=mini+(maxi-mini)*rnd(1) ! next j
next i
end def
' multiply all elements of a matrix with scalair s
'
def mat_scal (n,m,mat(,),s)
for i=1 to n
for j=1 to m ! mat(i,j)=s*mat(i,j) ! next j
next i
end def
' copy matrix from() into matrix into()
'
def mat_copy (n,m,from(,),into(,))
for i=1 to n
for j=1 to m ! into(i,j)=from(i,j) ! next j
next i
end def
' produce the transpose of matrix mat() into matrix matt()
'
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
' add matb() to mata()
'
def mat_plus (n,m,mata(,),matb(,))
for i=1 to n
for j=1 to m ! mata( i,j)=mata(i,j)+matb(i,j) ! next j
next i
end def
' subtract matb() from mata()
'
def mat_min (n,m,mata(,),matb(,))
for i=1 to n
for j=1 to m ! mata( i,j)=mata(i,j)-matb(i,j) ! next j
next i
end def
' produce product of mata() and matb() giving matc()
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
' produce the inverse of square matrix a() giving matrix ainv()
'
def mat_inv (nvar,a(,),ainv(,))
dim w(25,50)
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)
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
end def
' print a nxm matrix at position x,y left upper corner
' each element having total length l and d decimals
' n$ is an text, printed to identify the matrix (may be empty)
'
def mat_out (n$,n,m,mat(,),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
for j=1 to m
a$=n2a$(mat(i,j),l,d) ! draw text a$ at x+12*l*(j-1),ys+25*i
next j
next i
end def
' produce product of a matrix and a vector vin(), giving vout()
' the matrix has size nxm, the vin() has size m, vout() has size n
'
def mat_vec (n,m,mat(,),vin(),vout())
dim v(20)
for i=1 to n
tot=0
for j=1 to m ! tot=tot+mat( i,j)*vin(j) ! next j
v(i)=tot
next i
for i=1 to m ! vout(i)=v(i) ! next i
end def
' eigenvalues and eigenvectors of a 2x2 matrix
' matrix has real coefficients
' matrix is passed as "mat"
' function returns 0 if no real eigenvalues are found,
' else the function returns 1
' 2 eigenvalues are returned in the vector "lab"
' 2 eigenvectors are returned as colums in the matrix "ev"
' eigenvectors are normalized to a length of 1
' library "matlib" is needed
'
def eigen_2 (mat(,),ev(,),lab())
dim x(2)
eigen_2=0
discr=(mat(1,1)-mat(2,2))^2+4*mat(1,2)*mat(2,1)
if discr<0 then return
eigen_2=1
discr=sqrt(discr) ! s=mat(1,1)+mat(2,2)
lab(1)=(s+discr)/2 ! lab(2)=(s-discr)/2
for j=1 to 2
x(1)=1
if mat(1,2) then
x(2)=(lab(j)-mat(1,1))/mat(1,2)
else
x(2)=(lab(j)-mat(2,1))/mat(2,2)
end if
vec_norm(2,x)
ev(1,j)=x(1) ! ev(2,j)=x(2)
next j
end def
' Find largest eigenvalue with eigenvector for nxn matrix,
' using the simplest "power method".
' No results in case of complex or multiple eigenvalues, the
' function will then return a value of 0.
' If the iteration converges, the function returns the eigenvalue
' and the accompanying eigenvector in the vector "ev"
'
def eigen_n (n,mat(,),ev())
dim evo(n)
count=0 ! maxcount=100 ! eps=.00001
labo=1 ! vec_rnd(n,evo,-1,1)
do
mat_vec(n,n,mat,evo,ev)
lab=vec_len(n,ev)/vec_len(n,evo)
dif=abs(abs(lab)-abs(labo)) ! vec_norm(n,ev)
if dif>eps then
labo=lab ! vec_copy(n,ev,evo) ! count=count+1
end if
until dif<eps or count=maxcount
if ev(1)*evo(1)<0 then lab=-lab
if count=maxcount then eigen_n=0 else eigen_n=lab
end def
' this is a n-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
' rotate a 2-dimensional vector vin() deg degrees, giving vector vout()
'
def rot2(deg,vin(),vout())
x=vin(1) ! vout(1)=x*cos(deg)-vin(2)*sin(deg)
vout(2)=x*sin(deg)+vin(2)*cos(deg)
end def
' rotate the 3-dim. vector vin() deg degrees about the x-axes -> vout()
'
def rot_x(deg,vin(),vout())
vout(1)=vin(1)
y=vin(2) ! vout(2)=y*cos(deg)-vin(3)*sin(deg)
vout(3)=y*sin(deg)+vin(3)*cos(deg)
end def
' rotate the 3-dim. vector vin() deg degrees about the y-axes -> vout()
'
def rot_y(deg,vin(),vout())
vout(2)=vin(2)
x=vin(1) ! vout(1)=x*cos(deg)-vin(3)*sin(deg)
vout(3)=x*sin(deg)+vin(3)*cos(deg)
end def
' rotate the 3-dim. vector vin() deg degrees about the z-axes -> vout()
'
def rot_z(deg,vin(),vout())
vout(3)=vin(3)
x=vin(1) ! vout(1)=x*cos(deg)-vin(2)*sin(deg)
vout(2)=x*sin(deg)+vin(2)*cos(deg)
end def
' rotate vector vin() about all 3 axes, giving vector vout()
def rot_xyz(dgx,dgy,dgz,vin(),vout())
dim temp(3)
rot_x(dgx,vin,temp)
rot_y(dgy,temp,temp)
rot_z(dgz,temp,vout)
end def
' calculate the value of a polynomial for a give value of x
' n id the degree of the polynomial, hence n+1 coefficients must
' be passed: a0, a1, a2, ..... an, in that order
'
def poly (n,coef(),x)
res=coef(n+1)
for i=n to 1 step-1 ! res=res*x+coef(i) ! next i
poly=res
end def
' calculates the surface of a n-polygon, n>=3
' (divides the polygon in triangles and then uses surf_3 function)
' c() has size nx2 and contains the coordinates of the vertices.
'
def surf_n(n,c(,))
dim cr(20,2),u(2),v(2)
if n<3 then return
x1=c(1,1) ! y1=c(1,2) ! sum=0
for i=1 to n ! cr(i,1)=c(i,1)-x1 ! cr(i,2)=c(i,2)-y1 ! next i
for i=2 to n-1
u(1)=cr(i,1) ! u(2)=cr(i,2) ! v(1)=cr(i+1,1) ! v(2)=cr(i+1,2)
sum=sum+surf_3(u,v)
next i
surf_n=sum
end def
' calculates the surface of a triangle, produced by 2 vectors
'
def surf_3(a(),b())=.5*abs(a(1)*b(2)-a(2)*b(1))
' integer remainder of a/m
'
def mod(a,m)
d=a/m
mod=m*(d-int(d))
end def
' remainder of a/m
'
def rest (a,b)
q=a/b
rest=q-int(q)
end def
' value of pi
'
def pi=3.14159265
' value of 1 radian
'
def rad=180/pi
' formatting floats for graphics mode
'
def n2a$(num,lang,dec)
b$=" "
fac=10^dec
num$=int(fac*num)/fac
tot=lang-len(num$)
if tot<1 then tot=1
a$=substr$(b$,1,tot) & num$
n2a$=a$
end def
' pre-padding strings with blancs to a given width
'
def pre_pad$(w,a$)
sp$=" "
tot=w-len(a$)
pre_pad$=substr$(sp$,1,tot) & a$
end def
