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

/*
    buddhabrot.c
    Blake Shaw - 11/08
    Adapted from code from Paul Bourke
    
    To compile:
    gcc -Wall buddhabrot.c -o buddhabrot
    
    Example usage:
    ./buddhabrot 500 1000 100 buddha
*/

#define TRUE  1
#define FALSE 0
#define MAX(x,y) (x > y ? x : y)
#define MIN(x,y) (x < y ? x : y)

typedef struct {
   double x,y;
} XY;

// Prototypes
void WriteImage(char *,unsigned int *,int,int);
void WriteData(char *,unsigned int *,int,int);
int Iterate(int,double,double,int *,XY *);
void sigproc();
int doASAP;

int main(int argc,char **argv) 
{
   int i,t,tt,n,ix,iy,NX,NY,NMAX,TMAX;
   char imageFile[100];
   char dataFile[100];
   int particleCount = 0;
   double x,y;
   XY *xyseq = NULL;
    
    signal(SIGINT, sigproc);
    
    // Parse command line arguments
    if (argc == 5) {
        NX = atoi(argv[1]); //Resolution X
        NY = atoi(argv[1]); //Resolution Y
        NMAX = atoi(argv[2]); //Bailout
        TMAX = atoi(argv[3]); //Num iterations in millions
        printf("Size: %dx%d, Bailout: %d, File: %s\n", NX, NY, NMAX, argv[4]);
    } else {
        printf("./buddha res bailout numIts fileName\n");
        return(0);
    }
    
    sprintf(imageFile, "%s-%d-%d-%d.tga", argv[4], NX, NMAX, TMAX);
    sprintf(dataFile, "%s-%d-%d-%d.dat", argv[4], NX, NMAX, TMAX);
    
    // The density plot
    unsigned int *image = NULL;

    // Malloc space for the image and clear to black
    if ((image = (unsigned int *)malloc(NX*NY*sizeof(unsigned int))) == NULL) {
      fprintf(stderr,"Failed to malloc memory for the image\n");
      exit(-1);
    }
    
    // Initialize   
    for (i=0;i<NX*NY;i++){
      image[i] = 0;
    }

    // Malloc space for the sequence
    xyseq = (XY *)malloc(NMAX*sizeof(XY));
    
    srand48( (unsigned int)time( NULL ) );
    
    doASAP = 0;
    
    // Iterate
    for (t=0;t<TMAX;t++) {
        if(doASAP == 1){
            WriteImage(imageFile,image,NX,NY); 
            WriteData(dataFile,image,NX,NY);
            doASAP = 0;
        } else if(doASAP == 2){
            break;
        }
                
        for (tt=0;tt<1000000;tt++) {
              // Choose a random point in same range 
              x = 6 * drand48() - 3;
              y = 6 * drand48() - 3;
    
              // Determine state of this point, draw if it escapes 
              if (Iterate(NMAX,x,y,&n,xyseq)) {
                 for (i=0;i<n;i++) {
                    ix = 0.3 * NX * (xyseq[i].x + 0.5) + NX/2;
                    iy = 0.3 * NY * xyseq[i].y + NY/2;
                    if (ix >= 0 && iy >= 0 && ix < NX && iy < NY) {
                        image[iy*NX+ix]++;
                    }
                }
                particleCount++;
             }
        }
        printf("%d million iterations -- %d particles\n",t+1, particleCount); 
    }

    // Save the result    
    WriteImage(imageFile,image,NX,NY); 
    WriteData(dataFile,image,NX,NY);
    printf("Size: %dx%d, Bailout: %d, File: %s\n", NX, NY, NMAX, argv[4]);
    exit(0);
}

/*
   Write the buddha image to a data file. Triplets for sparse matrix in matlab
*/
void WriteData(char *fname,unsigned int *image,int width,int height)
{
   int i;
   FILE *fptr;

   // Write the file
   fprintf(stderr,"Writing \"%s\"\n",fname);
   if ((fptr = fopen(fname,"w+")) == NULL) {
      fprintf(stderr,"Failed to open output file\n");
      exit(0);
   }

   for (i=0;i<width*height;i++) {
      fprintf(fptr, "%u\n",image[i]);
   }
   fclose(fptr);
}



/*
   Write the buddha image to a minimal TGA file.
   Can be opened with gimp, PhotoShop, etc.
*/
void WriteImage(char *fname,unsigned int *image,int width,int height)
{
   int i;
   float ramp,biggest=0,smallest;
   FILE *fptr;
    
   printf("Writing out image file...\n");
    
   // Find the largest density value
   for (i=0;i<width*height;i++)
      biggest = MAX(biggest,image[i]);
   smallest = biggest;
   for (i=0;i<width*height;i++)
      smallest = MIN(smallest,image[i]);
   fprintf(stderr,"Density value range: %g to %g\n",smallest,biggest);

   // Write the image
   fprintf(stderr,"Writing \"%s\"\n",fname);
   if ((fptr = fopen(fname,"wb")) == NULL) {
      fprintf(stderr,"Failed to open output file\n");
      exit(0);
   }
   
   // TGA header, endian independent
   putc(0,fptr);  /* Length of ID */
   putc(0,fptr);  /* No colour map */
   putc(2,fptr); /* uncompressed RGB  */
   putc(0,fptr); /* Index of colour map entry */
   putc(0,fptr);
   putc(0,fptr); /* Colour map length */
   putc(0,fptr);
   putc(0,fptr); /* Colour map size */
   putc(0,fptr); /* X origin */
   putc(0,fptr);
   putc(0,fptr); /* Y origin */
   putc(0,fptr);
   putc((width & 0x00ff),fptr); /* X width */
   putc((width & 0xff00) / 256,fptr);
   putc((height & 0x00ff),fptr); /* Y width */
   putc((height & 0xff00) / 256,fptr);
   putc(24,fptr);                      /* 24 bit bitmap     */
   putc(0x00,fptr);

   // Raw uncompressed bytes
   for (i=0;i<width*height;i++) {
      ramp = 2*(image[i] - smallest) / (biggest - smallest);
      if (ramp > 1)
         ramp = 1;
      ramp = pow(ramp,0.5);
      fputc((int)(ramp*255),fptr);
      fputc((int)(ramp*255),fptr);
      fputc((int)(ramp*255),fptr);
   }
   fclose(fptr);
}



/*
   Iterate the Mandelbrot and return TRUE if the point escapes
*/
int Iterate(int NMAX,double x0,double y0,int *n,XY *seq)
{
   int i;
   double x=0,y=0,xnew,ynew;

   *n = 0;
   for (i=0;i<NMAX;i++) {
      xnew = x * x - y * y + x0;
      ynew = 2 * x * y + y0;
      seq[i].x = xnew;
      seq[i].y = ynew;
      if (xnew*xnew + ynew*ynew > 10) {
         *n = i;
         return(TRUE);   
      }
      x = xnew;
      y = ynew;
   }

   return(FALSE);
}

void sigproc()
{ 
    char inputc;
    printf("\nEnter command\n");
    printf("(w)rite image, (q)uit and write: ");
    if (scanf("%c", &inputc) == 1){
        switch (inputc) {
            case 'w':
                doASAP = 1;
                break;
            case 'q':
                doASAP = 2;
                break;
            default:
                doASAP = 0;
        }
    } 
}