/*
Copyright (C) 1993 George Leonidas Coulouris

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

To do:
	Fix color scaling bug
*/

/* - Preprocessor ------------------------------------------------ */

/*
	USE881 is only set for Amiga users with SAS/C and MC68881s
	COLSPERCOMP <= 256
	NUMCOLS = COLSPERCOMP ** 3
*/

#define USE881

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

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

#ifdef USE881
#include <m68881.h>
#else
#include <math.h>
#endif

/* - Declarations ------------------------------------------------ */

unsigned short iterations;	/* maximum iterations */

double IterateC(double cx,double cy);

/* - Main -------------------------------------------------------- */

void main(int argc,char *argv[])
{

register unsigned short
	x,				/* pixel index */
	y,				/* raster index */
	width,				/* pixelmap width */
	height;				/* pixelmap height */

register unsigned long
	tmpr;				/* temporary radius */

unsigned char
	tmp;				/* temporary counter */

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.2 (24 Sep 93)\n");
printf("Copyright (C) 1993 George Leonidas Coulouris\n\n");

					/* check argument list */
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]);

					/* check arguments */
if ((xmin >= xmax) || (ymin >= ymax) )
	{
	printf("Bad coordinates.\n\n");
	exit(1);
	}

if (iterations <= 0)
	{
	printf("Bad iterations.\n\n");
	exit(1);
	}

if ( (width <= 0) || (height <= 0) )
	{
	printf("Bad image dimensions.\n\n");
	exit(1);
	}

					/* 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++)
	maxr = (corners[tmp] > corners[tmp+1]) ? corners[tmp] : corners[tmp+1];

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

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

/*
	The following four statements write the PPM header.
*/

fprintf(fp,"P6\n");
fprintf(fp,"%d\n",width);
fprintf(fp,"%d\n",height);
fprintf(fp,"%d\n",(COLSPERCOMP-1));

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

/*
	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);
	}

}

/* - Functions --------------------------------------------------- */

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

double
	oldzx = 0,			/* Initial z seeds */
	oldzy = 0,
	zx,
	zy,
	r;

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