# To unbundle, run this file
echo mandel.c
sed 's/^X//' >mandel.c <<'!'
X#include <u.h>
X#include <libc.h>
X#include <draw.h>
X#include <event.h>
X#include <cursor.h>
X#include "mandel.h"

static int nsteps = NSTEPSMIN, colbrush, imodmins, imodmaxs;
static int nstepsx, nstepsy;
static Fp stepx, stepy;
static Complex	minpnt = { -1.5, 1 },
X		maxpnt = { 0.5, -1 };
static Complex mins = {MIN, MIN}, maxs = {MAX, MAX};
static Complex pixarray[MASRESX][MASRESY];
static Point pmin, pmax;
static Rectangle sweeprect;
static Image *background;
static Image *palette[NCOLS];
static Image *imagebuf = nil;

X/*
X * takes a number which is the position of the
X * byte for that primary color in the 4 bytes of the color int
X * colbrush and a number to sum to that color.
X */
int
addcolor(int colpos, int desvcol, int colbrush)
X{
X	int colormsk, thiscol;

X	colormsk = 0xFF;
X	colormsk <<= colpos * 8;

X	thiscol = colbrush & colormsk;
X	thiscol >>= colpos * 8;

X	thiscol += desvcol;

X	thiscol <<= colpos * 8;
X	thiscol &= colormsk;

X	colbrush &= ~colormsk;
X	colbrush |= thiscol;
X	return colbrush;
X}

X#ifdef SLOW
Complex
compmul(Complex *a, Complex *b)
X{
X	Fp are = a->re, aim = a->im, bre = b->re, bim = b->im;

X	return (Complex){are*bre - aim*bim, are*bim + aim*bre};
X}

Complex
compsum(Complex *a, Complex *b)
X{
X	return (Complex){a->re + b->re, a->im + b->im};
X}

XFp
module(Complex *a)
X{
X	Fp re = a->re, im = a->im;

X	return sqrt(re*re + im*im);
X//	return hypot(a->re, a->im);
X}

Complex
iterate(Complex a, int numsteps)
X{
X	Fp modz;
X	Complex z = a, z2;

X	while (numsteps-- > 0) {
X		z2 = compmul(&z, &z);
X		z = compsum(&z2, &a);
X		modz = module(&z);
X		if ( /* (int) */ modz > MAX)
X			return (Complex){MAX, MAX};
X		if (modz < MIN)
X			return (Complex){MIN, MIN};
X	}
X	return z;
X}
X#else			/* SLOW */

X#define compmul(a, b) (Complex){(a)->re*(b)->re - (a)->im*(b)->im,\
X				(a)->re*(b)->im + (a)->im*(b)->re}
X#define compsum(a, b) (Complex){(a)->re + (b)->re, (a)->im + (b)->im}
X// #define module(a) (Fp)hypot((a)->re, (a)->im)
X#define module(a) (Fp)sqrt((a)->re*(a)->re + (a)->im*(a)->im)
X#endif			/* SLOW */

static void
drawcol(int isteps, int x, Image *pal, Point *pixelp, Complex *pixpntp)
X{
X	int y, numsteps, dodraw;
X	Fp modz, incy = stepy;
X	Complex z, z2;
X	Complex *pixp = pixarray[x];

X	pixelp->x++;			// for x = 1
X	pixelp->y = pmin.y;
X	pixpntp->re += stepx;		// for x = 1
X	pixpntp->im = minpnt.im;
X	for (y = 1; y < nstepsy - 1; y++) {
X		pixelp->y++;		// for y = 1
X		pixpntp->im -= incy;	// for y = 1
X		// about half of all points pass these tests.
X		// the module(pixp)!=MIN test is costly and never true.
X		++pixp;			// for y = 1
X		if (pixp->re != MAX && pixp->im != MAX
X		    /* && module(pixp) != MIN */) {
X#ifdef SLOW
X			*pixp = iterate(pixelpnt, isteps);
X			dodraw = ((int)module(pixp) < MAX);
X#else
X			z = *pixpntp;
X			for (numsteps = isteps; numsteps > 0; numsteps--) {
X				z2 = compmul(&z, &z);
X				z = compsum(&z2, pixpntp);
X				modz = module(&z);
X				if (modz > MAX) {
X					*pixp = maxs;
X					dodraw = 0;
X					break;
X				}
X				if (modz < MIN) {
X					*pixp = mins;
X					dodraw = 1;
X					break;
X				}
X			}
X			if (numsteps <= 0) {
X				*pixp = z;
X				dodraw = ((int)modz < MAX);
X			}
X			// module(pixp)==modz
X#endif				/* SLOW */
X		} else
X			dodraw = ((int)module(pixp) < MAX);
X		if (dodraw)
X			draw(imagebuf, Rect(pixelp->x, pixelp->y, pixelp->x + 1,
X				pixelp->y + 1), pal, nil, ZP);
X	}
X	if (0 && x%200 == 0)
X		draw(screen, screen->r, imagebuf, nil, screen->r.min);
X}

void
redraw(void)
X{
X	int x, y, isteps;
X	Point pixel;
X	Complex pixelpnt;
X	Complex *colp;
X	Image *pal;

X	// tile pixarray with MIN
X	for (x = 0; x < MASRESX; x++) {
X		colp = pixarray[x];
X		for (y = MASRESY; y-- > 0; )
X			*colp++ = mins;
X	}

X	draw(screen, screen->r, background, nil, ZP);
X	if (imagebuf != nil)
X		freeimage(imagebuf);
X	if (screen->r.max.x - screen->r.min.x > MASRESX)
X		screen->r.max.x = screen->r.min.x + MASRESX;
X	if (screen->r.max.y - screen->r.min.y > MASRESY)
X		screen->r.max.y = screen->r.min.y + MASRESY;
X	imagebuf = allocimage(display, screen->r, CMAP8, 0, DBlack);
X	if (imagebuf == 0)
X		exits("out of image memory");
X	assert(imagebuf);

X	/* create two points on the corners of the square*/
X	pmin = Pt(imagebuf->r.min.x, imagebuf->r.min.y);
X	pmax = Pt(imagebuf->r.max.x, imagebuf->r.max.y);

X	nstepsx = pmax.x - pmin.x;
X	nstepsy = pmax.y - pmin.y;
X	stepx = (maxpnt.re - minpnt.re) / nstepsx;
X	stepy = (minpnt.im - maxpnt.im) / nstepsy;
X	imodmins = module(&mins);
X	imodmaxs = module(&maxs);

X	for (isteps = 11; isteps < nsteps + 11; isteps++) {
X		pal = palette[(isteps*19) % NCOLS];
X		pixel = pmin;
X		pixelpnt = minpnt;
X		for (x = 1; x < nstepsx - 1; x++)
X			drawcol(isteps, x, pal, &pixel, &pixelpnt);
X		draw(screen, screen->r, imagebuf, nil, screen->r.min);
X	}
X}

void
eresized(int new)
X{
X	if (new && getwindow(display, Refnone) < 0)
X		exits("getwindow failed");
X	redraw();
X}

void
main(int argc, char *argv[])
X{
X	int colmap = MINCOLOR;
X	int key, isteps;
X	Mouse mev;

X	USED(argc, argv);
X	if (initdraw(nil, nil, "mandel") < 0)
X		exits("initdraw failed");
X	/*
X	 * I allocate  the color palette to draw
X	 */
X	for (isteps = 10; isteps < NCOLS; isteps++) {
X		palette[isteps] = allocimage(display, Rect(0, 0, 1, 1),
X			CMAP8, 1, colmap);
X		if (palette[isteps] == 0)
X			exits("out of image memory");

X		colmap = addcolor(OFFGREEN, 1, colmap);
X	}

X	background = allocimage(display, Rect(0, 0, 1, 1), CMAP8, 1, DBlack);
X	if (background == 0)
X		exits("out of image memory");

X	einit(Emouse|Ekeyboard);
X	redraw();
X	for (; ; ) {
X		if (ecanmouse()) {
X			mev = emouse();
X			switch (mev.buttons) {
X			case 4:
X				sweeprect = egetrect(3, &mev);
X				if (sweeprect.min.x == 0 && sweeprect.min.y == 0
X				    && sweeprect.max.x == 0 &&
X				    sweeprect.max.y == 0)
X					break;
X				minpnt.re += stepx * (sweeprect.min.x - pmin.x);
X				minpnt.im -= stepy * (sweeprect.min.y - pmin.y);
X				maxpnt.re = minpnt.re + stepx *
X					(sweeprect.max.x - pmin.x);
X				maxpnt.im = minpnt.im - stepy *
X					(sweeprect.max.y - pmin.y);
X				redraw();
X				break;
X			}
X		}
X		if (ecankbd()) {
X			key = ekbd();
X			switch (key) {
X			case DELETEKEY:
X				exits(nil);
X				break;
X			case LESSKEY:
X				nsteps -= 5;
X				redraw();
X				break;
X			case PLUSKEY:
X				nsteps += 5;
X				redraw();
X				break;
X			case RDRAWKEY:
X				minpnt = (Complex){-1.5, 1};
X				maxpnt = (Complex){0.5, -1};
X				redraw();
X				break;
X			}
X		}
X		sleep(100);
X	}
X}
!
echo mandel.h
sed 's/^X//' >mandel.h <<'!'
X// typedef double Fp;
typedef float Fp;
typedef struct Complex {
X	Fp re;
X	Fp im;
X} Complex;

enum {
X	MASRESX = 1024,
X	MASRESY = 768,
X	NCOLS = 0xFF,
X	NSTEPSMIN = 10,

X	MAX = 500,		// module to consider non-convergence
X	MIN = 0.001,		// module to consider convergence

X	MINCOLOR = DBlue,
X	OFFBLUE = 1,
X	MAXCOLOR = DGreen,
X	OFFGREEN = 2,
X	OFFRED = 1,

X	DELETEKEY = '\177',
X	PLUSKEY =  '+',
X	LESSKEY =  '-',
X	RDRAWKEY = 'r',
X};
!
echo mkfile
sed 's/^X//' >mkfile <<'!'
X</$objtype/mkfile

TARG=mandel
OFILES=mandel.$O
HFILES=mandel.h
BIN=$HOME/bin/$objtype
X# LIB=/386/lib/lib387.a

X</sys/src/cmd/mkone
!
echo readme
sed 's/^X//' >readme <<'!'
Mandel is a program to explore the mandelbrot set.  The mandelbrot set
is the set of c numbers in the complex plane on which y=z²+c starting
with z=0 and iterating by making y=z converges.

Quoting Roger Penrose freely, "why can't we explore mathematical
objects as we explore mountains or rivers?".

Mandel can be zoomed in by using the right button mouse in the same
way you resize a window on rio.  You can also make the number of steps
greater or less the with + and - keys.  To redraw the original picture
use the key r.

Suggestions will be welcomed at paurea@gsyc.escet.urjc.es.
X						Saludos,
X								paurea

playing with RESX and RESY on mandel.h will use less memory if you are
on lower resolutions or plan to run mandelbrot in a little window.
!