#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <sys/time.h>

#define NX 2000
#define NY 2000
#define N 1000
#define MILLION 1000000
#define SCALE (NX / 7)

int main(int argc,char **argv)
{
   int i,j,m,n,c,ix,iy,biggest=0,index;
   unsigned short iv;
   double x=1,y=1,x1,y1;
   double a[25],b[25],v;
   char fname[64];
   int *image = NULL;
   FILE *fptr;

   if (argc < 2) {
      fprintf(stderr,"Usage: %s m\n",argv[0]);
      exit(-1);
   }
   m = atoi(argv[1]);

   // Create the image 
   image = malloc(NX*NY*sizeof(int));
   for (i=0;i<NX*NY;i++)
      image[i] = 0;
   srand(time(NULL));

   // Initialisation
   for (i=0;i<m;i++) {
      a[2*i] = cos(2*M_PI * (i+1) / (double)m);
      b[2*i] = sin(2*M_PI * (i+1) / (double)m);
   }
   for (i=0;i<m;i++) {
      a[2*i+1] = 0.5*(cos(2*M_PI * (i+1) / (double)m) + cos(2*M_PI * i / (double)m));
      b[2*i+1] = 0.5*(sin(2*M_PI * (i+1) / (double)m) + sin(2*M_PI * i / (double)m));
   }
   
   // Iterate
   for (n=0;n<N;n++) {
      if ((n % (N/50)) == 0)
         fprintf(stderr,"%d\n",n);
      for (j=0;j<MILLION;j++) {
         c = rand() % (2*m);
         x1 = (x/(x*x + y*y))/3 + a[c];
         y1 = (-y/(x*x + y*y))/3 + b[c];
         x = x1 / (x1*x1 + y1*y1);
         y = y1 / (x1*x1 + y1*y1);
         ix = x * SCALE + NX/2;
         iy = y * SCALE + NY/2;
         if (ix >= 0 && iy >= 0 && ix < NX && iy < NY) {
            index = iy*NX+ix;
            image[index]++;
            if (image[index] > biggest)
               biggest = image[index];
         }
      }
   }
   fprintf(stderr,"Biggest: %d\n",biggest);

   // Write the image to a raw file 
   sprintf(fname,"star_%d.raw",m);
   if ((fptr = fopen(fname,"w")) == NULL) {
      fprintf(stderr,"Unable to open output file\n");
      exit(0);
   }
   for (i=0;i<NX*NY;i++) {
      v = image[i] / (float)biggest; // 0 ... 1
      iv = 65535*pow(v,0.5);
      fwrite(&iv,sizeof(short),1,fptr);
   }
   fclose(fptr);

   exit(0);
}