mandelbrot-plotutils.a68
1 COMMENT
2
3 @section Synopsis
4
5 Plot a part of the Mandelbrot set.
6
7 Plot (part of) the Mandelbrot set, that are points z[0] in the complex plane
8 for which the sequence z[n+1] := z[n] ** 2 + z[0] (n >= 0) is bounded.
9
10 This program requires a68g to be built with the GNU plotutils library.
11
12 COMMENT
13
14 PRAGMAT need plotutils PRAGMAT
15
16 INT pix = 300, max iter = 256, REAL zoom = 0.33 / pix;
17 [-pix : pix, -pix : pix] INT plane;
18 COMPL ctr = 0.05 I 0.75 # center of set #;
19
20 # Compute the length of an orbit. #
21
22 PROC iterate = (COMPL z0) INT:
23 BEGIN COMPL z := z0, INT iter := 1;
24 WHILE (iter +:= 1) < max iter # not converged # AND ABS z < 2 # not diverged #
25 DO z := z * z + z0
26 OD;
27 iter
28 END;
29
30 # Compute set and find maximum orbit length. #
31
32 INT max col := 0;
33 FOR x FROM -pix TO pix
34 DO FOR y FROM -pix TO pix
35 DO COMPL z0 = ctr + (x * zoom) I (y * zoom);
36 IF (plane [x, y] := iterate (z0)) < max iter
37 THEN (plane [x, y] > max col | max col := plane [x, y])
38 FI
39 OD
40 OD;
41
42 # Make a plot. #
43
44 FILE plot;
45 INT num pix = 2 * pix + 1;
46 open (plot, "mandelbrot.png", stand draw channel);
47 make device (plot, "png", whole (num pix, 0) + "x" + whole (num pix, 0));
48 FOR x FROM -pix TO pix
49 DO FOR y FROM -pix TO pix
50 DO INT col = (plane [x, y] > max col | max col | plane [x, y]);
51 REAL c = sqrt (1- col / max col); # sqrt to enhance contrast #
52 draw colour (plot, c, c, c);
53 draw point (plot, (x + pix) / (num pix - 1), (y + pix) / (num pix - 1))
54 OD
55 OD;
56 close (plot)
