typedef struct
{ 
	char 	name[20];
	double 	mass;
	double 	x_pos,y_pos,z_pos;
	double 	x_vel,y_vel,z_vel;
} ObjectT;

/*
** Gravity_Proc : The Intergalactic Warping Calculating Routine
**                for FAST calculating of gravitational motions
**
** USAGE:-ObjectT object[] = array of type ObjectT, it contains
**        the necassary information for the calculations:
**
**        	typedef struct
**        	{ 
**			char 	naam[20]; <- not used in the routine
**			double 	mass;
**			double 	x_pos,y_pos,z_pos;
**			double 	x_vel,y_vel,z_vel;
**		} ObjectT;
**
**       -int objects = the funtion needs to know how much objects
**	  there are in the system.
**
**	 -double G = The gravity-constant = 6.67e-11
**
**	 -double dT = The 'infinitisemale' time-interval
**
**	 -int iterations = the number of iterations to be done
**			   previously always 1, but this routine is
**			   really FAST, so I warned you.
*/ 

void Gravity_Proc(ObjectT object[], int objects, double G, double dT, int iterations)
{
/* ObjectT type
**	{ 
**	char 	naam[20];
**	double 	mass;
**	double 	x_pos,y_pos,z_pos;
**	double 	x_vel,y_vel,z_vel;
**	}
*/

int i;
int first,second;
double r2,r;
double f,fx,fy,fz;
double dvxm,dvym,dvzm;

for (i=0; i<iterations; i++)
   {
	/*
	** The purpose of the way of using this 2 loops
	** is to avoid double calculations: e.v.
	** Calculating the force from the earth to the sun
	** and then again from the sun to the earth.
	*/ 

	for (first=0; first<objects; first++)
	{							/* +3 */
	for (second=first+1; second<objects; second++)
	{							/* +4 */

	/* Distance */
	r2= (object[first].x_pos - object[second].x_pos)
	   *(object[first].x_pos - object[second].x_pos)	
	   +(object[first].y_pos - object[second].y_pos)
	   *(object[first].y_pos - object[second].y_pos)
	   +(object[first].z_pos - object[second].z_pos)
	   *(object[first].z_pos - object[second].z_pos);
	r=sqrt(r2);

	/* Forces */
	f=G*(object[first].mass * object[second].mass)/r2;
	fx=f*(object[first].x_pos - object[second].x_pos)/r;
	fy=f*(object[first].y_pos - object[second].y_pos)/r;
	fz=f*(object[first].z_pos - object[second].z_pos)/r;

	/* Accelarations */
	dvxm=fx*dT;
	dvym=fy*dT;
	dvzm=fz*dT;
	object[first].x_vel -= (dvxm/object[first].mass);
	object[first].y_vel -= (dvym/object[first].mass);
	object[first].z_vel -= (dvzm/object[first].mass);
	object[second].x_vel += (dvxm/object[second].mass);
	object[second].y_vel += (dvym/object[second].mass);
	object[second].z_vel += (dvzm/object[second].mass);
	}							/* -4 */

	}							/* -3 */

	/* 
	** Don't make the 'spiral' miscalculation.
	** The velocity and the displacement has to be 
	** calculated for all objects separately.
	** (A spiral is a circular enclosing orbit)
	*/

	/* Displacement */
	for (first=0; first<objects; first++)
	{
	object[first].x_pos +=object[first].x_vel*dT;
	object[first].y_pos +=object[first].y_vel*dT;
	object[first].z_pos +=object[first].z_vel*dT;
	}
   }
}