[9fans] Little program to explore mandelbrot sets

[paurea@dei.inf.uc3m.es]

[-- Attachment #1: message body text --]
[-- Type: text/plain, Size: 73 bytes --]

It is a bit slow, but suggestions for *simple* solutions are welcome.


[-- Attachment #2: mandel.h --]
[-- Type: application/octet-stream, Size: 622 bytes --]

int addcolor( int colpos, int desvcol, int colbrush);
void redraw(void);
typedef struct Complex {
                     double re;
                     double im;
               } Complex;
Complex compmul(Complex a, Complex b);
Complex compsum(Complex a, Complex b);
Complex iterate(Complex a, int nsteps);
double module (Complex a);
enum{
	MASRESX=1024,
	MASRESY=768,
	DELETEKEY=127,
	NCOLS=0xFF,
	NSTEPSMIN=10,
	MAX=500,  // module to consider non convergence
	MIN=0.001, //module to consider convergence
	MINCOLOR=DBlue,
	OFFBLUE=1,
	MAXCOLOR=DGreen,
	OFFGREEN=2,
	OFFRED=1,
	PLUSKEY=43,
	LESSKEY=45,
	RDRAWKEY=114,
	};

[-- Attachment #3: mandel.c --]
[-- Type: application/octet-stream, Size: 4367 bytes --]

#include <u.h>
#include <libc.h>
#include <draw.h>
#include  <event.h>
#include <cursor.h>
#include <mandel.h>

/* 
mk
mk clean
window mandel
*/

/*
more complete:
void
main(int argc, char* argv[])
*/

Rectangle sweeprect;
Image *palette[NCOLS];
int colbrush;
Point pmin,pmax;
double stepx,stepy;
int nstepsx, nstepsy;
Complex minpnt={-1.5,1},maxpnt={0.5,-1};
Image *background;
int nsteps=NSTEPSMIN;
Image *imagebuf=nil;
Complex pixarray[MASRESX][MASRESY];

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

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

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

	thiscol+=desvcol;

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

	colbrush&=~colormsk;
	colbrush|=thiscol;
	return(colbrush);
}

Complex 
compmul (Complex a, Complex b)
{
	Complex resul;

	resul.re=a.re*b.re-a.im*b.im;
	resul.im=a.re*b.im+a.im*b.re;
	return(resul);
}

Complex 
compsum (Complex a, Complex b)
{
	Complex resul;
	resul.re=a.re+b.re;
	resul.im=a.im+b.im;
	return(resul);
}

double
module (Complex a)
{
	return(sqrt(a.re*a.re+a.im*a.im));
}

Complex
iterate(Complex a, int numsteps)
{	
	Complex z=a;
	Complex res;
	int i;
	for(i=1;i<numsteps;i++){
		z=compsum(compmul(z,z),a);
		if((int)module(z)>MAX){
			res.re=res.im=MAX;
			return(res);
		}	
	if(module(z)<MIN){
			res.re=res.im=MIN;
			return(res);
		}	
		
	}
	res=z;
	return(res);
}

void
	redraw(void)
{
	int isteps,i,j,k,l;
	Point pixel;
	Complex pixelpnt;
	Complex pastpix;

	for(k=0;k<1024;k++)
		for(l=0;l<768;l++){
			pixarray[k][l].re=MIN;
			pixarray[k][l].im=MIN;
		}

	draw(screen,screen->r,background,nil,ZP);
	if(imagebuf!=nil)   
		freeimage(imagebuf);
	imagebuf = allocimage(display, screen->r, CMAP8, 0, DBlack);
			if( imagebuf==0 )
				exits("out of image memory");
	assert(imagebuf);

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

	nstepsx=pmax.x-pmin.x;
	nstepsy=pmax.y-pmin.y;
	stepx=(maxpnt.re-minpnt.re)/nstepsx;
	stepy=(minpnt.im-maxpnt.im)/nstepsy;
	for(isteps=11;isteps<nsteps+11;isteps++){
		for(i=1;i<nstepsx-1;i++){
			for(j=1;j<nstepsy-1;j++){
				pastpix=pixarray[i][j];	
				pixel.x=pmin.x+i;
				pixel.y=pmin.y+j;
				pixelpnt.re=i*stepx+minpnt.re;
				pixelpnt.im=-j*stepy+minpnt.im;
				if(pastpix.re!=MAX && pastpix.im!=MAX){
					if(module(pastpix)!=MIN)
						pixarray[i][j]=iterate(pixelpnt,isteps);
				}
				if(module(pixarray[i][j])<=MAX)
					draw(imagebuf,Rect(pixel.x,pixel.y,pixel.x+1,pixel.y+1),palette[(isteps*19)%NCOLS],nil,ZP);
			}
		}
	draw(screen,screen->r,imagebuf,nil,screen->r.min);
	}
	
}

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

}
 
void
main()
{
	int colmap=MINCOLOR;
	int key;
	Mouse mev;
	int isteps;
	
	if (initdraw(nil, nil, "mandel") < 0)
		exits("initdraw failed");
/*
I allocate  the color palette to draw 
*/
	for(isteps=10;isteps<NCOLS;isteps++){
			palette[isteps] = allocimage(display, Rect(0,0,1,1), CMAP8, 1, colmap);
			if( palette[isteps]==0 )
				exits("out of image memory");

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

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

	einit(Emouse|Ekeyboard);

redraw();
	
	for (;;) {
		if(ecanmouse()){
			mev = emouse();
			switch (mev.buttons){
				case 4:					
					sweeprect=egetrect(3,&mev);
					if((sweeprect.min.x==0) && (sweeprect.min.y==0) 
						&& (sweeprect.max.x==0) && (sweeprect.max.y==0)){
						break;
					}
					minpnt.re=minpnt.re+stepx*(sweeprect.min.x-pmin.x);
					minpnt.im=minpnt.im-stepy*(sweeprect.min.y-pmin.y);
					maxpnt.re=minpnt.re+stepx*(sweeprect.max.x-pmin.x);
					maxpnt.im=minpnt.im-stepy*(sweeprect.max.y-pmin.y);
					redraw();
				break;
			}
		}
		if(ecankbd()){		
			key= ekbd();
			switch(key){
				case DELETEKEY:
					exits(nil);
				break;
				case LESSKEY:
					nsteps-=5;
					redraw();
				break;
				case PLUSKEY:
					nsteps+=5;
					redraw();
				break;
				case RDRAWKEY:
					minpnt.re=-1.5;
					minpnt.im=1;
					maxpnt.re=0.5;
					maxpnt.im=-1;
					redraw();
				break;
			}
		}
	sleep(100);
	}
}



[-- Attachment #4: readme --]
[-- Type: application/octet-stream, Size: 752 bytes --]

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.
						Saludos,
								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.

[-- Attachment #5: mkfile --]
[-- Type: application/octet-stream, Size: 121 bytes --]

</$objtype/mkfile
TARG=mandel
OFILES=mandel.$O
BIN=.
HFILES=mandel.h
mandel: install
	touch mandel

</sys/src/cmd/mkone


[-- Attachment #6: .signature --]
[-- Type: text/plain, Size: 99 bytes --]



-- 
                 Saludos,
                         Gorka

"Curiosity sKilled the cat"