#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include <sys/types.h>
#include <time.h>
#include "paulslib.h"
#include "bitmaplib.h"

#define NX 500
#define NY 500
#define N 1000000
#define SCALE (NX / 6)

int main(int argc,char **argv)
{
   int i,n,ix,iy;
   int type,s;
   double x=0.5,y=0,x1,y1,xa,ya,r,w=0;
   BITMAP *image,white={255,255,255},black={0,0,0};
   BITMAP colour;
   char fname[64];
   FILE *fptr;

   if (argc < 2) {
      fprintf(stderr,"Usage: %s type\n",argv[0]);
      exit(-1);
   }
   type = atoi(argv[1]);
   if (type < 0 || type > 2) {
      fprintf(stderr,"Type shold be 0, 1, or 2\n");
      exit(-1);
   }

   /* Create the image */
   image = CreateBitmap(NX,NY);
   EraseBitmap(image,NX,NY,white);
   srand(time(NULL));

   for (n=0;n<N;n++) {
      if ((n % (N/10)) == 0)
         fprintf(stderr,"%d\n",n);
      switch (type) {
      case 0:
         xa = -4 * x + cos(4 * w);
         ya = -4 * y + sin(4 * w);
         break;
      case 1:
         xa = -4 * x + 1;
         ya = -4 * y;
         break;
      case 2:
         xa = -4 * x + cos(w);
         ya = -4 * y + sin(w);
         break;
      }
      r = sqrt(xa*xa + ya*ya);
      w = atan2(ya,xa);
      switch (type) {
      case 0:
         if (rand() % 2 == 0) {
            x1 = - cos(w) - sqrt(r) * cos(w/2) / 2;
            y1 = - sin(w) + sqrt(r) * sin(w/2) / 2;
         } else {
            x1 = - cos(w) + sqrt(r) * cos(w/2) / 2;
            y1 = - sin(w) - sqrt(r) * sin(w/2) / 2;
         }
         break;
      case 1:
         s = 2 * (rand() % 2) - 1;
         if (rand() % 2 == 0) {
            x1 = - s * cos(w) - sqrt(r) * cos(w/2) / 2;
            y1 =   s * sin(w) - sqrt(r) * sin(w/2) / 2;
         } else {
            x1 = - s * cos(w) + sqrt(r) * cos(w/2) / 2;
            y1 =   s * sin(w) + sqrt(r) * sin(w/2) / 2;
         }
         break;
      case 2:
         if (rand() % 2 == 0) {
            x1 = - cos(2*w) - sqrt(r) * cos(w/2) / 2;
            y1 = - sin(2*w) + sqrt(r) * sin(w/2) / 2;
         } else {
            x1 = - cos(2*w) + sqrt(r) * cos(w/2) / 2;
            y1 = - sin(2*w) - sqrt(r) * sin(w/2) / 2;
         }
      }
      x = x1;
      y = y1;
      if (n < 100)
         continue;
      ix = x * SCALE + NX/2;
      iy = y * SCALE + NY/2;
      DrawPixel(image,NX,NY,ix,iy,black);
   }

   /* Write the image to a ppm file */
   sprintf(fname,"roger1_%d.tif",type);
   if ((fptr = fopen(fname,"w")) == NULL) {
      fprintf(stderr,"Unable to open bitmap file\n");
      exit(0);
   }
   WriteBitmap(fptr,image,NX,NY,5);
   fclose(fptr);
   DestroyBitmap(image);
}