#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "paulslib.h"
#include "bitmaplib.h"
#include "complexlib.h"

/*
   Crude and inefficient, but functional, little program to calculate the
   Mandelbrot-like sets including deep zooms.
*/

void WriteImage(int,int,int,int,char *,int);
int Iterate(double,double,long);

#define NX 650
#define NY 650
#define NA 4

int main(int argc,char **argv) 
{
   int i,j,ax,ay,n,nmax=0;
   double x,y,x0 = -0.5,y0 = 0,zoom = 1;
   FILE *fptr;
   char fname[128];
   int *grid,gridmin,gridmax;
   long itermax = 100;

   /* Command line stuff */
   if (argc < 4) {
      fprintf(stderr,"Usage: %s x0 y0 zoom\n",argv[0]);
      exit(0);
   }
   x0 = atof(argv[1]);
   y0 = atof(argv[2]);
   zoom = atof(argv[3]);
   itermax *= zoom;
   fprintf(stderr,"----------------------------------------------------\n");

   /* Create the storage */
   grid = malloc((NX+1)*(NY+1)*sizeof(int));
   for (i=0;i<(NX+1)*(NY+1);i++)
      grid[i] = 0;
   gridmin = itermax;
   gridmax = 0;

   /* Form the image */
   for (j=0;j<=NY;j++) {
      for (i=0;i<=NX;i++) {
         nmax = 0;
         for (ax=0;ax<NA;ax++) {
            for (ay=0;ay<NA;ay++) {
               x = 4 * (i + ax/(double)NA - NX/2) / (NX * zoom) + x0;
               y = 4 * (j + ay/(double)NA - NY/2) / (NY * zoom) + y0;
               n = Iterate(x,y,itermax);
               //nmax = MAX(n,nmax);
               nmax += n;
            }
         }

         /* Shell display */
         if (i % (NX/70) == 0 && j % (NY/25) == 0) {
            if (nmax > 4*gridmax/5)
               fprintf(stderr,"*");
            else if (nmax > 3*gridmax/5)
               fprintf(stderr,"=");
            else if (nmax > 2*gridmax/3)
               fprintf(stderr,":");
            else if (nmax > gridmax/5)
               fprintf(stderr,".");
            else
               fprintf(stderr," ");
         }

         /* Update the grid */
         grid[j*(NX+1)+i] = nmax;
         gridmin = MIN(gridmin,nmax);
         gridmax = MAX(gridmax,nmax);
      }

      /* Shell display */
      if (j % (NY/25) == 0)
         fprintf(stderr,"\n");
   }

   /* Stats */
   fprintf(stderr,"x: %g, y:%g, zoom: %g\n",x0,y0,zoom);
   fprintf(stderr,"\nRange: %d -> %d\n",gridmin,gridmax);

   /* Autoscale */
   for (i=0;i<(NX+1)*(NY+1);i++)
      grid[i] -= gridmin;
   gridmax -= gridmin;
   gridmin = 0;
   if (gridmax > 255) {
      for (i=0;i<(NX+1)*(NY+1);i++)
         if (grid[i] > 255)
            grid[i] = 255;
      gridmax = 255;
   }

   /* Save to a raw greyscale file */
   if (gridmax - gridmin > 10) {
      sprintf(fname,"%g_%g_%g.raw",x0,y0,zoom);
      for (i=0;i<strlen(fname);i++)
         if (fname[i] == '-')
            fname[i] = 'm';
      fptr = fopen(fname,"w");
      for (i=0;i<(NX+1)*(NY+1);i++) {
         n = 255 - 255 * (grid[i] - gridmin) / (gridmax - gridmin);
         fputc(n,fptr);
      }
      fclose(fptr); 
   } else {
      fprintf(stderr,"Warning, image not saved\n");
      fprintf(stderr,"Image possibly entirely within the set\n");
      fprintf(stderr,"or entirely outside the set.\n");
      fprintf(stderr,"Increase \"itermax\" and/or antialiasing\n");
   }

   exit(0);
}

/*
   Iterate the equation
*/
int Iterate(double x0,double y0,long imax)
{
   COMPLEX p={0.0,0.0},pa,pb,p2,pnew,c;
   int samecounter = 0;
   int i;
   
   c.real = x0;
   c.imag = y0;
   p = c;
   for (i=0;i<imax;i++) {
      p2 = Cmult(p,p);

      pa = p2; pa.real += 1;
      pa = Cmult(pa,c);
      pa = Cmult(pa,p2);

      pb = p2; pb.real -= 1;
      pb = Cmult(pb,pb);
      
      pnew = Cdiv(pa,pb);

      if (pnew.real*pnew.real + pnew.imag*pnew.imag > 1e10)
         return(imax);
      if (ABS(pnew.real-p.real) < 0.00001 && ABS(pnew.imag-p.imag) < 0.00001) 
         samecounter++;
      else
         samecounter = 0;
      if (samecounter > 4)
         return(i);

      p.real = pnew.real;
      p.imag = pnew.imag;
   }

   return(imax);
}