/* wobblemovie.c        Oliver Knill, Austin, TX, February 8, 2000           */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>                               // initialisation of rand    

#define N_BODIES  120                           // number of bodies to evolve
#define N_MAX     1000                          // number of frames in movie

#define FILE_INDEX(x) ((x) < 10) ? "00" : (((x) < 100) ? "0" : ""), (x)
#define ERR "/dev/null"

typedef float particlearray[N_BODIES][3];       // position and 
typedef float momentumarray[N_BODIES][3];       // momentum of bodies
typedef float distancearray[N_BODIES][N_BODIES];// distance array of particles

particlearray position;
particlearray momentum;
distancearray distance;

/*---------------------------------------------------------------------------*/
/* function random number generator                                          */
/*---------------------------------------------------------------------------*/
void initialize_number_generator(void) {
  srand((unsigned) time(NULL));
}

float R(void) {
  return(1.0*random()/RAND_MAX);   // Random number in [0,1]
}

float pow2(float a) { return(a*a); }

/*---------------------------------------------------------------------------*/
/* function update_distance                                                  */
/*---------------------------------------------------------------------------*/
void update_distances(void){
int i,j;
  for (i=0; i<N_BODIES; i++) for (j=0; j<N_BODIES; j++) {
     distance[i][j] = sqrt(pow2(position[i][0]-position[j][0])
                         + pow2(position[i][1]-position[j][1])
                         + pow2(position[i][2]-position[j][2]) );
  }
}

/*---------------------------------------------------------------------------*/
/* function initialize_array                                                 */
/*---------------------------------------------------------------------------*/
void initialize_array(void){
int i,j;
  for (i=0; i<N_BODIES; i++) for (j=0; j<3; j++) {
    position[i][j]=2.0*(R()-0.5);
    momentum[i][j]=1.0*(R()-0.5);
  }
  update_distances();
}

/*---------------------------------------------------------------------------*/
/* function evolve_particles                                                 */
/*---------------------------------------------------------------------------*/
void evolve_particles(void){
int k,l;
int t; 
float dt=0.002;                      // time step 
int T=10;                            // number of Euler steps until update

for (t=0; t<T; t++) {
  for (k=0; k<N_BODIES; k++) for (l=0; l<3; l++) {
        position[k][l] = position[k][l]+dt*momentum[k][l];
        momentum[k][l] = momentum[k][l]-dt*position[k][l];
  }
} 
update_distances();
}

/*---------------------------------------------------------------------------*/
/* function `make_filename'                                                  */
/*---------------------------------------------------------------------------*/
char *make_filename(int n, char *ext){
  char *basename = "wobble";
  char *filename;

  filename = (char *) malloc((strlen(basename)+strlen(ext) + 5)*sizeof(char));
  sprintf(filename,"%s%s%i.%s",basename,FILE_INDEX(n),ext);
  return(filename);
}

/*---------------------------------------------------------------------------*/
/* Function `make_wobble'                                                    */
/*---------------------------------------------------------------------------*/
void make_wobble(int i, char *run){
  char *giffile, *povfile;
  FILE *pov;
  float f=0.2;                              // treshold for drop
  float rr=0.3;                             // radius for drop
  
  povfile = make_filename(i,"pov");
  giffile = make_filename(i,"gif");
  
  printf("Writing Povray file %s\n",povfile);

  pov = fopen(povfile,"w");
  fprintf(pov,"camera { location <0,0,-5> up y right (4/3)*x look_at <0,0,5> }\n"
	  "light_source  { <-100,50,-200> color rgb <1,1,1> } \n" 
	  "light_source  { <100,-50,-200> color rgb <1,1,0> } \n" 
          "#macro d(xx,yy,zz,rr) \n"
          "  sphere { <xx,yy,zz>,rr,rr } \n"
          "#end  \n"
	  "blob { \n"
	  "threshold "
	  "%f\n"
          "#declare i=0; \n",f);
  for (i=0; i<N_BODIES; i++) {
         fprintf(pov, "d(%f,%f,%f,%f) \n",
            position[i][0],position[i][1],position[i][2],rr);
  }
  fprintf(pov,"pigment { color rgbt <1.0,1.0,0.0,0.3>}"
          "finish  { phong 1 }\n" 
          "} \n"); 
  fclose(pov);

  printf("Running povray on %s",povfile);
  sprintf(run,"povray good.ini %s;",povfile);
  system(run);
}

/*---------------------------------------------------------------------------*/
/* Function `main'                                                           */
/*---------------------------------------------------------------------------*/

int main() {
  int i = 0, n_max = N_MAX;
  char run[256];

  initialize_number_generator();
  initialize_array();

  for(i = 1; i < n_max; i++) {
    make_wobble(i,run);
    evolve_particles(); 
  }
  printf("All done.\n");

  return 0; 
}