#ifndef lint
static char *sccsid = "%W% (UKC) %G%";
#endif  lint

/*
 * Compute Mandelbrot fractal pattern.
 *
 * we keep this routine in a separate file so that it may be optimised,
 * since this is where the program spends all its time.
 */

#include "mandel.h"

short		pixel[MAX_NPIXEL];	/* Stores one pixel row		*/
double		gap[MAX_NPIXEL];	/* Real axis position		*/
extern	int	hist[];			/* Pixel density histogram	*/
extern	char	line[];			/* General text work area	*/
extern	FILE	*pixfd;			/* Pixel file (binary)		*/

/*
 * Compute the Mandelbrot set.
 */
process()
{
	register int	count;		/* Inner-loop counter		*/
	register int	j;		/* Column counter		*/
	register int	i;		/* Row counter			*/
	double	z_real, z_imag;		/* Pixel accumulator		*/
	double	z2_real;		/* z_real ** 2			*/
	double	z2_imag;		/* z_imag ** 2			*/
	double	c_real, c_imag;		/* Pixel position (constant)	*/
	double	float_pixels;		/* To computer pixel position	*/

	/*
	 * Precompute the position of each pixel on the real axis.
	 * This loop should not be "unrolled" to a succession of
	 * additions as that would lose accuracy.  Given what
	 * follows, the cost isn't excessive.
	 */
	float_pixels = npixel;
	for (j = 0; j < npixel; j++)	/* Precompute column (real) gap	*/
	    gap[j] = start_real + (side * j / float_pixels);
	for (i = 1; i <= npixel; i++) {
	    c_imag = start_imag - (side * i / float_pixels);
	    for (j = 0; j < npixel; j++) {
		/*
	 	 * Compute one point (pixel) of the Mandelbrot set:
		 * 1. Set z to 0 + 0i and
		 *    set c to the pixel location (real, imag)
		 * 2. Perform step 3 until either
		 *    1. count reaches the selected number of iterations or
		 *    2. the "size" of z exceeds 2.0, where size is
		 *       defined as sqrt(z_real**2 + z_imag**2)
		 *       (we don't bother with the square root.)
		 * 3. z = z**2 + c;
		 *    count = count + 1;
		 *    size = size of z as defined above.
		 * 4. Set pixel[i,j] to the number of iterations (count).
		 */
		c_real = gap[j];
		z_real = c_real;
		z_imag = c_imag;
		for (count = 0; count < niter; count++) {
		    if (((z2_real = (z_real * z_real))
		       + (z2_imag = (z_imag * z_imag))) > 4.0)
			break;
		    z_imag = (z_real * z_imag * 2.0) + c_imag;
		    z_real = z2_real - z2_imag + c_real;
		}
		pixel[j] = count;		/* Store the pixel	*/
		hist[count] += 1;		/* Gray-scale histogram	*/
	    }					/* All columns this row	*/
	    fwrite(pixel, npixel * sizeof (pixel[0]), 1, pixfd);
	    if (ferror(pixfd)) {
		perror("pixel file write error");
		return;
	    }
	}					/* All rows in picture	*/
}