THE CODE TO BE MODIFIED:

#include

#include

#include

#include

using namespace std;

int gravder(double x[], double t, double par[], double deriv[]);

int rk4(double q[], double t, double dt,

int (*derfunc)(double q[], double t, double par[], double deriv[]),

double par[]);

int gravder(double q[], double t, double par[], double deriv[]){

// deriv = dq/dt -- for the first two elements it is

// the velocity and the second two it's the force.

// par is a 1-d array where par[0] is GM

// t is not really needed but we include it so that

// we can use this with our general rk4 function.

deriv[0] = q[2]; //derivative of x = v_x

deriv[1] = q[3];

// the next two lines are the only different things

// for N2L problems.

deriv[2] = - par[0] * q[0]/pow(sqrt( q[0]* q[0] + q[1]*q[1]),3);

deriv[3] = - par[0] * q[1]/pow(sqrt( q[0]* q[0] + q[1]*q[1]),3);

return 0;

}

int rk4(double q[], double t, double dt,

int (*derfunc)(double q[], double t, double par[], double deriv[]),

double par[]){

double t12;

double qtmp[4];

double k1[4], k2[4], k3[4], k4[4];

int c1,c2,c3,c4;

/* First eval k1 */

c1 = derfunc(q,t,par,k1); // k1 = derfunc(q,t)

t12 = t + dt/2.;

for(int n=0;n

qtmp[n] = q[n] + k1[n]*dt/2.;

}

/* Now k2 */

c2 = derfunc(qtmp,t12,par,k2);

for(int n=0;n

qtmp[n] = q[n] + k2[n]*dt/2.;

}

/* Now k3 */

c3 = derfunc(qtmp,t12,par,k3);

for(int n=0;n

qtmp[n] = q[n] + k3[n]*dt;

}

t12 = t + dt;

/* Now k4 */

// c4 = derfunc(qtmp,t + dt,par,k4);

c4 = derfunc(qtmp,t12,par,k4);

for(int n=0;n

q[n] += (k1[n] + 2.*k2[n] + 2.*k3[n] + k4[n])*dt/6.0;

}

return c1+c2+c3+c4;

}

int main(int argc, char *argv[]){

double q[4]; // defines an array of length 4

// all arrays in C++ start counting at 0

double dt = 0.001, t = 0.0;

double tol = 1e-4;

double pi = acos(-1.0);

double GM = 4. * pi * pi;

double par[1];// array of parameters

int rval;

int nsteps= 30000; /* Should give three orbits*/

double energy;

ofstream fp;

double ell;

//int result = -1;

par[0] = GM;

// for(int i=0;i

// cout

// }

// return -2;

if(argc

cout

cout

cout

return -1;

//also fine if result were declared: return result;

}

else{

// atof(bbb) converts the string bbb to a floating-point number

// atoi(iii) converts the string iii to an integer

q[0] = atof(argv[1]); // 1.0;// x

q[1] = atof(argv[2]); // 0.0;// y

q[2] = atof(argv[3]); // 0.0;// x vel

q[3] = atof(argv[4]); // sqrt(GM);// y vel

}

fp.open("orbit_rk4.dat");

double kinetic = 0.5*(q[2]*q[2] + q[3]*q[3]);

double potential = - GM/sqrt(q[0]*q[0] + q[1]*q[1]);

energy = kinetic + potential;

ell = q[0]*q[3] - q[1]*q[2];

cout

fp

for(int i = 0 ; i

t = i*dt;

rval = rk4(q, t, dt, gravder, par);

kinetic = 0.5*(q[2]*q[2] + q[3]*q[3]);

potential = - GM/sqrt(q[0]*q[0] + q[1]*q[1]);

ell = q[0]*q[3] - q[1]*q[2];

energy = kinetic + potential;

// cout

fp

}

fp.close();

/*Plot orbit*/

FILE *gnuplot = popen("gnuplot", "w");

fprintf(gnuplot,"set out 'orbit_rk4.ps' ");

fprintf(gnuplot, "set term post land ");

fprintf(gnuplot, "plot 'orbit_rk4.dat' u 2:3 w l ");

fprintf(gnuplot, "plot 'orbit_rk4.dat' u 1:7 w l ");

pclose(gnuplot);

/*Plot energy*/

gnuplot = popen("gnuplot", "w");

fprintf(gnuplot,"set out 'orbit_energy_rk4.ps' ");

fprintf(gnuplot, "set term post land ");

fprintf(gnuplot, "plot 'orbit_rk4.dat' u 1:4 w l, 'orbit_rk4.dat' u 1:5 w l, 'orbit_rk4.dat' u 1:6 w l ");

pclose(gnuplot);

system((char *)"ps2pdf orbit_energy_rk4.ps");

system((char *)"ps2pdf orbit_rk4.ps");

return rval; //rval here is always 0.

}

kepler.cppl Modify our program orbit.rk4.cpp so that it simulates the orbit of a planet about the Sun and graphically displays the orbit. Allow for user-specified values of the initial distance from the sun and initial velocity using command-line arguments. Your program should also numerically (a) determine the period of the planet; (b) measure the semimajor axis and semiminor axis; (c) calculate the eccentricity of the orbit using Eq. (1) in problem 1 (d) calculate the eccentricity of the orbit using 62 and compare with part (2c) (e) verify Kepler's third law