/* Swarm trail generator.

   ./swarm <width> <height> > output.pgm
                                         */
#include<math.h>
#include<stdio.h>
#include<stdlib.h>
#include<time.h>
#ifdef _WIN32
   #include<fcntl.h>
   #include<io.h>
#endif

#define PI 3.14159265358979323846264338327950288419716939937510

/* Number of follower points. */
#define FOLLOWERS 512

/* Wait this many iterations for the followers' paths to converge. */
#define STABLIZE_STEPS 2

/* Number of movements steps to run in each simulation step. */
#define SIMULATION_SUBSTEPS 100


typedef struct { double x, y; } XY;

static void MoveToward(double max_velocity,
                       double acceleration,
                       XY target,
                       XY *position,
                       XY *velocity)
{
   double dx = target.x - position->x;
   double dy = target.y - position->y;
   double d = hypot(dx, dy);

   if( d < 1e-6 )
      d = 1e-6;
   dx *= acceleration / d;
   dy *= acceleration / d;

   velocity->x += dx;
   velocity->y += dy;
   d = hypot(velocity->x, velocity->y);
   if( d > max_velocity )
   {
      velocity->x *= max_velocity / d;
      velocity->y *= max_velocity / d;
   }
   position->x += velocity->x;
   position->y += velocity->y;
}

int main(int argc, char **argv)
{
   int width, height, short_edge, i, j, k, x, y;
   unsigned char *pixels;
   XY target, leader_position, leader_velocity;
   XY follower_position[FOLLOWERS], follower_velocity[FOLLOWERS];
   double max_velocity, acceleration, a;

   if( argc != 3 ||
       (width = atoi(argv[1])) <= 0 ||
       (height = atoi(argv[2])) <= 0 )
   {
      return printf("%s <width> <height>\n", *argv);
   }
   if( width >= 0x8000 || height >= 0x8000 )
      return !puts("Output size too large.");

   if( (pixels = (unsigned char*)calloc(width * height, 1)) == NULL )
      return !puts("Not enough memory");

   #ifdef _WIN32
      setmode(STDOUT_FILENO, O_BINARY);
   #endif

   srand(time(NULL));

   /* Initialize random leader position and velocity. */
   leader_position.x = (double)rand() / RAND_MAX * width;
   leader_position.y = (double)rand() / RAND_MAX * height;
   short_edge = width < height ? width : height;
   max_velocity = short_edge * 0.004;
   acceleration = max_velocity * 0.1;
   a = (double)rand() / RAND_MAX * PI * 2.0;
   leader_velocity.x = max_velocity * cos(a);
   leader_velocity.y = max_velocity * sin(a);

   /* Initialize random follower positions. */
   for(i = 0; i < FOLLOWERS; i++)
   {
      follower_position[i].x = (double)rand() / RAND_MAX * width;
      follower_position[i].y = (double)rand() / RAND_MAX * height;
      follower_velocity[i].x = 0.0;
      follower_velocity[i].y = 0.0;
   }

   /* Run swarm simulation. */
   for(i = 0; i < STABLIZE_STEPS + 1 + short_edge / 8; i++)
   {
      /* Generate new target position, preferring one that is
         some distance away from current leader position.     */
      for(j = 0; j < 20; j++)
      {
         target.x = (double)rand() / RAND_MAX * width;
         target.y = (double)rand() / RAND_MAX * height;
         if( hypot(target.x - leader_position.x, target.y - leader_position.y)
               > short_edge * 0.5 )
         {
            break;
         }
      }
      for(j = 0; j < SIMULATION_SUBSTEPS; j++)
      {
         /* Accelerate leader toward target. */
         MoveToward(max_velocity,
                    acceleration,
                    target,
                    &leader_position,
                    &leader_velocity);

         /* Leader point is never drawn, we only observe the movement
            of the leader through followers.                          */

         /* Accelerate follower toward leader. */
         for(k = 0; k < FOLLOWERS; k++)
         {
            MoveToward(max_velocity * 3.7,
                       acceleration * 3.7,
                       leader_position,
                       &follower_position[k],
                       &follower_velocity[k]);

            /* Draw point. */
            if( i > STABLIZE_STEPS )
            {
               x = round(follower_position[k].x);
               y = round(follower_position[k].y);
               if( x >= 0 && x < width && y >= 0 && y < height &&
                   pixels[y * width + x] < 255 )
               {
                  pixels[y * width + x]++;
               }
            }
         }
      }
   }

   /* Normalize and write pixels. */
   j = 1;
   for(i = 0; i < width * height; i++)
      j = j < (int)pixels[i] ? (int)pixels[i] : j;
   printf("P5\n%d %d\n255\n", width, height);
   for(i = 0; i < width * height; i++)
      fputc((int)pixels[i] * 255 / j, stdout);
   free(pixels);
   return 0;
}