Entry 5296

// Fourth-order Runge-Kutta integrator for solving N-body newtonian problem
//
// This is only a first approach to the problem.

#include <cmath>
#include <cstdio>

/*
  Class vector3: tridimensional vectors and their operators
*/

struct vector3 {
	double x, y, z;
	vector3() {
		x = y = z = 0;
	}
	vector3(double X, double Y, double Z) {
		x = X;
		y = Y;
		z = Z;
	}
	vector3 operator+ (vector3 a) {
		vector3 sum(x + a.x, y + a.y, z + a.z);
		return sum;
	}
	vector3 &operator+= (vector3 a) {
		x += a.x;
		y += a.y;
		z += a.z;
		return *this;
	}
	vector3 operator- () {
		vector3 op(-x, -y, -z);
		return op;
	}
	vector3 operator- (vector3 a) {
		vector3 sum(x - a.x, y - a.y, z - a.z);
		return sum;
	}
	vector3 &operator-= (vector3 a) {
		x -= a.x;
		y -= a.y;
		z -= a.z;
		return *this;
	}
	vector3 operator* (double k) {
		vector3 res(x * k, y * k, z * k);
		return res;
	}
	vector3 &operator*= (double k) {
		x *= k;
		y *= k;
		z *= k;
		return *this;
	}
	vector3 operator/ (double k) {
		vector3 res(x / k, y / k, z / k);
		return res;
	}
	friend double abs(vector3);
	friend double cubeabs(vector3);
	friend double dot(vector3, vector3);
	friend vector3 cross(vector3, vector3);
};

inline vector3 makevec(double a) {
	vector3 ret(a, a, a);
	return ret;
}

inline vector3 makevec(double a, double b, double c) {
	vector3 ret(a, b, c);
	return ret;
}

inline double abs(vector3 a) {
	return sqrt(a.x * a.x + a.y * a.y + a.z * a.z);
}

inline double dot(vector3 a, vector3 b) {
	return (a.x * b.x + a.y * b.y + a.z * b.z);
}

inline double cubeabs(vector3 a) {
	double sqabs = a.x * a.x + a.y * a.y + a.z * a.z;
	return (sqabs * sqrt(sqabs));
}

inline vector3 cross(vector3 a, vector3 b) {
	vector3 ret(a.y*b.z - a.z*b.y, a.z*b.x - a.x*b.z, a.x*b.y - a.y*b.x);
	return ret;
}

/*
  class solvebody: numerically integrate N-body problem
*/
class solvebody {
	/* Parameters of the integrator */
	int n;       // number of bodies
	char steps;  // number of calls to step mod 2 (error estimation)
	double time; // time
	double tol;  // step size [TODO automatic stepsize with tolerance]

	/* Body data */
	double *gm;  // product of the mass by the gravitational constant
	vector3 *r0; // initial positions for the bodies 0..n-1 
	vector3 *v0; // initial velocities for the bodies 0..n-1

	/* State information */
	vector3 *k1, *k2, *k3, *k4; // Constants used in Runge-Kutta method
	vector3 *epsk2, *epsk3, *epsk4; // Constants used in error estimation

	public:
	solvebody() {
		n = steps = 0;
		tol = time = 0.0;
		gm = NULL;
		r0 = v0 = k1 = k2 = k3 = k4 = epsk2 = epsk3 = epsk4 = NULL;
	}

	// Constructor
	/* Argument list:
		(double)itime    initial time
		(int)n           number of bodies
		(double [n])gm   product of the mass by the gravitational constant
		(vector3 [n])r0  initial positions
		(vector3 [n])v0  initial velocities
		(double)tol      tolerance (see documentation) [TODO explain when done]
		Negative values for gm have special meanings. If gm = -0.5, the body is
		not affected by gravitational forces. If gm = -n, the body is fixed with
		relation to the body numbered (n - 1).
    */
	solvebody(double, int, double *, vector3 *, vector3 *, double = 0.01);

	// Step the Runge-Kutta integrator
	/* Argument list:
		(vector3 **)r	(*r)[0..n-1] is set to the radius vector list
		(vector3 **)v	(*v)[0..n-1] is set to the velocity list
		(vector3 **)eps	(*eps)[0..n-1] are the estimated errors in 
	   Returns: the time for which the state vector has been computed
	*/
	double step();

	/* Clears the error estimate and step count */
	double clearerr();

	/* Destructor method */
	~solvebody() {
		delete [] k1;
	}
	friend int main(int, char **);
};

// Constructor
solvebody::solvebody(double TIME, int N, double *GM, vector3 *R0, vector3 *V0,
	                  double TOL) {
	time = TIME;
	n = N;
	gm = GM;
	r0 = R0;
	v0 = V0;
	tol = TOL;
	k1 = new vector3[8*N];
	k2    = k1 +     N;
	k3    = k1 + 2 * N;
	k4    = k1 + 3 * N;
	epsk2 = k1 + 4 * N;
	epsk3 = k1 + 5 * N;
	epsk4 = k1 + 6 * N;
}

// Step the Runge-Kutta integrator
double solvebody::step( ) {
	double &h = tol;
	// Zero coefficients	
	for(int i = 0; i != n; ++i)
		k1[i] = k2[i] = k3[i] = k4[i] = makevec(0);
	// Initial point
	for(int i = 0; i != n; ++i) {
		for(int j = i + 1; j != n; ++j) {
			vector3 relpos = r0[j] - r0[i];
			double crelpos = cubeabs(relpos);
			if(crelpos != 0.0) {
				k1[i] += relpos * gm[j] / crelpos;
				k1[j] += -relpos * gm[i] / crelpos;
			}
		}
		k1[i] *= h;
	}
	// Double step size and calculate
	if(!steps) {
		for(int i = 0; i != n; ++i)
			epsk2[i] = epsk3[i] = epsk4[i] = makevec(0);
	}
	// Midpoint, first evaluation
	for(int i = 0; i != n; ++i) {
		for(int j = i + 1; j != n; ++j) {
			vector3 relpos = (v0[j] - v0[i]) * h / 2 + r0[j] - r0[i];
			double crelpos = cubeabs(relpos);
			if(crelpos != 0.0) {
				k2[i] += relpos * gm[j] / crelpos;
				k2[j] += -relpos * gm[i] / crelpos;
			}
		}
		k2[i] *= h;
	}
	// Midpoint, second evaluation
	for(int i = 0; i != n; ++i) {
		for(int j = i + 1; j != n; ++j) {
			vector3 relpos = (v0[j] - v0[i] + (k1[j] - k1[i]) / 2) * h / 2 +
			                   r0[j] - r0[i];
			double crelpos = cubeabs(relpos);
			if(crelpos != 0.0) {
				k3[i] += relpos * gm[j] / crelpos;
				k3[j] += -relpos * gm[i] / crelpos;
			}
		}
		k3[i] *= h;
	}
	// Final point
	for(int i = 0; i != n; ++i) {
		for(int j = i + 1; j != n; ++j) {
			vector3 relpos = (v0[j] - v0[i] + (k2[j] - k2[i]) / 2) * h +
			                   r0[j] - r0[i];
			double crelpos = cubeabs(relpos);
			if(crelpos != 0.0) {
				k4[i] += relpos * gm[j] / crelpos;
				k4[j] += -relpos * gm[i] / crelpos;
			}
		}
		k4[i] *= h;
	}
	// Advance state vectors
	for(int i = 0; i != n; ++i) {
		r0[i] += (v0[i] + (k1[i] + k2[i] + k3[i]) / 6) * h;
		v0[i] += (k1[i] + k4[i]) / 6 + (k2[i] + k3[i]) / 3;
	}
	steps = (steps + 1) % 2;
	return (time += h);
}

int main(int argc, char *argv[]) {
	if(argc < 2) {
		fprintf(stderr, "usage: newton data [output]\n");
		return 0;
	}
	FILE *fin = fopen(argv[1], "r");
	if(fin == NULL) {
		fprintf(stderr, "Could not open file %s for input\n", argv[1]);
		return 0;
	}
	FILE *fout = stdout;
	if(argc > 2) {
		fout = fopen(argv[2], "w+");
		if(fout == NULL) {
			fprintf(stderr,
			        "Could not open file %s for output; using stdout instead\n",
			        argv[2]);
			fout = stdout;
		}
	}
	double g, step, time;
	int n;
	fscanf(fin, "%lf\n%lf\n%lf\n%d\n", &g, &step, &time, &n);
	fprintf(stderr, "Gravitational constant = %1.4E\n", g);
	fprintf(stderr, "Step size              = %lf\n", step);
	fprintf(stderr, "Initial time           = %1.9E\n", time);
	fprintf(stderr, "Number of bodies       = %d\n",  n);
	if(n < 2 || n > 320) {
		fprintf(stderr, "Number of bodies outside expected range\n");
		return 0;
	}
	char *dummy = new char[256];
	double *gm = new double[n];
	vector3 *pos = new vector3[n];
	vector3 *vel = new vector3[n];
	for(int i = 0; i != n; ++i) {
		fscanf(fin, "%s\n", dummy);
		fscanf(fin, "%lf\n%lf;%lf;%lf\n%lf;%lf;%lf|\n", &gm[i], &pos[i].x,
		       &pos[i].y, &pos[i].z, &vel[i].x, &vel[i].y, &vel[i].z);
		gm[i] *= g;
	}
	solvebody solver(time, n, gm, pos, vel, step);
	for(;;) {
		fprintf(fout, "%1.8E\n", solver.time);
		for(int i = 0; i != n; ++i) {
			fprintf(fout, "%d 0; %+1.11E; %+1.11E; %+1.11E\n",i, solver.r0[i].x,
			        solver.r0[i].y, solver.r0[i].z);
			fprintf(fout, "%d 1; %+1.11E; %+1.11E; %+1.11E\n",i, solver.v0[i].x,
			        solver.v0[i].y, solver.v0[i].z);
		}
		solver.step();
	}
	delete [] gm;
	delete [] pos;
	delete [] vel;
	delete [] dummy;
	fclose(fin);
	if(fout != stdout)
		fclose(fout);
	return 0;
}