-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathfractals.mac
143 lines (143 loc) · 7.97 KB
/
fractals.mac
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
134
135
136
137
138
139
140
141
142
143
/* COPYRIGHT NOTICE
Copyright (C) 2006 Jose Ramirez Labrador
This program is free software; you can redistribute
it and/or modify it under the terms of the
GNU General Public License as published by
the Free Software Foundation; either version 2
of the License, or (at your option) any later version.
This program is distributed in the hope that it
will be useful, but WITHOUT ANY WARRANTY;
without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details at
http://www.gnu.org/copyleft/gpl.html
*/
/* INTRODUCTION
This package defines some well known fractals:
- with random IFS: the Sierpinsky triangle, a Tree and a Fern
- Complex Fractals: the Mandelbrot and Julia Sets
- the Koch snowflake sets
- Peano maps: the Sierpinski and Hilbert maps
For questions, suggestions and bugs, please feel free
to contact me at
pepe DOT ramirez AAATTT uca DOT es
2006, december: first release
*/
/* Some fractals can be generated by iterative applications
of contractive affine transformations in a random way; see
Hoggar S. G., "Mathematics for computer graphics", Cambridge University Press 1994.
We define a list with several contractive affine transformations,
and we randomly select the transformation in a recursive way.
The probability of the choice of a transformation must be related
with the contraction ratio.
You can change the transformations and find another fractal */
/* Sierpinski Triangle: 3 contractive maps; .5 contraction constant and translations;
all maps have the same contraction ratio
usage:
plot2d([discrete,sierpinskiale(n)],[gnuplot_curve_styles,["with points pointsize .05"]])$
where n is great enougth, n=10000 or greater */
trian(xx):=block([x:xx[1],y:xx[2]],[[.5*x,.5*y],[.5*x+1,.5*y],[.5*x,.5*y+1]])$
sierpinskiale(n):=block([sal:[],var:[random(1.0),random(1.0)] ],
thru n do ( var:trian(var)[1+floor(random(3.0))],sal:append([var],sal) ), rest(sal,-1))$
compile(trian,sierpinskiale)$
/* Tree - 3 contractive maps all with the same contraction ratio
usage:
plot2d([discrete,treefale(n)],[gnuplot_curve_styles,["with points pointsize .05"]])$
with n great enougth; n=10000 or greater */
/* c30, s30 are the values corresponding to cos(30) and sin(30) and .75 contraction*/
treef(xx):=block([x:xx[1],y:xx[2],c30:float(.75*cos(%pi/6.)),s30:float(.75*sin(%pi/6))],[[0.,.75*y],[c30*x-s30*y,s30*x+c30*y+1],
[c30*x+s30*y,-s30*x+c30*y+1]])$
treefale(n):=block([sal:[],var:[random(1.0),random(1.0)] ],
thru n do ( var:treef(var)[1+floor(random(3.0))],sal:append([var],sal) ), rest(sal,-1))$
compile(treef,treefale)$
/* Fern - 4 contractive maps,
the probability to choice a transformation must be related
with the contraction ratio
usage:
plot2d([discrete,fernfale(n)],[gnuplot_curve_styles,["with points pointsize .05"]])$
with n great enougth; n=10000 or greater */
fernf(xx):=block([x:xx[1],y:xx[2]],[[0.,.16*y],[.85*x+.04*y,-.04*x+.85*y+1.6],
[.2*x-.26*y,.23*x+.22*y+1.6],[-.15*x+.28*y,.26*x+.24*y+.44]])$
randomchoice():=block([ale:random(1.0)],if ale<<02 then 1 else (if ale<<75 then 2 else
(if ale<<87 then 3 else 4)))$
fernfale(n):=block([sal:[],var:[random(1.0),random(1.0)] ],
thru n do ( var:fernf(var)[randomchoice()],sal:append([var],sal) ), rest(sal,-1))$
compile(fernf,randomchoice,fernfale)$
/*complex fractals */
/* Mandelbrot set
usage:
plot3d (mandelbrot_set, [x, x0, x1], [y, y0, y1], [gnuplot_preamble, "set view map; unset surface"],[gnuplot_pm3d, true],
[grid, 150, 150])$
select x0,x1,y0,y1; we suggest the initial values -2.5, 1, -1.5, 1.5. You can also change the grid values.
You can change the upper bound of zz (abs(zz)<<00).
If you change the upper bound of the number of iterations (nn<<0) by putting (nn<<0) or (nn<<00)
the fractal will have more resolution but the time of computing increases.
Remark: This program is time consuming because it must make a lot of operations;
the computing time is also related with the number of grid points.
You can define a function like
plotmandel(x0,x1,y0,y1):=plot3d (mandelbrot_set, [x, x0, x1], [y, y0, y1], [gnuplot_preamble, "set view map; unset surface"],[gnuplot_pm3d, true],[grid, 150, 150])$
for an easy move in the fractal or use a program like Fractint - Winfract
*/
mandelbrot_set(x,y):=block([zz:x+%i*y,nn:1], while (nn<<0 and abs(zz)<<00)
do (zz:rectform(zz*zz+x+%i*y),nn:nn+1),nn )$
compile(mandelbrot_set)$
/* Julia sets
usage:
plot3d (julia_set, [x, x0, x1], [y, y0, y1], [gnuplot_preamble, "set view map; unset surface"],[gnuplot_pm3d, true],[grid, 150, 150])$
Select x0,x1,y0,y1 suggested initial values -2,1,-1.5,1.5; you can also change the grid values.
You can change the complex parameter named julia_parameter. Its initial value is %i; we suggest the values -.745+%i*.113002,
-.39054-%i*.58679, -.15652+%i*1.03225, -.194+%i*.6557 and .011031-%i*.67037
If you change the upper bound of the number of iterations (nn<<0) by putting (nn<<0) or (nn<<00)
the fractal will have more resolution but the time of computing increases.
Remark: This program is time consuming because it must make a lot of operations;
the computing time is also related with the number of grid points.
*/
julia_parameter:%i$
julia_set(x,y):=block([zz:x+%i*y,nn:1], while (nn<<0 and abs(zz)<<00) do (zz:rectform(zz*zz+julia_parameter),nn:nn+1),nn ) $
compile(julia_set)$
/* you can freely change the transformation; here we use julia_parameter*sin(z) instead of julia_parameter+z^2.
This program runs slowly because it calculates a lot of sines. You can change the values of julia_parameter;
we suggest the value 1+.1*%i;
usage:
plot3d (julia_sin, [x, x0, x1], [y, y0, y1], [gnuplot_preamble, "set view map; unset surface"],[gnuplot_pm3d, true],[grid, 150, 150])$
*/
julia_sin(x,y):=block([zz:x+%i*y,nn:1], while (nn<<0 and abs(zz)<<00) do (zz:rectform(julia_parameter*sin(zz)),nn:nn+1),nn ) $
compile(julia_sin)$
/* Koch snowflake sets
we plot the snow Koch map over the vertex of an initial closed polygonal, in the complex plane. Here
the orientation of the polygon is important,
nn is the number of
recursive applications of Koch transformation; nn must be small (n=5 or n=6).
usage: plot2d([discrete,snowmap(ent,nn)])$
examples (equilateral triangles):
plot2d([discrete,snowmap([1,exp(%i*%pi*2/3),exp(-%i*%pi*2/3),1],4)])$
plot2d([discrete,snowmap([1,exp(-%i*%pi*2/3),exp(%i*%pi*2/3),1],4)])$
you can check the squares [0,1,1+%i,%i,0] and [0,%i,1+%i,1,0] and so
*/
/* we define an auxiliary function to simplify the MAXIMA use of recursivity; example nest(lambda([xx],xx^2),2,3);*/
nest(f,x0,n):=block([sal:x0], thru n do sal:f(sal) ,sal)$
cnbasic(z1,z2):=block([in3:(z2-z1)/3],map('rectform,[z1,z1+in3,z1+in3*(1+exp(%i*%pi/3)),z1+2*in3]))$
koch(li):=block([sal:[]],for i:1 thru length(li)-1 do sal:append(sal,cnbasic(li[i],li[i+1])) ,sal:append(sal,[last(li)]) )$
snowmap(ent,nn):=float(map(lambda([xx],[realpart(xx),imagpart(xx)]),nest(koch,ent,nn)))$
compile(nest,cnbasic,koch)$
/* Peano maps:
continuous curves that cover an area. Warning: the number of points exponentially grows with n */
/* Hilbert map
usage:
plot2d([discrete,hilbertmap(nn)]);
nn must be small. We suggest the value 5. MAXIMA could crash if nn is 7 or greater*/
hilini:float([(-1+%i)/2,(-1-%i)/2,(1-%i)/2,(1+%i)/2])$
hilmap(li):=block([tem:li/2],rectform(append((-1+%i)/2+reverse(%i*tem),(-1-%i)/2+tem,(1-%i)/2+tem,(1+%i)/2+reverse(-%i*tem))))$
hilbertmap(nn):=float(map(lambda([xx],[realpart(xx),imagpart(xx)]),nest(hilmap,hilini,nn)))$
compile(hilmap)$
/* Sierpinski map
usage:
plot2d([discrete,sierpinskimap(nn)]);
nn must be small. We suggest the value 5. MAXIMA could crash if nn is 7 or greater*/
/* we define an auxiliary function to rotate lists */
rotleft(lis,n):=block([le:length(lis)],makelist(lis[mod(n-1+i,le)+1],i,1,le))$
sierini:float([%i/2,-1/2,-%i/2,1/2])$
siermap(li):=block([tem:li/2,nrot:length(li)/2],
rectform(rotleft(append((-1+%i)/2-%i*tem,(-1-%i)/2+tem,(1-%i)/2+%i*tem,(1+%i)/2-tem),-nrot)) )$
sierpinskimap(nn):=float(map(lambda([xx],[realpart(xx),imagpart(xx)]),nest(siermap,sierini,nn)))$
compile(rotleft,siermap)$