Re: Matrix Algebra

new topic     » goto parent     » topic index » view thread      » older message » newer message

I needed to quickly put together Vector Add/Subtract/Dot/X product.
Created the following but it's not error free; I am still working on=20
ironing out some details; the arccos/arcsin would probably be better=20
done through an external call to 'C' function, I'll wait for that
feature. Having done some work in APL, I think there is alot of value
to euphoria in matrix manipulations; I was impressed by my ability to
quickly generate a 64 billion 3D array of integers, something I would be
hardpressed to do in 'C'. The special ascii characters are: 224,225,231
for
alpha, beta, gamma unit angles on the vectors.

-WWW  (that's not world wide web, just my initials)

include get.e
include graphics.e
include matx.e

--Calculate magnitude of 3d vector
function SizeVect(sequence v)
    atom siz
    sequence r1
   =20
    r1 =3D power(v,2)
    siz =3D power(r1[1] + r1[2] + r1[3],.5)
    return siz
end function

--Print sequence vector in ijk format (with highlight)
procedure PrintVect(sequence v, integer r, integer c)
sequence sign
    sign =3D " + "
    position(r,c)
    printf(1,"%g",v[1])
    text_color(BRIGHT_CYAN)
    puts(1,"i")
    text_color(WHITE)
    if v[2] < 0 then
        sign =3D " - "
        v[2] =3D v[2] * -1
    end if
    printf(1,sign & "%g",v[2])
    text_color(BRIGHT_CYAN)
    puts(1,"j")
    text_color(WHITE)
    sign =3D " + "
    if v[3] < 0 then
        sign =3D " - "
        v[3] =3D v[3] * -1
    end if
   =20
    printf(1,sign & "%g",v[3])
    text_color(BRIGHT_CYAN)
    puts(1,"k")
    text_color(WHITE)
end procedure

--Provide data retrieval
procedure Xtern(atom in)
    integer fn
    fn =3D open("tst.in","w")
    printf(fn,"%g",in)
    close(fn)
end procedure

sequence ans, mat, mat1, mat2, wrk
atom NOTDONE, rc, scalar,  w3, w4
atom t

NOTDONE =3D 1
while NOTDONE do
clear_screen()
position(20,1) =20
puts(1,"Enter 3D vector A sequence: {i,j,k}  0 to end->")
ans =3D get(0)
if atom(ans[2]) then
    exit
end if
rc =3D length(ans[2])
if rc =3D 0 then
    NOTDONE =3D 0
else

mat1 =3D ans[2]
w3 =3D SizeVect(mat1)
position(2,1)
puts(1,"Vector A is: ")
PrintVect(mat1,2,15)
printf(1," mag=3D%g",w3)
printf(1," =E0=3D%g =E1=3D%g =E7=3D%g",{arccos(mat1[1]/w3),arcsin(mat1[2]=
/w3),
                                arccos(mat1[3]/w3)})

position(21,1)
puts(1,"Enter 3D vector B sequence ->")
ans =3D get(0)
mat2 =3D ans[2]
w4 =3D SizeVect(mat2)

position(3,1)
puts(1,"Vector B is: ")
PrintVect(mat2,3,15)
printf(1," mag=3D%g",w4)
printf(1," =E0=3D%g =E1=3D%g =E7=3D%g",{arccos(mat2[1]/w4),arcsin(mat2[2]=
/w4),
                                arccos(mat2[3]/w4)})

-- Vector manipulations -----
-- Vector Add
mat =3D mat1 + mat2
position(5,1)
puts(1,"^A + ^B =3D ")
PrintVect(mat,5,15)
=20
-- Vector Subtract -------
mat =3D mat1 - mat2
position(6,1)
puts(1,"^A - ^B =3D ")
PrintVect(mat,6,15)

-- Dot Product ------
mat =3D mat1 * mat2
scalar =3D mat[1] + mat[2] + mat[3]
position(7,1)
puts(1,"^A =F8 ^B =3D ")  -- graphic dot is 248 decimal
printf(1,"%g",scalar)

-- Cross product ----
wrk =3D repeat(0,6)     -- initialize work matrix
wrk[1] =3D mat1[1] * mat2[2]
wrk[2] =3D mat1[1] * mat2[3] * -1
wrk[3] =3D mat1[2] * mat2[1] * -1
wrk[4] =3D mat1[2] * mat2[3]
wrk[5] =3D mat1[3] * mat2[1]
wrk[6] =3D mat1[3] * mat2[2] * -1

mat[1] =3D wrk[4] + wrk[6]  -- new i
mat[2] =3D wrk[2] + wrk[5]  -- new j
mat[3] =3D wrk[1] + wrk[3]  -- new k

position(8,1)
puts(1,"^A X ^B =3D ")
PrintVect(mat,8,15)

t =3D scalar / (w3 * w4)
Xtern(t)
position(9,1)
printf(1,"Angle =E9 =3D %g degrees", arccos(scalar / (w3 * w4)))

rc =3D wait_key()
end if
end while
clear_screen()
puts(1,"--The End--")

---- Include file matx.e -----
global constant PI        =3D  3.14159265
--constant PI2      =3D  1.570796326794896619231322    -- pi / 2
--constant PI3      =3D  1.047197551196597746154214   -- pi / 3
--constant PI4      =3D  0.7853981633974483096156608  -- pi / 4
--constant SQRTPI   =3D  1.772453850905516027298167   -- sqrt(pi)
--constant SQRT2PI  =3D  2.506628274631000502415765   -- sqrt(2.pi)
--constant TWOPI    =3D  6.283185307179586476925287   -- 2.pi
--constant LNPI     =3D  1.144729885849400174143427   -- ln(pi)
--constant LOGPI    =3D  0.4971498726941338543512683  -- log(pi)
--constant LOGE     =3D  0.4342944819032518276511289  -- log(e)
--constant LN10     =3D  2.302585092994045684017991   -- ln(10)
global constant E         =3D  2.718281828459045235360287   -- exp(1)
global constant ONERAD   =3D 57.295779513082320876798155   -- 1 rad in =F8
global constant ONEDEG   =3D  0.017453292519943295769237   -- 1=F8 in rad
sequence cos_angls
sequence sin_angls
atom cos_init
atom sin_init

cos_init =3D 0
sin_init =3D 0

--Generate an integer sequence
-- behaves as APL interval function
global function interval(atom lim)
    sequence n

    n =3D {}
    puts(1,"Sequence initialized\n")
    for x=3D1  to lim do
        n =3D append(n,x)
    end for
    return n
end function

global function stripCR(sequence in)
    sequence wrk
    wrk =3D in[1..(length(in)-1)]
    return wrk
end function

--arc trig functions do a onetime build of sin/cos table
-- then extrapolates to return angle in degrees=20
-- compute arccos, return angle in degrees whose cosine is theta
global function arccos(atom in)=20
atom sel, wrk1, wrk2, wrk3
if cos_init =3D 0 then
    cos_init =3D 9
    cos_angls =3D repeat(0.0,360)
    for x=3D1  to 360  do
        cos_angls[x] =3D cos(x * ONEDEG)
        if x =3D 90 or x =3D 270 then
            cos_angls[x] =3D 0
        end if
    end for
end if
sel =3D 1=20
while in < cos_angls[sel] do
    sel =3D sel + 1
end while=20
sel =3D sel - 1
wrk1 =3D in - cos_angls[sel]
wrk2 =3D cos_angls[sel+1]-cos_angls[sel]
wrk3 =3D sel + wrk1 / wrk2
return wrk3
end function

-- compute arcsin, return angle in degrees whose sine is theta
global function arcsin(atom in)=20
atom sel, wrk1, wrk2, wrk3
if sin_init =3D 0 then
    sin_init =3D 9
    sin_angls =3D repeat(0.0,360)
    for x=3D1  to 360  do
        sin_angls[x] =3D sin(x * ONEDEG)
        if x =3D 90 or x =3D 360 then
            sin_angls[x] =3D 0
        end if
    end for
end if
sel =3D 1=20
while in > sin_angls[sel] do
    sel =3D sel + 1
    if sin_angls[sel] < 0 then
        in =3D -in
    else
        if sin_angls[sel-1] <=3D 0 then
            in =3D -in
        end if     =20
    end if
end while=20
sel =3D sel - 1
if sel =3D 0 then
    sel =3D 359
end if
wrk1 =3D in - sin_angls[sel]
wrk2 =3D sin_angls[sel+1]-sin_angls[sel]
wrk3 =3D sel + wrk1 / wrk2
return wrk3
end function

new topic     » goto parent     » topic index » view thread      » older message » newer message

Search



Quick Links

User menu

Not signed in.

Misc Menu