THE CODE TO BE MODIFIED: #include <iostream> #include <fstream> #include <cmath>
ID: 3757209 • Letter: T
Question
THE CODE TO BE MODIFIED:
#include <iostream>
#include <fstream>
#include <cmath>
#include <cstdlib>
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<4;n++) {
qtmp[n] = q[n] + k1[n]*dt/2.;
}
/* Now k2 */
c2 = derfunc(qtmp,t12,par,k2);
for(int n=0;n<4;n++) {
qtmp[n] = q[n] + k2[n]*dt/2.;
}
/* Now k3 */
c3 = derfunc(qtmp,t12,par,k3);
for(int n=0;n<4;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<4;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<5;i++){
// cout << "argv[" << i << "] = "<<argv[i] << endl;
// }
// return -2;
if(argc < 5){
cout << "Usage: " << argv[0] << " x0 y0 vx0 vy0 " << endl;
cout << " x0 and y0 are the initial x and y positions..."<<endl;
cout << "For circular orbit use 1 0 0 6.28" << endl;
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 << t << " " << q[0] << " " << q[1] << " " << energy << endl;
fp << t << " " << q[0] << " " << q[1] << " " << energy << " " << kinetic << " " << potential << " " << ell << endl;
for(int i = 0 ; i < nsteps ; 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 << t << " " << q[0] << " " << q[1] << " " << energy << endl;
fp << t << " " << q[0] << " " << q[1] << " " << energy << " " << kinetic << " " << potential << " " << ell << endl;
}
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.Explanation / Answer
#include "stdafx.h"
#include <stdint.h>
#include <conio.h>
#include <math.h>
#include <iostream>
#include <conio.h>
#include <vector>
#include <fstream>
#include <Windows.h>
#include <string>
#include <vector>
using namespace std;
double Ax(double t, double x, double vx)
{
double wd = 2./3.;
double A = 1.15;
double gamma = 0.5;
double accel = -gamma*vx - sin(x) + A*sin(wd*t);
return accel;
}
double Vx(double t, double x,double vx)
{
double ret;
ret = vx;
return vx;
}
void main()
{
//our integration loops are nested inside of this initial condtion loop
//for a range of allowable conditions y crossings at our plane will
//generated and recorded. Par for the course fists of ham will be used
//to generate a very dense output of both initial conditions and thus
//a dense sos plot. Plotting will be performed outside of the scope of this
//cpp program, in mathematica possibly if it plays nice and doesn't
//crash upon import.
//Opening Storage File for saving x vx phase space stuffs
ofstream plane,lambdak;
//plane.open("C:\Users\Uma Mobile\Dropbox\Spring 2013\380\hw6\Plots\phasedata.txt");
//lambdak.open("C:\Users\Uma Mobile\Dropbox\Spring 2013\380\hw6\Plots\lambdadata.txt");
plane.open("C:\Users\Uma\My Documents\My Dropbox\Spring 2013\380\hw6\Plots\phasedata.txt");
lambdak.open("C:\Users\Uma\My Documents\My Dropbox\Spring 2013\380\hw6\Plots\lambdadata.txt");
//Variable Declaration fiducial Orbit rk4
double x,vx;
double tf,t;
double pi = 3.1415926535898;
int count = 0, foolproof;
double k1x,k1vx;
double k2x,k2vx;
double k3x,k3vx;
double k4x,k4vx;
//Variable Declaration for Test Orbit rk4
double xy,vxy;
double d0;
double k1xy,k1vxy;
double k2xy, k2vxy;
double k3xy, k3vxy;
double k4xy, k4vxy;
int taui;
double dt;
double loggy,tempxy,tempvxy;
int counter;
vector<double> lambda;
vector<double> kk;
//Initiliazing Conditions for
d0 = 0.001;
x = 3;
vx = 0.0;
xy = x + .001;
vxy = vx;
taui = 1000;
double tau = (double) taui;
tf = 1000000.0;
t = 0.0;
foolproof = 0;
loggy = 0.0;
counter = 0;
dt = 0.01;
//First loop executes RK4 for trajectories over a Period of seconds = tf
while (t < tf)
{
t += dt;
foolproof += 1;
k1x = Vx(t, x, vx);
k1vx = Ax(t, x, vx);
k1xy = Vx(t, xy, vxy);
k1vxy = Ax(t, xy, vxy);
k2x = Vx(t, x + 0.5*dt*k1x, vx + 0.5*dt*k1vx);
k2vx = Ax(t, x + 0.5*dt*k1x, vx + 0.5*dt*k1vx);
k2xy = Vx(t, xy + 0.5*dt*k1xy, vxy + 0.5*dt*k1vxy);
k2vxy = Ax(t, xy + 0.5*dt*k1xy, vxy + 0.5*dt*k1vxy);
k3x = Vx(t, x + 0.5*dt*k2x, vx + 0.5*dt*k2vx);
k3vx = Ax(t, x + 0.5*dt*k2x, vx + 0.5*dt*k2vx);
k3xy = Vx(t, xy + 0.5*dt*k2xy, vxy + 0.5*dt*k2vxy);
k3vxy = Ax(t, xy + 0.5*dt*k2xy, vxy + 0.5*dt*k2vxy);
k4x = Vx(t, x + dt*k3x, vx + dt*k3vx);
k4vx = Ax(t, x + dt*k3x, vx + dt*k3vx);
k4xy = Vx(t, xy + dt*k3xy, vxy + dt*k3vxy);
k4vxy = Ax(t, xy + dt*k3xy, vxy + dt*k3vxy);
x = x + dt*(1.0/6.0)*(k1x + 2*k2x + 2*k3x + k4x);
vx = vx + dt*(1.0/6.0)*(k1vx + 2*k2vx + 2*k3vx + k4vx);
xy = xy + dt*(1.0/6.0)*(k1xy + 2*k2xy + 2*k3xy + k4xy);
vxy = vxy + dt*(1.0/6.0)*(k1vxy + 2*k2vxy + 2*k3vxy + k4vxy);
if (foolproof%(taui) == 0)
{
counter += 1;
kk.push_back(counter*d0*10);
tempxy = x + (d0*(xy -x))/sqrt((xy - x)*(xy - x)+ (vxy - vx)*(vxy - vx));
tempvxy = vx + (d0*(vxy- vx))/sqrt((vxy - vx)*(vxy - vx) + (xy - x)*(xy - x));
lambda.push_back(log((sqrt(((xy - x)*(xy - x))+ (vxy - vx)*(vxy - vx)))/d0));
plane << x << ' ' << vx << ' ' << xy << ' ' << vxy << endl;
xy = tempxy;
vxy = tempvxy;
}
}
double sum,out1,out2;
for (int i = 0; i < lambda.size(); i++)
{
sum = 0;
for ( int j = 0; j < i ; j++)
{
sum += (lambda[j])/(kk[i]*tau);
}
out1 = sum;
out2 = kk[i];
lambdak << out2 << ' ' << out1 << ' ' << endl;
}
plane.close();
lambdak.close();
cout << "Maths has computed"<< endl;
//program exits at keyboard input
//_getch();
}
Related Questions
Navigate
Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.