/*
  born.c
  
  By Bill Kendrick
  For Astronomy 305 Assignment 3 - Stars
  
  kendrick@zippy.sonoma.edu
  http://zippy.sonoma.edu/kendrick/projects/astro305-3/
  
  Simple simulation of gas falling into itself to form a sphere.
  All points are stored in two-dimensional, real number (not integer) space.

  Creates multiple GIF files which can easily be turned into an animated GIF.
  
  October 24, 1996 - October 28, 1996
  (Upgraded to PROPERLY simulate physics (oops!))
*/


/* Library #includes: */

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


/* Program option #defines: */

#define MAX 1000                                   /* Max. particles allowed */
#define PPMTOGIF "/usr/local/bin/pbmplus/ppmtogif" /* Path to "ppmtogif"     */


/* MAIN: */

int main(int argc, char * argv[])
{
  int seed, num_particles, frames_wanted;     /* User options              */
  int frame_spacing;
  float part_x[MAX], part_y[MAX];             /* Particle positions        */
  float part_xm[MAX], part_ym[MAX];           /* Particle momentums        */
  int part_m[MAX];                            /* Particle masses           */
  FILE * fi;                                  /* A stream for output files */
  char filename[100];                         /* Output filename buffer    */
  unsigned char bitmap[100][100][3];          /* 100 x 100 x 24bit buffer  */
  int x, y, rgb;                              /* Bitmap-manipulation vars  */
  int p, p2, p3, f, s, ok;                    /* Algorithm variables       */
  float fx, fy, new_x, new_y;
  
  
  /* Get the user options
     (the random number generation seed to use,
     the number of particles to simulate,
     the number of images to generate, and
     the number of frames between each image): */
  
  if (argc != 5)
    {
      fprintf(stderr, "Usage: %s [random seed] [particles] [frames wanted] ",
	      argv[0]);
      fprintf(stderr, "[frame spacing]\n");
      exit(1);
    }
  
  seed = atoi(argv[1]);
  num_particles = atoi(argv[2]);
  frames_wanted = atoi(argv[3]);
  frame_spacing = atoi(argv[4]);
  
  
  if (num_particles > MAX)
    {
      fprintf(stderr, "Most particles allowed is %d.\n", MAX);
      exit(1);
    }
  
  
  /* "Initialize" the particles: */
  
  srand(seed);
  
  for (p = 0; p < num_particles; p++)
    {
      do
	{
	  /* X and Y positions between 0 and 99.00: */
	  
	  part_x[p] = (float) (rand() % 10000) / 100;
	  part_y[p] = (float) (rand() % 10000) / 100;
	  
	  
	  /* Start with random momentum: */
	  
	  part_xm[p] = (((float) (rand() % 200)) / 100) - 1;
	  part_ym[p] = (((float) (rand() % 200)) / 100) - 1;
	  
	  
	  /* Masses between 1 and 5 "units": */
	  
	  part_m[p] = (rand() & 5) + 1;
	  
	  
	  /* Make sure particle isn't in the same place as another
	     particle! */
	  
	  ok = 1;
	  for (p2 = 0; p2 < p; p2++)
	    {
	      if (floor(part_x[p]) == floor(part_x[p2]) &&
		  floor(part_y[p]) == floor(part_y[p2]))
		ok = 0;
	    }
	}
      while (ok == 0);
    }
  
  
  /* Create the images:  (MAIN SIMULATION LOOP) */
  
  for (f = 0; f < frames_wanted; f++)
    {
      /* Blank the current bitmap to black: */
      
      for (y = 0; y < 100; y++)
	for (x = 0; x < 100; x++)
	  for (rgb = 0; rgb < 3; rgb++)
	    bitmap[y][x][rgb] = 0;
      
      
      /* Draw the particles onto the bitmap: */
      
      for (p = 0; p < num_particles; p++)
	for (rgb = 0; rgb < 3; rgb++)
	  {
	    y = (int) part_y[p];
	    x = (int) part_x[p];
	    
	    /* Luminosity depends on mass: */
	    
	    if (y >= 0 && y < 100 && x >= 0 && x < 100)
	      bitmap[y][x][rgb] = 25 * part_m[p] + 120;
	  }
      
      
      /* Dump the current frame into a file: */
      
      sprintf(filename, "%s > frame%.3d.gif", PPMTOGIF, f);
      
      fi = popen(filename, "w");
      
      if (fi == NULL)
	{
	  /* Totally abort on a bad file stream!!! */
	  
	  perror(filename);
	  exit(1);
	}
      
      
      /* Send PPM header: */
      
      fprintf(fi, "P6\n100\n100\n255\n");
      
      
      /* Send our bitmap's data (raw): */
      
      for (y = 0; y < 100; y++)
	for (x = 0; x < 100; x++)
	  for (rgb = 0; rgb < 3; rgb++)
	    fputc(bitmap[y][x][rgb], fi);
      
      pclose(fi);
      
      
      /* Move "frame_spacing" more frames... */
      
      for (s = 0; s < frame_spacing; s++)
	{
	  /* Move the particles: */
	  
	  for (p = 0; p < num_particles; p++)
	    {
	      for (p2 = 0; p2 < num_particles; p2++)
		{
		  if (p != p2)          /* Don't let us influence ourselves! */
		    {
		      /* VERY POOR MATH to determine momentum of the current
			 particle, depending on the other particles. */
		      
		      fx = (float) ((part_x[p2] - part_x[p]) / (5000 /
								part_m[p2]));
		      fy = (float) ((part_y[p2] - part_y[p]) / (5000 /
								part_m[p2]));
		      
		      part_xm[p] = part_xm[p] + (fx / 100);
		      part_ym[p] = part_ym[p] + (fy / 100);
		    }
		}
	      
	      
	      /* Move particle depending on its momentum: */
	      
	      part_x[p] = part_x[p] + part_xm[p];
	      part_y[p] = part_y[p] + part_ym[p];
	      
	      
	      /* "Bounce" on collisions */
	      
	      ok = 1;
	      for (p2 = 0; p2 < p; p2++)
		{
		  if (floor(part_x[p]) == floor(part_x[p2]) &&
		      floor(part_y[p]) == floor(part_y[p2]))
		    {
		      part_x[p] = part_x[p] - part_xm[p];
		      part_y[p] = part_y[p] - part_ym[p];
		      part_xm[p] = 0;
		      part_ym[p] = 0;
		    }
		}
	    }
	}
    }
  
  return(0);
}