/* 12.010 Homework 3, Question 3.
N-body orbit integration problem
Program to solve the N-body gravitational problem using Runge-Kutta
numerical integration.
Arguments for the program are passed through runsting:
% HW03_03
*/
#include
#include
#include
#include
#include
#include "NBodyC.h"
/* Prototype definitions */
void read_nbody() ;
void report_ICS(double , double [][max_body],double [][max_body], char *) ;
void int_step(double [][max_body], double [][max_body], double [][max_body], double [][max_body], double , double) ;
void get_accel(double , double [][max_body], double [][max_body], double [][max_body]) ;
void eval_step(double , double [][max_body], double [][max_body], double , double *, int) ;
void animate(double , double , double [][max_body]);
/* Start Main program */
int main(argc, argv)
int argc ; // Number of arguments
char *argv[] ; // Pointer to arguments
{
double t ; // Current time (days)
double step ; // Current step size (days)
double new_step ; // Updated step size that will yield desired accuracy
double poss[2][max_body], vels[2][max_body] ; /* Position and velocity at start
of integration step */
double posc[2][max_body], velc[2][max_body] ; /* Position and velocity at end
of integration step when done as one step */
double posa[2][max_body], vela[2][max_body] ; /* Position and velocity at end
of first half/step integration */
double posb[2][max_body], velb[2][max_body] ; /* Position and velocity at end
of second half/step integration (posb should = posc
within tolerance.*/
double minstep ; // Minimum step size needed during integaration
int sc ; // Counter for steps size calculaiton so that we don't get stuck
int step_not_OK ; // set non-zero while step size is still being evaluated
int no_step_increase ; // Set non-zero to stop step from inceasing.
int ierr ; // Number of entries read with sscanf (0 means none)
int i ; // Loop counter
/* Startby decoding the runstring arguments */
if( argc > 1 ) {
strcpy(infile,argv[1]);
} else {
printf("HW02_03:N-body problem. Correct runstring:\n");
printf("% HW02_03 \n");
printf("Default duration is 515 days, default tolerance 1.e-6\n");
exit(-1);
}
/* See if other arguments are passed */
if( argc > 2 ) {
ierr = sscanf(argv[2],"%lf", &dur);
if ( ierr < 1 ) {
printf("HW02_03:N-body problem. Correct runstring:\n");
printf("% HW02_03 \n");
printf("Default duration is 515 days, default tolerance 1.e-6\n");
printf("Error reading the duration\n");
exit(-2);
}
} else {
dur = 515.0 ;
}
/* Check tolerance See if passed */
if( argc > 3 ) {
ierr = sscanf(argv[3],"%lf", &tol);
if ( ierr < 1 ) {
printf("HW02_03:N-body problem. Correct runstring:\n");
printf("% HW02_03 \n");
printf("Default duration is 515 days, default tolerance 1.e-6\n");
printf("Error reading the tolerance\n");
exit(-3);
}
} else {
tol = 1.e-6 ;
}
printf("File %s\n",infile);
printf("Run %6.1f days, %10.2e tol\n",dur, tol);
/* Now read in the data about the bodies to be integrated */
read_nbody() ;
/* Start the integration */
step = 1.0 ; minstep = step ;
t = 0.0;
/* Initialize starting position and velocity */
for ( i = 0 ; i < num_body ; i++ ) {
poss[0][i] = posi[0][i] ; poss[1][i] = posi[1][i];
vels[0][i] = veli[0][i] ; vels[1][i] = veli[1][i];
}
/* Report initial conditions */
report_ICS(t, poss, vels, "Initial Conditions");
animate(t,step, poss);
/* Now start the integrator */
while ( t < dur ) {
/* Nake sure that the next step will nor overshoot the end
of the time (this can happen because are are varying the
step size during the run. */
no_step_increase = 0;
if( t + step > dur ) {
step = dur - t ;
no_step_increase = 1;
}
/* Iterate on step size to get the value needed to maintain the
requested tolerance */
step_not_OK = 1; sc = 0 ;
while ( step_not_OK != 0 && sc < 15 ) {
sc++;
int_step( posc, velc, poss, vels, step, t);
// Now take two steps with half size
int_step( posa, vela, poss, vels, step/2, t );
int_step( posb, velb, posa, vela, step/2, t+step/2 );
/* Compare the two results and see if the meet tolerances or if
results too good can we increase step */
eval_step( t, posc, posb, step, &new_step, sc) ;
// Make sure we have not increased step at wrong time
if( no_step_increase != 0 && new_step > step ) {
new_step = step;
}
// If the new_step is the same as step then we are done.
if( new_step == step ) {
step_not_OK = 0;
} else {
step = new_step ;
}
}
// save smallest step size
if( step < minstep ) {minstep = step;}
// OK step hase been correctly taken, update state and time
t = t + step ;
for ( i = 0 ; i < num_body ; i++ ) {
poss[0][i] = posb[0][i] ; poss[1][i] = posb[1][i];
vels[0][i] = velb[0][i] ; vels[1][i] = velb[1][i];
}
/* Use more man's graphics with vt100 sequences */
animate(t, step, poss) ;
}
/* Output the final values */
report_ICS(t, poss, vels, "Final conditions");
printf("Smallest step size needed %9.6f days\n",minstep);
printf("Closest approach distance was %14.5e km Bodies %d and %d\n",closest, close_body[0], close_body[1]);
}
void read_nbody() {
/* Function to read the initial conditions of the bodies
Routine to read the input data file.
File format lines are
Mass (kg) Xpos (km) Ypos (km) Xvel (km/s) Yvel (km/s)
*/
int n ; // Counter for body number
int jerr ; // number of entries read from data line
char line[128] ; // line from file
FILE *fp, *fopen() ; // File pointer
/* Open the file and readlines */
fp = fopen(infile,"r");
if ( fp == NULL ) {
printf("\nError opening file %s\n",infile);
exit(-1);
}
/* Now read file */
line[0] = 1;
n = 0 ; // Counter for number of bodies
while ( fgets(line, 128, fp) != NULL && n < max_body-1 ) {
if( line[0] != 0 && line[0] == ' ' ) {
// decode the line
jerr = sscanf(line,"%lf %lf %lf %lf %lf", &mass[n], &posi[0][n], &posi[1][n], &veli[0][n], &veli[1][n]);
if( jerr < 5 ) {
// Not enough values read
printf("Only %d values read in line \n%s\n",jerr,line);
exit(-1);
}
n++;
}
}
if( n == max_body ) {
printf("HW03_03: Maximum number of bodies %d reached\n",n);
}
num_body = n;
}
void report_ICS(double t, double posr[][max_body], double velr[][max_body], char *type) {
/* Routine to report status of bodies at time t ond of type type */
int i ; // Loop variable
printf("+ 12.010 HW 03 Q 03: %s At time %14.5f days for %d bodies\n",type, t, num_body);
printf("Body Mass (kg) PosX (km) PosY (km) VelX (km/s) VelY (km/s)\n");
for ( i = 0 ; i < num_body ; i++ ) {
printf("%2d %18.6e %18.7e %18.7e %18.7e %18.7e \n",i+1, mass[i], posr[0][i], posr[1][i], velr[0][i], velr[1][i]);
}
}
void int_step( double pose[][max_body], double vele[][max_body], double poss[][max_body], double vels[][max_body], double step, double t) {
/* Routine to integrate 1 time step of duration step in days */
double SS ; // Step size in seconds
double P[2][max_body], V[2][max_body], A[2][max_body] ; /* Positon, velcity and acceleration at
inter[2][max_body]mediate points */
double k1[2][max_body], k2[2][max_body], k3[2][max_body], k4[2][max_body]; /* Intermediate
integration points in Runge-Kutta integration.*/
int i ; // Loop counter
// Step in seconds
SS = step*86400.0 ;
// Get the accelerations of all bodies at start
get_accel(t, poss, vels, A) ;
for( i = 0 ; i < num_body ; i ++ ) {
k1[0][i]= SS*A[0][i];
k1[1][i]= SS*A[1][i];
// Step system to mid-point
P[0][i] = poss[0][i] + SS*vels[0][i]/2 + SS*k1[0][i]/8;
P[1][i] = poss[1][i] + SS*vels[1][i]/2 + SS*k1[1][i]/8;
V[0][i] = vels[0][i] + k1[0][i]/2 ;
V[1][i] = vels[1][i] + k1[1][i]/2 ;
}
// Get acceleration at this new point
get_accel(t+step/2, P, V, A);
for( i = 0 ; i < num_body ; i ++ ) {
k2[0][i]= SS*A[0][i];
k2[1][i]= SS*A[1][i];
// Step system to mid-point
P[0][i] = poss[0][i] + SS*vels[0][i]/2 + SS*k1[0][i]/8;
P[1][i] = poss[1][i] + SS*vels[1][i]/2 + SS*k1[1][i]/8;
V[0][i] = vels[0][i] + k2[0][i]/2 ;
V[1][i] = vels[1][i] + k2[1][i]/2 ;
}
// Next value at mid point
get_accel(t+step/2, P, V, A);
for( i = 0 ; i < num_body ; i ++ ) {
k3[0][i]= SS*A[0][i];
k3[1][i]= SS*A[1][i];
// Step system
P[0][i] = poss[0][i] + SS*vels[0][i] + SS*k3[0][i]/2;
P[1][i] = poss[1][i] + SS*vels[1][i] + SS*k3[1][i]/2;
V[0][i] = vels[0][i] + k3[0][i] ;
V[1][i] = vels[1][i] + k3[1][i] ;
}
// Now do final end point
get_accel(t+step, P, V, A);
for( i = 0 ; i < num_body ; i ++ ) {
k4[0][i]= SS*A[0][i];
k4[1][i]= SS*A[1][i];
// Step system to final point
pose[0][i] = poss[0][i] + SS*(vels[0][i] + (k1[0][i]+k2[0][i]+k3[0][i])/6);
pose[1][i] = poss[1][i] + SS*(vels[1][i] + (k1[1][i]+k2[1][i]+k3[1][i])/6);
vele[0][i] = vels[0][i] + (k1[0][i]+2*k2[0][i]+2*k3[0][i]+k4[0][i])/6;
vele[1][i] = vels[1][i] + (k1[1][i]+2*k2[1][i]+2*k3[1][i]+k4[1][i])/6;
}
}
void get_accel( double t, double P[][max_body], double V[][max_body], double A[][max_body]) {
/* Routine to compute accelerations */
int i, j ; // Counters over the bodies.
double dP[2][max_body] ; // Position difference for each body
double R[max_body] ; // Distances bewteen current body and others
/* Loop over each body */
for ( i = 0; i < num_body ; i++ ) {
A[0][i] = 0.0 ; A[1][i] = 0.0 ;
for ( j = 0 ; j < num_body ; j++ ) {
if ( i != j ) {
dP[0][j] = P[0][j] - P[0][i] ;
dP[1][j] = P[1][j] - P[1][i] ;
R[j] = sqrt(dP[0][j]*dP[0][j]+dP[1][j]*dP[1][j]);
/* Now get vector sum of accelerations */
A[0][i] = A[0][i] + G_univ*mass[j]*dP[0][j]/(R[j]*R[j]*R[j]);
A[1][i] = A[1][i] + G_univ*mass[j]*dP[1][j]/(R[j]*R[j]*R[j]);
//printf("ACC %d %d %e %e %e\n",i,j,dP[0][j],R[j],A[0][i]);
//printf("G Mass %e %e \n",G_univ,mass[j] );
}
}
}
}
void eval_step( double t, double posc[][max_body], double posb[][max_body], double step, double *new_step, int sc) {
/* Routine to see if step size meetds tolerances */
double dPe[2][max_body] ; // Difference in position estimates
double R ; // distance betweeb bodies
double Rmin ; // Smallest difference
double err ; // Position error in km
double minerr = 1.e20 , maxerr = 0.0 ; // Min and max relative error
int i, j ; // Loop counters
/* Get the difference in positions and see how low */
for ( i = 0 ; i < num_body ; i++ ) {
dPe[0][i] = posc[0][i] - posb[0][i];
dPe[1][i] = posc[1][i] - posb[1][i];
}
/* See if meets tolerance */
for ( i = 0 ; i < num_body ; i++ ) {
err = sqrt(dPe[0][i]*dPe[0][i]+dPe[1][i]*dPe[1][i]);
Rmin = 1.e20 ;
for ( j = 0 ; j < num_body ; j++ ) {
if( j != i ) {
R = sqrt((posb[0][j]-posb[0][i])*(posb[0][j]-posb[0][i]) +
(posb[1][j]-posb[1][i])*(posb[1][j]-posb[1][i])) ;
if( R < Rmin ) {Rmin = R;}
if( Rmin < closest ) {
closest = Rmin;
close_body[0] = i+1; close_body[1] = j+1;
}
}
}
/* Now see if this smallest or largest error */
if( err/Rmin < minerr) {minerr = err/Rmin ; }
if( err/Rmin > maxerr) {maxerr = err/Rmin ; }
}
/* See if we have meet tolerance */
if( sc > 11 ) {
printf("'UpateStep Time %12.6f %3d Size $9.6f Errs %13.3e %13.3e Tol %13.3e\n", t, sc, step, maxerr, minerr, tol);
}
*new_step = step ;
if( maxerr < tol/40 && step < 1.0 ) {
*new_step = step*2;
}
if( maxerr > tol && step > 1.0/65536.0 ) {
*new_step = step/2;
}
}
void animate( double t, double step, double pos[][max_body]) {
/* Routine to animate the changes in position using vt100 escape sequences */
double maxx, maxy, minx,miny ; // Min and max coordintes in initial state
double maxa, mina ; // Max and min of both x and y
static double scalex, scaley ; // Scales to convert x and y to screen coordinates
static double offx, offy ; // Offsets for origin in X and Y
int SWX = 85, SHY = 45 ; // Screen width and height
int IX, IY ; // Coorinates in screen coordinatea
int i, j ; // Loop counter
static char sym[max_body] ; // Symbols for each body
char screen[86][45] ; // 45-lines of 85-characters to output (extra char for null)
static char mv45[6] ; // Escape sequence to move up 51 lines (extra char for null)
/* Initialize when t = 0; */
if( t == 0 ) {
/* First call; so don't move cursor back up and compute size we need
Get o or . based on mass */
for ( i = 0 ; i < num_body ; i ++ ) {
if( mass[i] > 1.e30 ) {
sym[i] = 'O' ;
} else if( mass[i] > 1.e28 ) {
sym[i] = 'o' ;
} else {
sym[i] = '.' ;
}
}
/* Now get the scale */
maxx = pos[0][0]; minx = pos[0][0];
maxy = pos[1][0]; miny = pos[1][0];
for ( i = 1 ; i < num_body ; i ++ ) {
maxx = fmax(maxx, pos[0][i]); minx = fmin(minx, pos[0][i]);
maxy = fmax(maxy, pos[1][i]); miny = fmin(miny, pos[1][i]);
}
maxa = fmax(maxx, maxy); mina = fmin(minx, miny);
/* Set Scale in X and Y */
scalex = SWX/(maxa-mina) * 0.95;
scaley = SHY/(maxa-mina) * 0.95;
offx = minx * 1.05 ;
offy = (miny - 1.5e8 ) * 1.05;
sprintf(mv45,"%c[47A",(char) 27) ;
} else {
/* Move the cursor back 52-lines */
printf("%s\n",mv45);
}
/* Output for all calls: */
printf("TIME %13.6f days, Step %10.6f days",t,step);
/* OK form up the output lines */
for( i = 0 ; i < 45 ; i++ ) {
for ( j = 0 ; j < 85 ; j ++ ) {
screen[j][i] = ' ';
}
screen[86][i] = 0;
}
for ( i = 0 ; i < num_body ; i ++ ) {
IX = (int) ((pos[0][i]-offx)*scalex+0.5) - 1 ;
IY = 45-(int) ((pos[1][i]-offy)*scaley+0.5) -1 ;
if( IX >= 0 && IX < 85 &&
IY >= 0 && IY < 45 ) {
if( screen[IX][IY] == ' ' ) {
screen[IX][IY] = sym[i];
}
}
}
/* write lines to screen */
for ( i = 0 ; i < 45 ; i++ ) {
for ( j = 0 ; j < 85 ; j++ ) {
printf("%c", screen[j][i]);
}
printf("\n");
}
}