/* 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"); } }