\documentclass{article}
\usepackage{axiom}
\begin{document}
\title{\$SPAD/src/input mountain.input}
\author{The Axiom Team}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{License}
<<license>>=
--Copyright The Numerical Algorithms Group Limited 1994.
@
<<*>>=
<<license>>

-- Draw a fractal mountain

)clear all

-- compile the functions
)set function compile on

-- Generate Gaussian random numbers
-- Algorithm by Richard Voss from "The Science of Fractal Images", pg. 77

-- function to convert a number into machine floating point
sf f == f::DFLOAT

Nrand := 4
Arand := 2**26 - 1
GaussAdd := sqrt(sf(3.0) * Nrand)
GaussFac := sf(2.0) * GaussAdd/((sf Nrand) * (sf Arand))

-- generate a random number
Gauss() ==
  sum := sf 0.0
  for i in 1..Nrand repeat
    sum := sum + random()$INT
  GaussFac * sum - GaussAdd

-- Generate fractal mountains.

-- Algorithms by Richard Voss from "The Science of Fractal Images", pg. 100

sfHalf  := sf 0.5
sfThree := sf 3.0
sfFour  := sf 4.0

f3(delta,x0,x1,x2) == (x0+x1+x2)/sfThree + delta*Gauss()
f4(delta,x0,x1,x2,x3) == (x0+x1+x2+x3)/sfFour + delta*Gauss()

-- perform midpoint subdivision
MidPointFM(maxLevel, sigma, H) ==
  N := 2**maxLevel
  delta := sigma
  arraySize := (N+1)
  X:IARRAY2(DFLOAT,0,0) := new(arraySize, arraySize, sf 0.0)
  setelt(X, 0, 0, delta*Gauss())
  setelt(X, 0, N, delta*Gauss())
  setelt(X, N, 0, delta*Gauss())
  setelt(X, N, N, delta*Gauss())
  D := N
  d := N quo 2
  for stage in 1..maxLevel repeat
    delta := delta*(sfHalf**(sfHalf*H))
    for x in d..(N-d) by D repeat
      for y in d..(N-d) by D repeat
        setelt(X, x, y, f4(delta, elt(X,x+d,y+d), elt(X,x+d,y-d),
                           elt(X, x-d, x+d), elt(X, x-d, y-d)))
    for x in 0..N by D repeat
      for y in 0..N by D repeat
        setelt(X, x, y, elt(X,x,y) + delta*Gauss())
    delta := delta*(sfHalf**(sfHalf*H))
    for x in d..(N-d) by D repeat
      setelt(X,x,0, f3(delta, elt(X,x+d,0), elt(X,x-d,0), elt(X,x,d)))
      setelt(X,x,N, f3(delta, elt(X,x+d,N), elt(X,x-d,N), elt(X,x,N-d)))
      setelt(X,0,x, f3(delta, elt(X,0,x+d), elt(X,0,x-d), elt(X,d,x)))
      setelt(X,N,x, f3(delta, elt(X,N,x+d), elt(X,N,x-d), elt(X,N-d,x)))
    for x in d..(N-d) by D repeat
      for y in D..(N-d) by D repeat
        setelt(X,x,y, f4(delta, elt(X,x,y+d), elt(X,x,y-d), 
                         elt(X,x+d,y), elt(X,x-d,y)))
    for x in D..(N-d) by D repeat
      for y in d..(N-d) by D repeat
        setelt(X,x,y, f4(delta, elt(X,x,y+d), elt(X,x,y-d), 
                         elt(X,x+d,y), elt(X,x-d,y)))
    for x in 0..N by D repeat
      for y in 0..N by D repeat
        setelt(X,x,y, elt(X,x,y) + delta*Gauss())
    for x in d..(N-d) by D repeat
      for y in d..(N-d) by D repeat
        setelt(X,x,y, elt(X,x,y) + delta*Gauss())
    D := D quo 2
    d := d quo 2
  X


sfZero := sf 0
Sigma := sf 7

-- function passed to the draw 
tableVal(x: DFLOAT, y:DFLOAT):DFLOAT ==
  free table, xIndex, yIndex, rowSize
  val := elt(table, xIndex, yIndex)
  xIndex := xIndex + 1
  if xIndex > rowSize then (xIndex := 0; yIndex := yIndex + 1)
  val < sfZero => sfZero
  val

-- draw a mountain with maxLevel subdivisions with Haussdorf dimension H
-- the number of subdivisions of the mountain is 2**maxLevel, so you 
-- probably should keep maxLevel <= 8.  Also 0 < H <= 1.  The closer
-- H is to one, the smoother the mountain will be.
drawMountain(maxLevel, H) ==
  free table, xIndex, yIndex, rowSize
  table := MidPointFM(maxLevel, Sigma, H)
  N := 2**maxLevel
  xIndex := 0
  yIndex := 0
  rowSize := N
  draw(tableVal, -20..20, -20..20,
    var1Steps == N, var2Steps == N, title == "Fractal Mountain")

drawMountain(3, sf 0.95)
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}