1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
|
\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}
|
