$SPAD/src/input mountain.input

The Axiom Team December 3, 2016

Abstract

Contents

2

1 License

— license —

--Copyright The Numerical Algorithms Group Limited 1994.

———-

— * —

\getchunk{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)

———-

References^{[1] nothing}

6