share/fractals/fractals.mac

   1 /*               COPYRIGHT NOTICE
   2
   3 Copyright (C) 2006 Jose Ramirez Labrador
   4
   5 This program is free software; you can redistribute
   6 it and/or modify it under the terms of the
   7 GNU General Public License as published by
   8 the Free Software Foundation; either version 2
   9 of the License, or (at your option) any later version.
  10
  11 This program is distributed in the hope that it
  12 will be useful, but WITHOUT ANY WARRANTY;
  13 without even the implied warranty of MERCHANTABILITY
  14 or FITNESS FOR A PARTICULAR PURPOSE. See the
  15 GNU General Public License for more details at
  16 http://www.gnu.org/copyleft/gpl.html
  17 */
  18
  19
  20
  21
  22 /*             INTRODUCTION
  23
  24 This package defines some well known fractals:
  25
  26 - with random IFS: the Sierpinsky triangle, a Tree and a Fern
  27 - Complex Fractals: the  Mandelbrot and Julia Sets
  28 - the Koch snowflake sets
  29 - Peano maps: the Sierpinski and Hilbert maps
  30
  31
  32 For questions, suggestions and  bugs, please feel free
  33 to contact me at
  34
  35 pepe DOT ramirez AAATTT uca DOT es
  36
  37 2006, december: first release
  38
  39 */
  40
  41
  42 /* Some fractals can be generated by iterative applications
  43 of contractive affine transformations in a random way; see
  44
  45 Hoggar S. G., "Mathematics for computer graphics", Cambridge University Press 1994.
  46
  47 We define a list with several contractive affine transformations,
  48 and  we randomly select the transformation in a recursive way.
  49 The probability of the choice of a transformation must be related
  50 with the contraction ratio.
  51
  52 You can change the transformations and find another fractal   */
  53
  54
  55 /* Sierpinski Triangle: 3 contractive maps; .5 contraction constant and translations;
  56 all maps have the same contraction ratio
  57 usage:
  58
  59 plot2d([discrete,sierpinskiale(n)],[gnuplot_curve_styles,["with points pointsize .05"]])$
  60
  61 where  n is great enough, n=10000 or greater */
  62
  63 trian(xx):=block([x:xx[1],y:xx[2]],[[.5*x,.5*y],[.5*x+1,.5*y],[.5*x,.5*y+1]])$
  64 sierpinskiale(n):=block([sal:[],var:[random(1.0),random(1.0)] ],
  65 thru n do ( var:trian(var)[1+floor(random(3.0))],sal:append([var],sal) ), rest(sal,-1))$
  66 compile(trian,sierpinskiale)$
  67
  68 /* Tree - 3 contractive maps all with the same contraction ratio
  69 usage:
  70
  71 plot2d([discrete,treefale(n)],[gnuplot_curve_styles,["with points pointsize .05"]])$
  72
  73 with n great enough; n=10000 or greater */
  74
  75
  76 /* c30, s30 are the values corresponding to cos(30) and sin(30) and .75 contraction*/
  77
  78 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],
  79 [c30*x+s30*y,-s30*x+c30*y+1]])$
  80
  81 treefale(n):=block([sal:[],var:[random(1.0),random(1.0)] ],
  82 thru n do ( var:treef(var)[1+floor(random(3.0))],sal:append([var],sal) ), rest(sal,-1))$
  83 compile(treef,treefale)$
  84
  85
  86
  87 /* Fern - 4 contractive maps,
  88 the probability to choice a transformation must be related
  89 with the contraction ratio
  90
  91 usage:
  92
  93 plot2d([discrete,fernfale(n)],[gnuplot_curve_styles,["with points pointsize .05"]])$
  94
  95 with n great enough; n=10000 or greater */
  96
  97
  98 fernf(xx):=block([x:xx[1],y:xx[2]],[[0.,.16*y],[.85*x+.04*y,-.04*x+.85*y+1.6],
  99 [.2*x-.26*y,.23*x+.22*y+1.6],[-.15*x+.28*y,.26*x+.24*y+.44]])$
 100 randomchoice():=block([ale:random(1.0)],if ale<.02 then 1 else (if ale<.75 then 2 else
 101 (if ale<.87 then 3 else 4)))$
 102 fernfale(n):=block([sal:[],var:[random(1.0),random(1.0)] ],
 103 thru n do ( var:fernf(var)[randomchoice()],sal:append([var],sal) ), rest(sal,-1))$
 104 compile(fernf,randomchoice,fernfale)$
 105
 106
 107 /*complex fractals */
 108
 109 /* Mandelbrot set
 110
 111 usage:
 112 plot3d (mandelbrot_set, [x, x0, x1], [y, y0, y1], [gnuplot_preamble, "set view map; unset surface"],[gnuplot_pm3d, true],
 113 [grid, 150, 150])$
 114
 115  select x0,x1,y0,y1; we suggest the initial values  -2.5, 1, -1.5, 1.5.  You can also change the grid values.
 116
 117 You can change the upper bound of zz (abs(zz)<100).
 118 If you change the upper bound of the number of iterations (nn<30) by putting (nn<50) or (nn<100)
 119 the fractal will have more resolution but the time of computing increases.
 120
 121 Remark: This program is time consuming because it must make a lot of operations;
 122 the computing time is also related with the number of grid points.
 123
 124 You can define a function like
 125 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])$
 126 for an easy move in the fractal or use a program like Fractint - Winfract
 127 */
 128
 129 mandelbrot_set(x,y):=block([zz:x+%i*y,nn:1], while (nn<30 and abs(zz)<100)
 130 do (zz:rectform(zz*zz+x+%i*y),nn:nn+1),nn )$
 131 compile(mandelbrot_set)$
 132
 133
 134 /* Julia sets
 135
 136 usage:
 137 plot3d (julia_set, [x, x0, x1], [y, y0, y1], [gnuplot_preamble, "set view map; unset surface"],[gnuplot_pm3d, true],[grid, 150, 150])$
 138
 139  Select x0,x1,y0,y1 suggested initial values  -2,1,-1.5,1.5; you can also change the grid values.
 140
 141 You can change the complex parameter named julia_parameter. Its initial value is %i; we  suggest the  values -.745+%i*.113002,
 142 -.39054-%i*.58679, -.15652+%i*1.03225, -.194+%i*.6557 and .011031-%i*.67037
 143
 144
 145 If you change the upper bound of the number of iterations (nn<30) by putting (nn<50) or (nn<100)
 146 the fractal will have more resolution but the time of computing increases.
 147
 148 Remark: This program is time consuming because it must make a lot of operations;
 149 the computing time is also related with the number of grid points.
 150
 151 */
 152
 153 julia_parameter:%i$
 154 julia_set(x,y):=block([zz:x+%i*y,nn:1], while (nn<30 and abs(zz)<100) do (zz:rectform(zz*zz+julia_parameter),nn:nn+1),nn ) $
 155 compile(julia_set)$
 156
 157 /* you can freely change the transformation; here  we use julia_parameter*sin(z) instead of julia_parameter+z^2.
 158 This program runs slowly  because it calculates a lot of sines. You can change the values of julia_parameter;
 159 we suggest the value 1+.1*%i;
 160
 161 usage:
 162 plot3d (julia_sin, [x, x0, x1], [y, y0, y1], [gnuplot_preamble, "set view map; unset surface"],[gnuplot_pm3d, true],[grid, 150, 150])$
 163
 164 */
 165
 166 julia_sin(x,y):=block([zz:x+%i*y,nn:1], while (nn<30 and abs(zz)<100) do (zz:rectform(julia_parameter*sin(zz)),nn:nn+1),nn ) $
 167 compile(julia_sin)$
 168
 169
 170 /* Koch snowflake sets
 171 we plot the snow Koch map over the vertex of an initial closed  polygonal, in the complex plane. Here
 172 the orientation of the polygon  is important,
 173 nn is the number of
 174 recursive applications of Koch transformation; nn must be small (n=5 or n=6).
 175
 176 usage: plot2d([discrete,snowmap(ent,nn)])$
 177
 178 examples (equilateral triangles):
 179 plot2d([discrete,snowmap([1,exp(%i*%pi*2/3),exp(-%i*%pi*2/3),1],4)])$
 180
 181 plot2d([discrete,snowmap([1,exp(-%i*%pi*2/3),exp(%i*%pi*2/3),1],4)])$
 182
 183 you can check the squares [0,1,1+%i,%i,0] and [0,%i,1+%i,1,0] and so
 184 */
 185
 186
 187 /* we define an auxiliary function to simplify the MAXIMA use of recursivity; example nest(lambda([xx],xx^2),2,3);*/
 188
 189 nest(f,x0,n):=block([sal:x0], thru n do sal:f(sal) ,sal)$
 190
 191 cnbasic(z1,z2):=block([in3:(z2-z1)/3],map('rectform,[z1,z1+in3,z1+in3*(1+exp(%i*%pi/3)),z1+2*in3]))$
 192
 193 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)]) )$
 194
 195 snowmap(ent,nn):=float(map(lambda([xx],[realpart(xx),imagpart(xx)]),nest(koch,ent,nn)))$
 196 compile(nest,cnbasic,koch)$
 197
 198
 199
 200 /* Peano maps:
 201 continuous curves that cover an area. Warning: the number of points exponentially grows with n */
 202
 203 /* Hilbert map
 204 usage:
 205
 206 plot2d([discrete,hilbertmap(nn)]);
 207
 208  nn must be small. We suggest the value 5. MAXIMA could crash if nn is 7 or greater*/
 209
 210 hilini:float([(-1+%i)/2,(-1-%i)/2,(1-%i)/2,(1+%i)/2])$
 211 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))))$
 212 hilbertmap(nn):=float(map(lambda([xx],[realpart(xx),imagpart(xx)]),nest(hilmap,hilini,nn)))$
 213 compile(hilmap)$
 214
 215 /* Sierpinski map
 216 usage:
 217 plot2d([discrete,sierpinskimap(nn)]);
 218
 219  nn must be small. We suggest the value 5. MAXIMA could crash if nn is 7 or greater*/
 220
 221 /* we define an auxiliary function to rotate lists */
 222 rotleft(lis,n):=block([le:length(lis)],makelist(lis[mod(n-1+i,le)+1],i,1,le))$
 223
 224
 225 sierini:float([%i/2,-1/2,-%i/2,1/2])$
 226 siermap(li):=block([tem:li/2,nrot:length(li)/2],
 227 rectform(rotleft(append((-1+%i)/2-%i*tem,(-1-%i)/2+tem,(1-%i)/2+%i*tem,(1+%i)/2-tem),-nrot)) )$
 228
 229 sierpinskimap(nn):=float(map(lambda([xx],[realpart(xx),imagpart(xx)]),nest(siermap,sierini,nn)))$
 230 compile(rotleft,siermap)$