#include <math.h>
#include <stdio.h>
#include "mandel.h"

/* Mandelbrot Set Generator                                           */
/* From page 8, Micro Cornucopia, #39, Jan-Feb 1988                   */
/* mandel(-2.0,0.5,-1.25,1.25) for the picture on page 7.             */
/* modified as per page 49, Micro Cornucopia #43, Setp-Oct 1988.      */
/*                                                                    */
/* modified nx,ny control and delta values and incrementation to be   */
/* consistent with surface plotting routine.  scaling also added.     */

#define sqr(x) (x*x)
#define BIGNUM 1.7e38

const int max_size = 4;                                    /* usually */

void mandel(float Pmin, float Pmax, float Qmin, float Qmax,int max_its,
	    int nx, int ny, float **z,float *zmin,float *zmax)

{
  int color, row, col;
  float P, Q, modulus, deltaP, deltaQ,
	Xcur, Xlast, Ycur, Ylast;
  float scalefact,temp;

/* scale the vertical values so it will look nice as a surface:       */

  scalefact = (Pmax - Pmin);
  if ((Qmax - Qmin) > scalefact) scalefact = Qmax - Qmin;
  scalefact = scalefact / max_its;

  *zmin =  BIGNUM;
  *zmax = -BIGNUM;

  deltaP = (Pmax - Pmin) / (nx-1);   /* determine real axis increment */
  deltaQ = (Qmax - Qmin) / (ny-1);   /* determine imaginary increment */

  P = Pmin;
  for (col = 0; col < nx; col++)  /* look at each point in the region */
    {
      Q = Qmin;
      printf("col = %d\n",col);
      for (row = 0; row < ny; row++)
        {
          Xlast = Ylast = modulus = 0.0;
          color = 0;
          while ((modulus < max_size) && (color < max_its))
              /* loop until function blows up or until max_its */
            {
              /* find components of next Z */
              Xcur = sqr(Xlast) - sqr(Ylast) + P;
              Ycur = 2 * Xlast * Ylast + Q;
              color++;
              Xlast = Xcur;                          /* update last Z */
              Ylast = Ycur;
              modulus = sqr(Xcur) + sqr(Ycur);      /* find size of Z */
	    }
	    temp = (float)color * scalefact;
	    z[col][row] = temp;
	    if (temp < *zmin) *zmin = temp;
	    if (temp > *zmax) *zmax = temp;

	    Q += deltaQ;
        }
        P += deltaP;
    }
}  /* void mandel() */


