$SPAD/src/input images8a.input

The Axiom Team December 3, 2016

Abstract

Contents

2

1 License

— license —

--Copyright The Numerical Algorithms Group Limited 1994.

———-

— * —

\getchunk{license}
-- Sierpinsky's Tetrahedron

-- Sierpinsky's Tetrahedron
-- Bring DH matrices into the environment
)set expose add con DenavitHartenbergMatrix
-- Create Affine transformations (DH matrices) that transform
-- a given triangle into another given triangle
-- tri2tri(t1, t2) returns a DHMATRIX which transforms t1 to t2,
-- where t1 and t2 are the vertices of two triangles in 3-space.
tri2tri(t1: List Point DoubleFloat, t2: List Point DoubleFloat): _
 DHMATRIX(DoubleFloat) ==
  n1 := triangleNormal(t1)
  n2 := triangleNormal(t2)
  tet2tet(concat(t1, n1), concat(t2, n2))
-- tet2tet(t1, t2) returns a DHMATRIX which transforms t1 to t2,
-- where t1 and t2 are the vertices of two tetrahedrons in 3-space.
tet2tet(t1: List Point DoubleFloat, t2: List Point DoubleFloat): _
 DHMATRIX(DoubleFloat) ==
  m1 := makeColumnMatrix t1
  m2 := makeColumnMatrix t2
  m2 * inverse(m1)
-- put the vertices of a tetrahedron into matrix form
makeColumnMatrix(t) ==
  m := new(4,4,0)$DHMATRIX(DoubleFloat)
  for x in t for i in 1..repeat
    for j in 1..3 repeat
      m(j,i) := x.j
    m(4,i) := 1
  m
-- return a vector normal to the given triangle, whose length
-- is the square root of the area of the triangle
triangleNormal(t) ==

  a := triangleArea t
  p1 := t.2 - t.1
  p2 := t.3 - t.2
  c := cross(p1, p2)
  len := length(c)
  len = 0 => error "degenerate triangle!"
  c := (1/len)*c
  t.1 + sqrt(a) * c
-- compute the are of a triangle using Heron's formula
triangleArea t ==
  a := length(t.2 - t.1)
  b := length(t.3 - t.2)
  c := length(t.1 - t.3)
  s := (a+b+c)/2
  sqrt(s*(s-a)*(s-b)*(s-c))
-- set up the coordinates of the corners of a tetrahedron
x1:DoubleFloat := sqrt(2.0@DoubleFloat/3.0@DoubleFloat)
x2:DoubleFloat := sqrt(3.0@DoubleFloat)/6
p1 := point [-0.5@DoubleFloat, -x2, 0.0@DoubleFloat]
p2 := point [0.5@DoubleFloat, -x2, 0.0@DoubleFloat]
p3 := point [0.0@DoubleFloat, 2*x2, 0.0@DoubleFloat]
p4 := point [0.0@DoubleFloat, 0.0@DoubleFloat, x1]
-- The base of the tetrahedron
baseTriangle  := [p2, p1, p3]
-- The "middle triangle" inscribed in the base of the tetrahedon
mt  := [0.5@DoubleFloat*(p2+p1), 0.5@DoubleFloat*(p1+p3), 0.5@DoubleFloat*(p3+p2)]
-- The bases of the triangles of the subdivided tetrahedon
bt1 := [mt.1, p1, mt.2]
bt2 := [p2, mt.1, mt.3]
bt3 := [mt.2, p3, mt.3]
bt4 := [0.5@DoubleFloat*(p2+p4), 0.5@DoubleFloat*(p1+p4), 0.5@DoubleFloat*(p3+p4)]
-- create the transformations which bring the base of the tetrahedron
-- to the bases of the subdivided tetrahedron
tt1 := tri2tri(baseTriangle, bt1)
tt2 := tri2tri(baseTriangle, bt2)
tt3 := tri2tri(baseTriangle, bt3)
tt4 := tri2tri(baseTriangle, bt4)
-- Draw a Sierpinsky tetrahedron with n levels of recursive subdivision
drawPyramid(n) ==
  s := create3Space()$ThreeSpace DoubleFloat
  dh := rotatex(0.0@DoubleFloat)
  drawPyramidInner(s, n, dh)

makeViewport3D(s, "Sierpinsky Tetrahedron")$VIEW3D
-- Recursively draw a Sierpinsky tetrahedron
drawPyramidInner(s, n, dh) ==
n = 0 => makeTetrahedron(s, dh, n)
-- draw the 4 recursive pyramids
drawPyramidInner(s, n-1, dh * tt1)
drawPyramidInner(s, n-1, dh * tt2)
drawPyramidInner(s, n-1, dh * tt3)
drawPyramidInner(s, n-1, dh * tt4)
-- draw a tetrahedron into the given space, with the given color,
-- transforming it by the given DH matrix
makeTetrahedron(sp, dh, color) ==
w1 := dh*p1
w2 := dh*p2
w3 := dh*p3
w4 := dh*p4
polygon(sp, [w1, w2, w4])
polygon(sp, [w1, w3, w4])
polygon(sp, [w2, w3, w4])
void()
drawPyramid 4
-- Antoine's Rings
-- Draw Antoine's Necklace
-- Thanks to Matt Grayson (formerly at IBM's T.J Watson Research Center)
-- for the idea.
-- Bring DH matrices into the environment
)set expose add con DenavitHartenbergMatrix
-- The current transformation for drawing a sub-ring
torusRot: DHMATRIX(DoubleFloat)
-- Draw Antoine's Rings with n levels of recursive subdivision.
-- The number of subrings is 10**n.
drawRings(n) ==
s := create3Space()$ThreeSpace DoubleFloat
-- create an identity transformation
dh:DHMATRIX(DoubleFloat) := identity()
drawRingsInner(s, n, dh)
makeViewport3D(s, "Antoine's Necklace")
-- Recursively draw Antoine's Necklace.
drawRingsInner(s, n, dh) ==
n = 0 =>
  drawRing(s, dh)

  void()
t := 0.0@DoubleFloat             -- the current angle around the ring
p := 0.0@DoubleFloat             -- the angle of the subring from the plane
tr := 1.0@DoubleFloat            -- the amount to translate the subring
inc := 0.1@DoubleFloat           -- translation increment
-- subdivide the ring into 10 linked rings
for i in 1..10 repeat
  tr := tr + inc
  inc := -inc
  dh' := dh * rotatez(t) * translate(tr, 0.0@DoubleFloat, 0.0@DoubleFloat) *
         rotatey(p) * scale(0.35@DoubleFloat, 0.48@DoubleFloat, 0.4@DoubleFloat)
  drawRingsInner(s, n-1, dh')
  t := t + 36.0@DoubleFloat
  p := p + 90.0@DoubleFloat
void()
-- draw a single ring into the given subspace, transformed by the given
-- DHMATRIX.
drawRing(s, dh) ==
free torusRot
torusRot := dh
makeObject(torus, 0..2*%pi, 0..2*%pi, var1Steps == 6, space == s,
           var2Steps == 15)
-- Parameterization of a torus, transformed by the DHMATRIX in torusRot.
torus(u ,v) ==
cu := cos(u)/6
torusRot * point [(1+cu)*cos(v), (1+cu)*sin(v), (sin u)/6]
drawRings 2
-- Scherk's Minimal surface
-- Defined by:
--    exp(z) * cos(x) = cos(y)
-- See: A comprehensive Introduction to Differential Geometry, Vol. 3,
--  by Michael Spivak, Publish Or Persih, Berkeley, 1979, pp 249-252.
-- Off set for a single piece of Scherk's minimal surface
(xOffset, yOffset):DoubleFloat
-- DrawScherk's minimal surface on an m by n patch.
drawScherk(m,n) ==
free xOffset, yOffset
space := create3Space()$ThreeSpace(DoubleFloat)
for i in 0..m-1 repeat
  xOffset := i*%pi
  for j in 0 .. n-1 repeat
    rem(i+j, 2) = 0 => 'iter
    yOffset := j*%pi
    drawOneScherk(space)

makeViewport3D(space, "Scherk's Minimal Surface")
-- The four patches which make up a single piece of Scherk's minimal surface.
scherk1(u,v) ==
x := cos(u)/exp(v)
point [xOffset + acos(x), yOffset + u, v, abs(v)]
scherk2(u,v) ==
x := cos(u)/exp(v)
point [xOffset - acos(x), yOffset + u, v, abs(v)]
scherk3(u,v) ==
x := exp(v) * cos(u)
point [xOffset + u, yOffset + acos(x), v, abs(v)]
scherk4(u,v) ==
x := exp(v) * cos(u)
point [xOffset + u, yOffset - acos(x), v, abs(v)]

-- We draw the surface by breaking it into 4 patches, and drawing them
-- into a single space.
drawOneScherk(s) ==
makeObject(scherk1, -%pi/2..%pi/2, 0..%pi/2,  space == s, _
           var1Steps == 28, var2Steps == 28)
makeObject(scherk2, -%pi/2..%pi/2, 0..%pi/2,  space == s, _
           var1Steps == 28, var2Steps == 28)
makeObject(scherk3, -%pi/2..%pi/2, -%pi/2..0, space == s, _
           var1Steps == 28, var2Steps == 28)
makeObject(scherk4, -%pi/2..%pi/2, -%pi/2..0, space == s, _
           var1Steps == 28, var2Steps == 28)
void()
drawScherk(3,3)
-- Ribbon Plot of [x**i for i in 1..5]
drawRibbons(flist, xrange, yrange) ==
sp := createThreeSpace()
num := # flist
yVar := variable yrange
y0:Float    := lo segment yrange
width:Float := (hi segment yrange - y0)/num
for f in flist for color in 1..num repeat
makeObject(f, xrange, yVar = y0..y0+width,
 var2Steps == 1, colorFunction == (x,y) +-> color,
  space ==sp)
y0:= y0 + width
vp := makeViewport3D(sp, "Ribbons")
drawStyle(vp, "shade")

outlineRender(vp, "on")
showRegion(vp, "on")
vp

)set message test off
drawRibbons([x**i for i in 1..5], x=-1..1, y=0..2)
)set message test on

———-

References^{[1] nothing}

9