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

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);
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:
• Csq for squaring a complex number
• Csize for determining 'the length' of a complex number (although the length would be the root of it)
The stopping criteria are either
• Z becoming too big (z > 2 or z-squared > 4)
• iteration count becoming too big
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:

• X ordinate
• Y coordinate
• threshold value for 'n'
• granularity value for a kind of zooming
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

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);
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

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);
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`     Long reals for more detail: Mandel04

Sooner or later you dive in to deep. You go down the mandelbrot alley too deep, exploring places where no man has gone before. A commandline like

`./mand03 -0.6735 -0.358 2000 600000000`
will zoom in to 1.16EE-08 precision and that is beyond the accuracy of a normal REAL number. You see that on screen when the smallest details are rectangles scattered over the place and around perimeters. The solution is to use LONGREALs for the calculations.

Below is the source for mand04.mod. This source is exactly 8 bytes longer than mand03.mod since in two places the word "REAL" is replaced by the word "LONGREAL". Look for yourself:

```MODULE mand04;

IMPORT	Args, Conv, Err, XYplane;

TYPE	Complex	= RECORD
real, imag	: LONGREAL
END;

VAR	z, z0, c		: Complex;
startR, startI		: LONGREAL;
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

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.LongRealVal (option);
Args.GetArg (2, option);
startI := Conv.LongRealVal (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);
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 mand04.
```
The results are astonishing. The speed penalty is marginal. The accuracy increase is incredible. The blurred image mentioned above is back with full detail. The pictures below were created with the commandlines
```jan@fluor:~/Oberon/mand\$ ./mand03 -0.6735 -0.3625 200 300000000
jan@fluor:~/Oberon/mand\$ ./mand04 -0.6735 -0.3625 200 300000000
```  `obc -o exefile mand0x.mod`