\documentclass{article}
\usepackage{axiom}
\begin{document}
\title{\$SPAD/src/input images8a.input}
\author{The Axiom Team}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{License}
<<license>>=
--Copyright The Numerical Algorithms Group Limited 1994.
@
<<*>>=
<<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
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}