/* L-14 MCS 572 Friday 10 Feb 2006: simple program for the n-body problem,
 * compile as "gcc -o /tmp/nbody nbody.c -lm" */

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

void simulate ( int n, int t, double dt );
/* simulates the n-body problem for t time steps of size dt */

void read_bodies ( int n, double m[n], double x[n], double y[n], double z[n] );
/* reads n masses and initial positions (x[i],y[i],z[i]), i=0,1,..,n-1 */

void write_bodies ( int n, double m[n], double x[n], double y[n], double z[n],
                    double vx[n], double vy[n], double vz[n] );
/* writes the n masses with their corresponding coordinates of their
 * positions and velocities */

double square_of_distance ( double x1, double y1, double z1,
                            double x2, double y2, double z2 );
/* returns the square of the distance between (x1,y1,z1) and (x2,y2,z2) */

void gforce ( int n, double m[n], double x[n], double y[n], double z[n],
              double ax[n], double ay[n], double az[n] );
/* computes for every body at (x[i],y[i],z[i]), i=0,1,..,n-1,
 * the components of the acceleration (ax[i],ay[i],az[i]) */

void move ( int n, double x[n], double y[n], double z[n],
            double vx[n], double vy[n], double vz[n], 
            double ax[n], double ay[n], double az[n], double dt );
/* given position (x[i],y[i],z[i]), velocities (vx[i],vy[i],vz[i]),
 * and acceleration (ax[i],ay[i],az[i]), for i=0,1,..,n-1, the position
 * and velocity of every body is updated for one time step of size dt */

int main ( int argc, char *argv[] )
{
   int n,t;
   double dt;

   printf("Give the number of bodies : "); scanf("%d",&n);
   printf("Give the number of time steps : "); scanf("%d",&t);
   printf("Give the size of one time step : "); scanf("%lf",&dt);
   simulate(n,t,dt);

   return 0;
}

void simulate ( int n, int t, double dt )
{
   double m[n],x[n],y[n],z[n],vx[n],vy[n],vz[n],ax[n],ay[n],az[n];
   int i;

   read_bodies(n,m,x,y,z);
   for(i=0; i<n; i++)
   {
      vx[i] = 0.0;
      vy[i] = 0.0;
      vz[i] = 0.0;
   }
   write_bodies(n,m,x,y,z,vx,vy,vz);
   for(i=0; i<t; i++)
   {
      gforce(n,m,x,y,z,ax,ay,az);
      move(n,x,y,z,vx,vy,vz,ax,ay,az,dt);
      write_bodies(n,m,x,y,z,vx,vy,vz);
   }
}

void read_bodies ( int n, double m[n], double x[n], double y[n], double z[n] )
{
   int i;

   for(i=0; i<n; i++)
   {
      printf("Give mass for body %d : ",i); scanf("%lf",&m[i]);
      printf("  and its x-coordinate : ");  scanf("%lf",&x[i]);
      printf("  and its y-coordinate : ");  scanf("%lf",&y[i]);
      printf("  and its z-coordinate : ");  scanf("%lf",&z[i]);
   }
}

void write_bodies ( int n, double m[n], double x[n], double y[n], double z[n],
                    double vx[n], double vy[n], double vz[n] )
{
   int i;

   for(i=0; i<n; i++)
   {
      printf("Body %d with mass %.3e is at (%.3e,%.3e,%.3e),\n",
             i,m[i],x[i],y[i],z[i]);
      printf("        with velocity vector (%.3e,%.3e,%.3e).\n",
             vx[i],vy[i],vz[i]);
   }

}

double square_of_distance ( double x1, double y1, double z1,
                            double x2, double y2, double z2 )
{
   double dx = x1 - x2;
   double dy = y1 - y2;
   double dz = z1 - z2;

   return dx*dx + dy*dy + dz*dz;
}

void gforce ( int n, double m[n], double x[n], double y[n], double z[n],
              double ax[n], double ay[n], double az[n] )
{
   int i,j;
   double r2,r,a;
   const double G = 1.0;   /* just to keep the numbers simple ... */

   for(i=0; i<n; i++)
   {
      ax[i] = 0.0;
      ay[i] = 0.0;
      az[i] = 0.0;
      for(j=0; j<n; j++)
         if(i!=j)
         { 
            r2 = square_of_distance(x[i],y[i],z[i],x[j],y[j],z[j]);
            r = sqrt(r2);
	    a = G*m[j]/r2;             /* common to all components */
            ax[i] += a*(x[j]-x[i])/r;  /* components of direction of force */
	    ay[i] += a*(y[j]-y[i])/r;
            az[i] += a*(z[j]-z[i])/r;
         }
   }
}

void move ( int n, double x[n], double y[n], double z[n],
            double vx[n], double vy[n], double vz[n],
            double ax[n], double ay[n], double az[n], double dt )
{
   int i;

   for(i=0; i<n; i++)
   {
      vx[i] += ax[i]*dt;   /* update velocities */
      vy[i] += ay[i]*dt;
      vz[i] += az[i]*dt;
      x[i] += vx[i]*dt;    /* update position */
      y[i] += vy[i]*dt; 
      z[i] += vz[i]*dt;
   }
}