Producing Mandelbrot plots

Back in '88 I read this article about producing Mandelbrot plots. I put it on hold for some weeks (let's say 1500 weeks) and now I finally got the inspiration to do it. I made some programs to produce the curly swirly plots. Using the XYplane module that comes with the obc compiler.


The algorithm

I will not bother you with what the set is and how to calculate it and all the other things you need not know. If you do insist and want to know, just consult wikipedia. what it boils down to is that complex numbers are presented in the complex plane as if they were cartesian coordinates x and y:

z = x + iy
And then we start an iteration on that complex number z according to
Z1 = Z0^2 + C
and keep on doing this until Z1 becomes a run-away (i.e. only gets bigger after each iteration). Some numbers are immediate runners. Some take 'n' iterations to become a runner. The mandelbrot set is about that 'n'.

First attempt : mand01

Below is the sourcecode for mand01.mod. It runs, but to enter new starting coordinates you need to edit the source and recompile. Good enough for a proof of concept.

MODULE mand01;

IMPORT	XYplane;

TYPE	Complex	= RECORD
		    real, imag	: REAL
		  END;

VAR	z, z0, c		: Complex;
	startR, startI		: REAL;
	n, x, y			: INTEGER;


PROCEDURE Csq (VAR  z : Complex);

VAR	res	: Complex;

BEGIN
  res.real := z.real * z.real - z.imag * z.imag;
  res.imag := 2 * z.real * z.imag;
  z := res
END Csq;


PROCEDURE Cadd (VAR z : Complex;  y : Complex);

VAR	res	: Complex;

BEGIN
  res.real := y.real + z.real;
  res.imag := y.imag + z.imag;
  z := res
END Cadd;


PROCEDURE Csize (VAR z : Complex) : INTEGER;

VAR	len	: INTEGER;

BEGIN
  len := ENTIER (z.real * z.real + z.imag * z.imag);
  RETURN len
END Csize;


BEGIN
  XYplane.Open;
  z0.real := 0.0;		z0.imag := 0.0;
  startR := -0.5;		startI := -1.25;
  FOR  x := 0 TO 639 DO
    c.real := startR + x / 200;
    FOR  y := 0  TO  479  DO
      c.imag := startI + y / 200;
      n := 0;
      z := z0;
      REPEAT
        INC (n);
	Csq (z);
	Cadd (z, c)
      UNTIL  (Csize (z) > 4) OR (n > 100);
      IF  n < 100  THEN  XYplane.Dot (x, y, 1)  ELSE  XYplane.Dot (x, y, 0)  END
    END
  END;
  REPEAT  UNTIL XYplane.Key () > ' '
END mand01.
   
I define a complex number as a record consisting of 2 REAL's. Now it becomes easy to access the real and imaginary parts of the complex number. I defined three complex number routines for betetr program structure: The stopping criteria are either If the iteration threshold is set to 100 and the 'n' is 101, then we consider this number to be in the mandelbrot set for threshold 100. The threshold is a kind of quality factor. Below are some examples of what can produce. Remember: each parameter change requires a recompile.
jan@fluor:~/Oberon/mand$ ./mand
jan@fluor:~/Oberon/mand$ jed mand01.mod 
jan@fluor:~/Oberon/mand$ obc -o man mand01.mod
jan@fluor:~/Oberon/mand$ ./man
jan@fluor:~/Oberon/mand$ jed mand01.mod 
jan@fluor:~/Oberon/mand$ obc -o man mand01.mod
jan@fluor:~/Oberon/mand$ ./man
jan@fluor:~/Oberon/mand$ jed mand01.mod 
jan@fluor:~/Oberon/mand$ obc -o man mand01.mod
jan@fluor:~/Oberon/mand$ ./man
jan@fluor:~/Oberon/mand$ 
   
Not very convenient.

Version 2, with command line support : mand02.mod

The first version made clear the algorithm worked. But having to recompile for each attempt (and finding nice patterns is purely trial and error) I made mand02 that accepts parameters on the command line:

The commandline is now something like:
./mand02 -0.73 0.32 100 10000
which sets the constant C to (-0.73 + 0.32 i) with maximal 100 iterations before we call it a day and a step size of 1/10000.

MODULE mand02;

IMPORT	Args, Conv, Err, XYplane;

TYPE	Complex	= RECORD
		    real, imag	: REAL
		  END;

VAR	z, z0, c		: Complex;
	startR, startI		: REAL;
	count, step,
	n, x, y			: INTEGER;
	option			: ARRAY 20 OF CHAR;


PROCEDURE Csq (VAR  z : Complex);

VAR	res	: Complex;

BEGIN
  res.real := z.real * z.real - z.imag * z.imag;
  res.imag := 2 * z.real * z.imag;
  z := res
END Csq;


PROCEDURE Cadd (VAR z : Complex;  y : Complex);

VAR	res	: Complex;

BEGIN
  res.real := y.real + z.real;
  res.imag := y.imag + z.imag;
  z := res
END Cadd;


PROCEDURE Csize (VAR z : Complex) : INTEGER;

VAR	len	: INTEGER;

BEGIN
  len := ENTIER (z.real * z.real + z.imag * z.imag);
  RETURN len
END Csize;


PROCEDURE Init;

BEGIN
  IF  Args.argc < 4  THEN
    Err.String ("Mandel X0 Y0 count step");	Err.Ln;
    startR := -0.5;
    startI := -1.25;
    count := 100;
    step := 200
  ELSE
    Args.GetArg (1, option);
    startR := Conv.RealVal (option);
    Args.GetArg (2, option);
    startI := Conv.RealVal (option);
    Args.GetArg (3, option);
    count := Conv.IntVal (option);
    Args.GetArg (4, option);
    step := Conv.IntVal (option)
  END;
  z0.real := 0.0;		z0.imag := 0.0;
  XYplane.Open
END Init;


BEGIN
  Init;
  FOR  x := 0 TO 639 DO
    c.real := startR + x / step;
    FOR  y := 0  TO  479  DO
      c.imag := startI + y / step;
      n := 0;
      z := z0;
      REPEAT
        INC (n);
	Csq (z);
	Cadd (z, c)
      UNTIL  (Csize (z) > 4) OR (n > count);
      IF  n > count  THEN  XYplane.Dot (x, y, 1)  ELSE  XYplane.Dot (x, y, 0)  END
    END
  END;
  REPEAT  UNTIL XYplane.Key () > ' '
END mand02.
   
Only minor changes in the algorithm (step and count instead of fixed values) and an Init section added. See it run, for the respective command lines:
      ./mand02a -0.72 0.3125 100 10000
      ./mand02a -0.71 0.3125 100 10000
      ./mand02a -0.71 0.3125 100 1000
      ./mand02 -0.73 0.3125 100 10000
      ./mand02 -0.73 0.33125 100 10000
      ./mand02 -0.73 0.33125 1000 10000
mand02 produces inverted colours from mando2a. For some images white on black looks better than the other way round. Remember, this is a technical website.

Look at the command lines and the pictures produced. Improving the quality factor completely 'ruins' the image.... Changing the R/I values has a major influence.

Version 3: center the starting point (mand03.mod)

Mand01 and -02 had the initial coordinates set to the lower left of the screen and all successive new points went to the right and to the top. This version put's the complex number in the center of the screen and moves around it. Now it is easier to see which effect parameter changes have on the image. Not many changes; just the looping constants have another name (X0 and Y0 instead of x and y) and some corrections are required later on. Here is the source:

MODULE mand03;

IMPORT	Args, Conv, Err, XYplane;

TYPE	Complex	= RECORD
		    real, imag	: REAL
		  END;

VAR	z, z0, c		: Complex;
	startR, startI		: REAL;
	count, step, X0, Y0,
	n, x, y			: INTEGER;
	option			: ARRAY 20 OF CHAR;


PROCEDURE Csq (VAR  z : Complex);

VAR	res	: Complex;

BEGIN
  res.real := z.real * z.real - z.imag * z.imag;
  res.imag := 2 * z.real * z.imag;
  z := res
END Csq;


PROCEDURE Cadd (VAR z : Complex;  y : Complex);

VAR	res	: Complex;

BEGIN
  res.real := y.real + z.real;
  res.imag := y.imag + z.imag;
  z := res
END Cadd;


PROCEDURE Csize (VAR z : Complex) : INTEGER;

VAR	len	: INTEGER;

BEGIN
  len := ENTIER (z.real * z.real + z.imag * z.imag);
  RETURN len
END Csize;


PROCEDURE Init;

BEGIN
  IF  Args.argc < 4  THEN
    Err.String ("Mandel X0 Y0 count step");	Err.Ln;
    startR := -0.71;
    startI := -0.36;
    count := 100;
    step := 40000
  ELSE
    Args.GetArg (1, option);
    startR := Conv.RealVal (option);
    Args.GetArg (2, option);
    startI := Conv.RealVal (option);
    Args.GetArg (3, option);
    count := Conv.IntVal (option);
    Args.GetArg (4, option);
    step := Conv.IntVal (option)
  END;
  z0.real := 0.0;		z0.imag := 0.0;
  XYplane.Open
END Init;


BEGIN
  Init;
  FOR  X0 := 0 TO 639 DO
    x := X0 - 320;
    c.real := startR + x / step;
    FOR  Y0 := 0  TO  479  DO
      y := Y0 - 240;
      c.imag := startI + y / step;
      n := 0;
      z := z0;
      REPEAT
        INC (n);
	Csq (z);
	Cadd (z, c)
      UNTIL  (Csize (z) > 4) OR (n > count);
      IF  n < count  THEN  XYplane.Dot (X0, Y0, 1)  ELSE  XYplane.Dot (X0, Y0, 0)  END
    END
  END;
  REPEAT  UNTIL XYplane.Key () > ' '
END mand03.
   
Compile it, as usual, with
obc -o mand03 mand03.mod
The pictures below were generated with
./mand03 -0.6735 -0.358 100 30000
and
./mand03 -0.6735 -0.358 200 30000

Downloads

You can download the images with right click etc. You can download the sources with cut and paste or via this link: mand03.mod There is a chance that your browser tries t play this file in your music player since 30 eyars ago MOD was a music format of Amiga. In that case download with Shift-Leftclick or with RightClick.

Page created 3 Nov 2018,