/*
Copyright (C) 1993 George Leonidas Coulouris

Program	name	: MandelBrute
Description	: Computes an image of (z**2 + c) in the complex plane
Last modified	: 18 Sep 93

To do:
	Implement bounds checking

Revision history:
	1.0 (18 Sep 93) - First fully working version
*/

/* ----------------------------------------------------- */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/*
Notes:
	COLSPERCOMP must not exceed 256.0
	NUMCOLS must equal COLSPERCOMP ** 3
*/

#define NUMCOLS		16777216	/* total number of colors */
#define COLSPERCOMP	256		/* colors per component */

unsigned short iterations;		/* maximum iterations */

double IterateC(double cx,double cy);	/* prototype for IterateC */

/* ----------------------------------------------------- */

void main(int argc,char *argv[])
{
register unsigned short x,y;		/* pixel and raster indices */

unsigned char tmp;			/* temporary counter */

unsigned short	width,			/* pixelmap width */
		height;			/* pixelmap height */

unsigned long tmpr;			/* temporary radius */

double	corners[4],			/* radii of corners */
	maxr,				/* maximum radius */
	scalefactor,			/* color scale factor */
	xmin,				/* real minimum */
	xmax,				/* real maximum */
	xfrac,				/* real increment */
	ymin,				/* imaginary minimum */
	ymax,				/* imaginary maximum */
	yfrac;				/* imaginary increment */

FILE *fp;

					/* banner */
printf("MandelBrute 1.0 (18 Sep 93)\n");
printf("Copyright (C) 1993 George Leonidas Coulouris\n\n");

					/* check arguments */
if (argc != 9)
	{
	printf("usage: MandelBrute <file> <xmin> <xmax> <ymin> <ymax> <iterations> <width> <height>\n\n");
	exit(0);
	}

					/* check & open file */
if ((fp = fopen(argv[1],"w+")) == NULL)
	{
	printf("File could not be opened.\n\n");
	exit(1);
	}

					/* parse arguments */
xmin		= atof(argv[2]);
xmax		= atof(argv[3]);
ymin		= atof(argv[4]);
ymax		= atof(argv[5]);
iterations	= atoi(argv[6]);
width		= atoi(argv[7]);
height		= atoi(argv[8]);

					/* compute corner radii */
corners[0] = IterateC(xmin,ymin);
corners[1] = IterateC(xmin,ymax);
corners[2] = IterateC(xmax,ymin);
corners[3] = IterateC(xmax,ymax);

					/* find maximum radius */
for(tmp=0;tmp<3;tmp++) if (corners[tmp] > corners[tmp+1]) maxr = corners[tmp];

					/* compute increments */
xfrac 		= (xmax-xmin) / width;
yfrac 		= (ymax-ymin) / height;

					/* compute color scale factor */
scalefactor	= (NUMCOLS-1) / maxr;

fprintf(fp,"P6\n");			/* ppm header */
fprintf(fp,"%d\n",width);		/* pixelmap width */
fprintf(fp,"%d\n",height);		/* pixelmap height */
fprintf(fp,"%d\n",((int)COLSPERCOMP-1));	/* maximum component value */

for(y=0;y<height;y++)			/* number crunching loops */
	{
	printf("Line %d, %d%%\r",y,(100*y/height));
	for(x=0;x<width;x++)
		{
		tmpr = (long)(scalefactor * IterateC( (xmin + x*xfrac) , (ymin + y*yfrac) ) );


/*
Notes:
	The following three statements extract the red, green, and blue
components of the pixel, and write them to the ppm file.
*/

		fputc( (int)((COLSPERCOMP-1) * tmpr / (NUMCOLS-1)) ,fp);

		fputc( (int)((tmpr/COLSPERCOMP) % COLSPERCOMP), fp);

		fputc( (int)(tmpr % COLSPERCOMP) ,fp);
		}
	}

					/* Close up and exit */
if (fclose(fp) == EOF)
	{
	printf("File could not be closed.\n");
	exit(1);
	}
else
	{
	printf("Operation complete.\n");
	exit(0);
	}

}

/* ----------------------------------------------------- */

double IterateC(double cx,double cy)
/*
Arguments	: cx, cy
Value		: 0 if c in Mandelbrot set, else escape radius
*/
{
register unsigned short i;		/* Iteration index */

double oldzx,oldzy,zx,zy,r;

oldzx	= 0.0;				/* Initial z seeds */
oldzy	= 0.0;

for(i=0;i<iterations;i++)
	{
	zx = oldzx * oldzx - oldzy * oldzy + cx;
	zy = 2.0 * oldzx * oldzy + cy;
	r = zx * zx + zy * zy;

	if (r > 4.0) return sqrt(r);
	else
		{
		oldzx = zx;
		oldzy = zy;
		}
	}
return 0;
}