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