  ## 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 in the complex plane
8  for which the sequence z[n+1] := z[n] ** 2 + z (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)

```