#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <errno.h>
#include <math.h>
#include <time.h>

/* The name of the program, for error messages. */
const char *progname;

/* The size of the EPS we will output in PostScript points. */
double width, height;

/* How distant particles can be from a particle in the crystal before they
 * latch on to it.  This decreases as the crystal grows. */
double range;

/* The initial value of 'range' and the ratio of line width to the current
 * range value. */
const double MAX_RANGE = 3;
const double WIDTH_TO_RANGE = 0.23;

/* The 'image' we build up is a collection of points, which are each (except
 * for the first) connected to another point.  We store these in a linked list
 * of dynamically allocated 'Point' objects.  The 'size' decreases as we
 * add points. */
struct Point {
   double x, y;
   double size;
   struct Point *connected_to;
   struct Point *next;
};
struct Point *points = 0;


/* Write an encapsulated PostScript file to the filehandle 'f'. */
void
write_eps (FILE *f)
{
   const struct Point *p;
   time_t t = time(NULL);

   fprintf(f, "%%!PS-Adobe-3.0 EPSF-3.0\n"
              "%%%%BoundingBox: 0 0 %d %d\n"
              "%%%%Creator: %s\n"
              "%%%%CreationDate: %s"
              "%%%%Pages: 1\n"
              "%%%%EndComments\n\n",
           (int) ceil(width), (int) ceil(height), progname, ctime(&t));

   fputs("%%BeginSetup\n"
         "1 setlinecap\n"
         "%%EndSetup\n\n"
         "%%Page: 1 1\n", f);

   for (p = points; p; p = p->next) {
      /* Ignore the single point which isn't connected to another. */
      if (p->connected_to) {
         /* Draw a single line for this point. */
         fprintf(f, "%f setlinewidth\n"
                    "newpath %f %f moveto %f %f lineto stroke\n",
                 p->size, p->x, p->y, p->connected_to->x, p->connected_to->y);
      }
   }

   fputs("\nshowpage\n"
         "%%EOF\n", f);
}

/* Return a random integer number in the range 0 to n-1. */
int
rand_int_over_range (int n)
{
   double dummy;
   return (int) (modf((double) rand() / RAND_MAX, &dummy) * n);
}

/* Return a random floating-point number in the range 0 to n. */
double
rand_double_over_range (double n)
{
   return (double) rand() / RAND_MAX * n;
}

/* Set the values pointed to by 'x' and 'y' to a random location along the
 * side of the image. */
void
random_start_position (double *x, double *y)
{
   /* Decide which side to put the new particle on.  The value of 'side'
    * is 0 for the top of the image, 1 for the right, 2 for the bottom
    * or 3 for the left. */
   int side = rand_int_over_range(4);

   *x = (side == 0 || side == 2) ? rand_double_over_range(width)
                                 : (side == 1) ? width - 1 : 0;
   *y = (side == 1 || side == 3) ? rand_double_over_range(height)
                                 : (side == 2) ? height - 1 : 0;
}

/* If (x, y) is close enough to any points in the crystal to latch on to it,
 * then return the closest point. */
struct Point *
next_to_crystal (double x, double y)
{
   struct Point *p;
   struct Point *closest = 0;
   double dist, closest_dist;
   double x_dist, y_dist;

   for (p = points; p; p = p->next) {
      x_dist = x - p->x;
      y_dist = y - p->y;
      dist = sqrt(x_dist * x_dist + y_dist * y_dist);

      if (!closest || dist < closest_dist) {
         closest = p;
         closest_dist = dist;
      }
   }

   return closest_dist < range ? closest : 0;
}

/* Shift the point (x, y) randomly. */
void
shift_point (double *x, double *y)
{
   *x += rand_double_over_range(range) - 0.5 * range;
   *y += rand_double_over_range(range) - 0.5 * range;

   if (*x < 0)
      *x = width;
   else if (*x > width)
      *x = 0;
   else if (*y < 0)
      *y = height;
   else if (*y > height)
      *y = 0;
}

/* Add a point to the start of the linked-list 'points', given its position
 * and the point it should be connected to. */
void
add_point (double x, double y, struct Point *connected)
{
   struct Point *new_p = malloc(sizeof(struct Point));
   if (!new_p) {
      fprintf(stderr, "%s: error allocating memory for point: %s\n",
              progname, strerror(errno));
      exit(2);
   }

   new_p->x = x;
   new_p->y = y;
   new_p->connected_to = connected;
   new_p->size = range * WIDTH_TO_RANGE;

   /* Attach the new point at the start of the list. */
   new_p->next = points;
   points = new_p;
}

/* Simulate a single particle starting at the edge of the image and moving
 * about randomly until it bumps into part of the crystal. */
void
do_particle ()
{
   double x, y;
   struct Point *p;

   random_start_position(&x, &y);

   while (!(p = next_to_crystal(x, y)))
      shift_point(&x, &y);

   add_point(x, y, p);
}


int
main (int argc, char **argv)
{
   int iterations;
   int i;
   struct Point *p, *p_temp;

   /* Make sure we get different random numbers each time we run this. */
   srand(time(0));

   /* Check and decode the command line arguments. */
   progname = argv[0];
   if (argc != 4) {
      fprintf(stderr, "Usage: %s WIDTH HEIGHT ITERATIONS > out.eps\n",
              progname);
      return 1;
   }

   width = atof(argv[1]);
   height = atof(argv[2]);
   iterations = atoi(argv[3]);

   /* Start with a single pixel set in the middle of the image. */
   add_point(width / 2, height / 2, 0);

   /* Simulate the particles as many times as requested. */
   for (i = 0; i < iterations; ++i) {
      range = MAX_RANGE * (1 - (double) i / iterations);
      range = MAX_RANGE / 10 + range / 10 * 9;

      do_particle();

      if (i && (iterations < 100 || i % (iterations / 100) == 0))
         fprintf(stderr, "Done %d/%d particles (%d%%).\n",
                 i, iterations, (int) ((i * 100.0) / iterations + 0.5));
   }

   write_eps(stdout);

   /* Deallocate the points, just for completeness. */
   for (p = points; p; p = p_temp) {
      p_temp = p->next;
      free(p);
   }

   return 0;
}