/*************************************************************/
/*Gravity Simulator*******************************************/
/***********************Scallops, inc.************************/
/************************Conroy Lee***************************/
/**************************7/31/91****************************/
/*************************************************************/

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

#define GRAV_MAX_NUM_OBJECTS 42       /* max number of objects in the system */
#define STRING_MAX_LENGTH 255

typedef struct _SystemRec *System;
typedef struct _ObjectRec *Object;

struct _ObjectRec                 /******an object in the system*****/
{
  double m;                          /* mass of an object                    */
  double *x_array;                   /* array of x-coordinates of object     */
  double *y_array;                   /* array of y-coordinates               */
  double *vx_array;                  /* array of x-component of velocity     */
  double *vy_array;                  /* array of y-component                 */
  double *ax_array;                  /* array of x-component of acceleration */
  double *ay_array;                  /* array of y-component                 */
};

struct _SystemRec                /*******a system******************/
{
  Object planets;
  double G;                           /* Gravitational Constant              */
  int n;                              /* number of objects                   */
  double dt;                          /* length of each time step            */
  double t;                           /* length of time to watch system      */
  int i_steps;                        /* number of time steps                */
};

int data_entry(System*);                               /* entering data      */
void init_planet(Object,int,int);                      /*   for a planet     */
void sum_accel(Object,struct _ObjectRec,int);
double dist(struct _ObjectRec,struct _ObjectRec,int);  /* dist bet. planets  */
void move_object(Object,double,int);                   /* move a planet      */
void run_simulation(struct _SystemRec);                /* run the model      */
void print_simulation(System);                         /* print results      */
void cdsPrintData(System);

main(int argc, char **argv)
{
  System mysys;
  data_entry(&mysys);
  run_simulation(*mysys);
  cdsPrintData(mysys);
}

void cdsPrintData(System thesys)
{
  int i,j;
  printf("ready\n");
  for (j=0;j<thesys->n;j++)
    {
      printf("$GravityPlotter insert line plot-%d -xydata {", j);
      for (i=0;i<thesys->i_steps-1;i++)
	printf("%lf %lf ", 
	       thesys->planets[j].x_array[i]/1e+9,
	       thesys->planets[j].y_array[i]/1e+9);
      printf("}\n");
    }
  printf("ready\n");
}

int data_entry(System *sysPtrPtr)
/* Scan variables in this order:          */
/*    G    Gravitational constant         */
/*    t    time to run system             */
/*    dt   length of time for time step   */
/*    n    number of objects in the sys.  */
{
  int i,n,i_steps;
  System sysPtr=*sysPtrPtr;

  sysPtr= (System) malloc (sizeof(struct _SystemRec));

  *sysPtrPtr=sysPtr;
  printf("#Gravitational constant?\n");
    scanf("%lf",&sysPtr->G);
    if (sysPtr->G <0)
      exit(1);
  printf("#Time to run system? (sec)\n");
    scanf("%lf",&sysPtr->t);
    if (sysPtr->t < 0)
      exit(1);
  printf("#Length of each time step: (sec)\n");
    scanf("%lf",&sysPtr->dt);
    if (sysPtr->dt <= 0)
      exit(1);
  sysPtr->i_steps = sysPtr->t / sysPtr->dt;
  i_steps= sysPtr->i_steps;
  printf("#How many objects are there?\n");
    scanf("%d",&sysPtr->n);
    n= sysPtr->n;
    if (n<0 || n>GRAV_MAX_NUM_OBJECTS)
      exit(1);
  sysPtr->planets= (Object) malloc (n* sizeof(struct _ObjectRec));
  for (i=0;i<n; i++)
    init_planet(sysPtr->planets+i,i,i_steps);
}

void init_planet(Object planet,int planet_num,int i_steps)
/* Scan variables in this order:              */
/*    x_array[i]: x position at time-step i   */
/*    y_array[i]: y position at time-step i   */
/*    vx_array[i]: x velocity at time-step i  */
/*    vy_array[i]: y velocity at time-step i  */

{
  planet->x_array= (double*) malloc (i_steps * (sizeof (double)));
  planet->y_array= (double*) malloc (i_steps * (sizeof (double)));
  planet->vx_array= (double*) malloc (i_steps * (sizeof (double)));
  planet->vy_array= (double*) malloc (i_steps * (sizeof (double)));
  planet->ax_array= (double*) malloc (i_steps * (sizeof (double)));
  planet->ay_array= (double*) malloc (i_steps * (sizeof (double)));
  printf("Mass of planet %d:\n",planet_num);
    scanf("%lf",&planet->m);
  printf("Initial X position of planet %d:\n",planet_num);
    scanf("%lf",&planet->x_array[0]);
  printf("Initial Y position of planet %d:\n",planet_num);
    scanf("%lf",&planet->y_array[0]);
  printf("Initial X velocity of planet %d:\n",planet_num);
    scanf("%lf",&planet->vx_array[0]);
  printf("Initial Y velocity of planet %d:\n",planet_num);
    scanf("%lf",&planet->vy_array[0]);
}

void run_simulation(struct _SystemRec MySystem)
{
  int i,j,k;
  double accelx,accely;
  i=1;
  while (i<MySystem.i_steps)                /* cycle through time steps */
    {
      MySystem.t += MySystem.dt;
      for (j=0;j<MySystem.n;j++)            /* cycle through planets    */
	{
	  i--;
	  MySystem.planets[j].ax_array[i]= 0;
	  MySystem.planets[j].ay_array[i]= 0;
	  for (k=0;k<j;k++)
	    sum_accel(&MySystem.planets[j],MySystem.planets[k],i);
	  for (k=j+1;k<MySystem.n;k++)
	    sum_accel(&MySystem.planets[j],MySystem.planets[k],i);
	  MySystem.planets[j].ax_array[i]*= MySystem.G;
	  MySystem.planets[j].ay_array[i]*= MySystem.G;
/*	  printf("%d\ttime-step %d\t%lf\t%lf\n",j,i,MySystem.planets[j].x_array[i],MySystem.planets[j].y_array[i]);	  */

	  i++;
	  move_object(&MySystem.planets[j],MySystem.dt,i);
	
	}
      i++;
    }
}

void sum_accel(Object Planet1, struct _ObjectRec Planet2, int i)
{
  double r=dist(*Planet1,Planet2,i);
  if (r!=0)
    {
      Planet1->ax_array[i]+= Planet2.m/r/r * (Planet2.x_array[i]-Planet1->x_array[i])/ r;
      Planet1->ay_array[i]+= Planet2.m/r/r * (Planet2.y_array[i]-Planet1->y_array[i])/ r;
    }
}

double dist(struct _ObjectRec Obj1, struct _ObjectRec Obj2, int i)
{
  return(sqrt((Obj1.x_array[i]-Obj2.x_array[i])*(Obj1.x_array[i]-Obj2.x_array[i]) + (Obj1.y_array[i]-Obj2.y_array[i])*(Obj1.y_array[i]-Obj2.y_array[i])));
}

void move_object(Object planet, double dt, int i)
{
  planet->vx_array[i]= planet->vx_array[i-1]+ planet->ax_array[i-1]*dt;
  planet->vy_array[i]= planet->vy_array[i-1]+ planet->ay_array[i-1]*dt;
  planet->x_array[i]= planet->x_array[i-1]+ planet->vx_array[i]*dt;
  planet->y_array[i]= planet->y_array[i-1]+ planet->vy_array[i]*dt;
}
/*
void print_simulation(System mySys)
{
  int i,j;
  printf("Planet\tStep\tX Position\t Y Position\n");
  for (i=0;i<mySys.i_steps;i++)
    for (j=0;j<mySys.n;j++)
      
  exit(1);
}
*/