Algol 68 Genie

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)
Deze website is gearchiveerd door de KB, nationale bibliotheek.

This website is archived by the National Library of the Netherlands.