|
This is working code written in Microsoft Visual FoxPro
5.0, and object oriented version of xBase.
* Wandering Turtle Defined IVM Tetrahedra
* by Kirby T. Urner, November 8, 1997
* coded in Visual FoxPro 5.0 (Microsoft)
* Last modified: December 26, 1997
*
* Three classes defined:
* Turtle an object moving in space, keeping track
* of its position in xyz and quadray coordinates
* as well as its distance from the origin
* Edger an object containing methods for operating on
* pairs of turtles as edge vertices
* Tetrahedron an object that keeping track of four vertices
* and knows how to compute the volume of a tetrahedron
* based on its six edge lengths
*
* Experiment:
*
* Let four turtles wander for i jumps, each time choosing one of
* the 12 IVM directions. Measure the 6 turtle interdistances and
* feed these to the tetrahedron's volume calculator.
*
* Let the turtles wander some more, and recompute, j time.
*
* Expectation:
*
* The tetrahedra will always have whole number volumes and will
* generally get bigger as the turtles tend to wander further apart
*
* Outcome:
*
* See Expectation.
*
#define root2 (2^.5)
otet = createobject("Tetrahedron")
? otet.ivmvol
* do walkit
release objects
return
procedure walkit
* run this experiment 5 times, having vertex turtles start
* from wherever they ended up last time
for j=1 to 5
* have the 4 vertex turtles wander -- 20 hops
i = 20
otet.init()
otet.oturtle1.randivm(i) && randivm is the 'hop i times' method
otet.oturtle2.randivm(i)
otet.oturtle3.randivm(i)
otet.oturtle4.randivm(i)
otet.calcivm()
? "Volume = " + str(otet.ivmvol,10,2)
endfor && repeat this volume experiment j times
endproc
define class turtle as custom
* class variables
distance = 0 && distance from origin
dimension xyzpos(3) && current position in xyz coords
dimension quadpos(4) && current position in quadrays
* when object is defined...
procedure init
this.reset()
endproc
* reset to origin
procedure reset
with this
store 0.0 to .xyzpos, .quadpos, .distance
endwith
endproc
* add a quad vector
procedure quadmov(a,b,c,d)
with this
.quadpos(1) = .quadpos(1) + a
.quadpos(2) = .quadpos(2) + b
.quadpos(3) = .quadpos(3) + c
.quadpos(4) = .quadpos(4) + d
.simplify()
.quad2xyz(a,b,c,d) && update xyz position
.distance = .quaddist() && update distance
endwith
endproc
* add an xyz vector
procedure xyzmov(x,y,z)
with this
.xyzpos(1) = .xyzpos(1) + x
.xyzpos(2) = .xyzpos(2) + y
.xyzpos(3) = .xyzpos(3) + z
.xyz2quad(x,y,z) && update quadray position
.distance = .xyzdist() && update distance
endwith
endproc
procedure quad2xyz(a,b,c,d)
with this
.xyzpos(1) = .xyzpos(1) + 1/root2 * (a - b - c + d)
.xyzpos(2) = .xyzpos(2) + 1/root2 * (a - b + c - d)
.xyzpos(3) = .xyzpos(3) + 1/root2 * (a + b - c - d)
endwith
endproc
procedure xyz2quad(x,y,z)
with this
.quadpos(1) = .quadpos(1) ;
+ 1/root2 * (iif(x>=0,x, 0)+iif(y>=0,y, 0)+iif(z>=0,z, 0))
.quadpos(2) = .quadpos(2) ;
+ 1/root2 * (iif(x>=0,0,-x)+iif(y>=0,0,-y)+iif(z>=0,z, 0))
.quadpos(3) = .quadpos(3) ;
+ 1/root2 * (iif(x>=0,0,-x)+iif(y>=0,y, 0)+iif(z>=0,0,-z))
.quadpos(4) = .quadpos(4) ;
+ 1/root2 * (iif(x>=0,x, 0)+iif(y>=0,0,-y)+iif(z>=0,0,-z))
.simplify()
endwith
endproc
* keep quadray coordinates in simplest form
procedure simplify
with this
local i
minval=.quadpos(1)
for i=1 to 4
minval = min(minval,.quadpos(i))
endfor
for i=1 to 4
.quadpos(i)=.quadpos(i)-minval
endfor
endwith
* compute distance using quadray distance formula
function quaddist
rtnval=0
with this
rtnval = ((3/2)*(.quadpos(1)^2 + .quadpos(2)^2 + ;
.quadpos(3)^2 + .quadpos(4)^2) - ;
( .quadpos(1)*.quadpos(2) + ;
.quadpos(1)*.quadpos(3) + ;
.quadpos(1)*.quadpos(4) ;
+ .quadpos(2)*.quadpos(3) + ;
.quadpos(2)*.quadpos(4) + .quadpos(3)*.quadpos(4)) )^0.5
endwith
return rtnval
endfunc
* compute distance using xyz distance formula
function xyzdist
rtnval=0
with this
rtnval = (.xyzpos(1)^2 + .xyzpos(2)^2 + .xyzpos(3)^2)^0.5
endwith
return rtnval
endfunc
* hop in one of 12 IVM directions 'nojumps' times (parameter)
procedure randivm(nojumps)
local pos2, pos0, rv(4)
for i = 1 to nojumps
rv = 1 && preset (a,b,c,d) to (1,1,1,1)
* pick a number between 1 and 4
pos2 = mod(1+int(1000*rand()),4)+1 && pos2 is where the 2 will be
pos0 = pos2 && pos0 is where the 0 will be
do while pos2=pos0
pos2 = mod(1+int(1000*rand()),4)+1
enddo
rv(pos2) = 2 && one of the coordinates becomes a 2
rv(pos0) = 0 && another, different coord becomes a 0
* the remaining 2 coordinates remain set at 1
this.quadmov(rv(1),rv(2),rv(3),rv(4))
endfor
endproc
* display current position information
procedure saypos
with this
? this.name
?
? "XYZ = (" + ltrim(str(.xyzpos(1),5,2)) + ", " ;
+ ltrim(str(.xyzpos(2),5,2)) + ", " ;
+ ltrim(str(.xyzpos(3),5,2)) + ")"
? "QUAD = (" + ltrim(str(.quadpos(1),5,2)) + ", " ;
+ ltrim(str(.quadpos(2),5,2)) + ", " ;
+ ltrim(str(.quadpos(3),5,2)) + ", " ;
+ ltrim(str(.quadpos(4),5,2)) + ")"
?
? "XYZDIST = "+ str(.xyzdist(),5,2)
? "QUADDIST= "+ str(.quaddist(),5,2)
endwith
endproc
enddefine
define class tetrahedron as custom
add object oturtle1 as turtle
add object oturtle2 as turtle
add object oturtle3 as turtle
add object oturtle4 as turtle
add object oedger as edger
dimension elengths(6) && edgelengths
a=0
b=0
c=0
d=0
e=0
f=0
a2=0
b2=0
c2=0
d2=0
e2=0
f2=0
xyzvol=0
ivmvol=0
procedure init
this.oturtle1.quadmov(1,0,0,0)
this.oturtle2.quadmov(0,1,0,0)
this.oturtle3.quadmov(0,0,1,0)
this.oturtle4.quadmov(0,0,0,1)
this.elengths=1
this.setlengths()
this.calcivm()
this.calcxyz()
endproc
procedure setverts(matrix)
with this
.oturtle1.reset()
.oturtle2.reset()
.oturtle3.reset()
.oturtle4.reset()
.oturtle1.quadmov(matrix(1,1),matrix(1,2),matrix(1,3),matrix(1,4))
.oturtle2.quadmov(matrix(2,1),matrix(2,2),matrix(2,3),matrix(2,4))
.oturtle3.quadmov(matrix(3,1),matrix(3,2),matrix(3,3),matrix(3,4))
.oturtle4.quadmov(matrix(4,1),matrix(4,2),matrix(4,3),matrix(4,4))
.setlengths()
.calcivm()
.calcxyz()
endwith
endproc
procedure setlengths
* a,b,c radiate from a common vertex with triangle
* def forming the opposite triangle, giving closed
* triangles abd, bce, acf, def. Pairs of opposite
* edges are ae, cd, bf.
with this
.oedger.setpos1(.oturtle1.quadpos(1),;
.oturtle1.quadpos(2),;
.oturtle1.quadpos(3),;
.oturtle1.quadpos(4))
.oedger.setpos2(.oturtle2.quadpos(1),;
.oturtle2.quadpos(2),;
.oturtle2.quadpos(3),;
.oturtle2.quadpos(4))
.elengths(1)=.oedger.diffdist
.oedger.setpos2(.oturtle3.quadpos(1),;
.oturtle3.quadpos(2),;
.oturtle3.quadpos(3),;
.oturtle3.quadpos(4))
.elengths(2)=.oedger.diffdist
.oedger.setpos2(.oturtle4.quadpos(1),;
.oturtle4.quadpos(2),;
.oturtle4.quadpos(3),;
.oturtle4.quadpos(4))
.elengths(3)=.oedger.diffdist
.oedger.setpos1(.oturtle2.quadpos(1),;
.oturtle2.quadpos(2),;
.oturtle2.quadpos(3),;
.oturtle2.quadpos(4))
.oedger.setpos2(.oturtle3.quadpos(1),;
.oturtle3.quadpos(2),;
.oturtle3.quadpos(3),;
.oturtle3.quadpos(4))
.elengths(4)=.oedger.diffdist
.oedger.setpos2(.oturtle4.quadpos(1),;
.oturtle4.quadpos(2),;
.oturtle4.quadpos(3),;
.oturtle4.quadpos(4))
.elengths(6)=.oedger.diffdist
.oedger.setpos1(.oturtle3.quadpos(1),;
.oturtle3.quadpos(2),;
.oturtle3.quadpos(3),;
.oturtle3.quadpos(4))
.elengths(5)=.oedger.diffdist
.a=.elengths(1)/2
.b=.elengths(2)/2
.c=.elengths(3)/2
.d=.elengths(4)/2
.e=.elengths(5)/2
.f=.elengths(6)/2
endwith
endproc
procedure calcxyz
this.setlengths()
this.power2()
this.xyzvol = 1/12 * (this.addopen()-this.addclosed()-this.addopposite())^0.5
return
procedure calcivm
this.setlengths()
this.power2()
this.ivmvol = ((this.addopen()-this.addclosed()-this.addopposite())/2)^0.5
return
* do all 2nd powering in one place, for simplicity
protected procedure power2
with this
.a2=.a*.a
.b2=.b*.b
.c2=.c*.c
.d2=.d*.d
.e2=.e*.e
.f2=.f*.f
endwith
return
* open triangles term
function addopen
local sumval
with this
sumval = .f2*.a2*.b2
sumval = sumval + .d2*.a2*.c2
sumval = sumval + .a2*.b2*.e2
sumval = sumval + .c2*.b2*.d2
sumval = sumval + .e2*.c2*.a2
sumval = sumval + .f2*.c2*.b2
sumval = sumval + .e2*.d2*.a2
sumval = sumval + .b2*.d2*.f2
sumval = sumval + .b2*.e2*.f2
sumval = sumval + .d2*.e2*.c2
sumval = sumval + .a2*.f2*.e2
sumval = sumval + .d2*.f2*.c2
endwith
return sumval
endfunc
* closed triangles term
function addclosed
local sumval
with this
sumval = .a2*.b2*.d2
sumval = sumval + .d2*.e2*.f2
sumval = sumval + .b2*.c2*.e2
sumval = sumval + .a2*.c2*.f2
endwith
return sumval
endfunc
* opposite triangles term
function addopposite
local sumval
with this
sumval = .a2*.e2*(.a2+.e2)
sumval = sumval + .b2*.f2*(.b2+.f2)
sumval = sumval + .c2*.d2*(.c2+.d2)
endwith
return sumval
endfunc
enddefine
define class edger as custom
* include two Turtle objects as part of class definition
add object vertex1 as turtle
add object vertex2 as turtle
vert1dist = 0 && length of v1
vert2dist = 0 && length of v2
diffdist = 0 && length of difference
angle = 0 && angle between v1 and v2 in radians
xyzArea = 0 && area between 2 vectors (parallelogram)
Area = 0 && area between 2 vectors (triangle) with
&& equiangular triangle 1 x 1 = 1
procedure init
this.vertex1.reset()
this.vertex2.reset()
endproc
procedure destroy
release object vertex1
release object vertex2
endproc
procedure setpos1(a,b,c,d)
with this
.vertex1.reset()
.vertex1.quadmov(a,b,c,d)
.getdiffdist()
.vert1dist = .vertex1.distance
if .vert1dist>0 and .vert2dist>0
.getangle()
.getArea()
.getxyzArea()
endif
endwith
endproc
procedure setpos2(a,b,c,d)
with this
.vertex2.reset()
.vertex2.quadmov(a,b,c,d)
.vert2dist = .vertex2.distance
.getdiffdist()
if .vert1dist>0 and .vert2dist>0
.getangle()
.getArea()
.getxyzArea()
endif
endwith
endproc
procedure getangle
local i, temp(4), vsum(4), abdist
with this
for i = 1 to 4
vsum(i) = (.vertex1.quadpos(i)/.vertex1.distance ;
- .vertex2.quadpos(i)/.vertex2.distance)
temp(i)=.vertex1.quadpos(i)
endfor
.vertex1.reset()
.vertex1.quadmov(vsum(1),vsum(2),vsum(3),vsum(4))
abdist =.vertex1.distance
.vertex1.reset()
.vertex1.quadmov(temp(1),temp(2),temp(3),temp(4))
endwith
.angle = 2 * asin(.5 * abdist)
endproc
procedure getXyzArea
.xyzArea = .vert1dist/2 * .vert2dist/2 * sin(.angle)
endproc
procedure getArea
.Area = (4/3)^0.5 * .vert1dist/2 * .vert2dist/2 * sin(.angle)
endproc
procedure getdiffdist
local i, temp(4), vsum(4)
with this
for i = 1 to 4
vsum(i) = .vertex1.quadpos(i) - .vertex2.quadpos(i)
temp(i)=.vertex1.quadpos(i)
endfor
.vertex1.reset()
.vertex1.quadmov(vsum(1),vsum(2),vsum(3),vsum(4))
.diffdist =.vertex1.distance
.vertex1.reset()
.vertex1.quadmov(temp(1),temp(2),temp(3),temp(4))
endwith
endproc
enddefine
Synergetics on the
Web |