Commit 74758d57 authored by s-jdoyle4's avatar s-jdoyle4
Browse files

Answer to 1 and 2 of HW implemented from git reposit orbit_bw.C

parent f11d2cb5
// orbit - Program to compute the orbit of a comet.
#include <iostream.h>
#include <fstream.h>
#include <assert.h>
#include <math.h>
void main() {
//* Set initial position and velocity of the comet.
double r0, v0;
cout << "Enter initial radial distance (AU): "; cin >> r0;
cout << "Enter initial tangential velocity (AU/yr): "; cin >> v0;
double r[2], v[2], accel[2];
r[0] = r0; r[1] = 0; v[0] = 0; v[1] = v0;
//* Set physical parameters (mass, G*M)
const double pi = 3.141592654;
double GM = 4*pi*pi; // Grav. const. * Mass of Sun (au^3/yr^2)
double param[1]; param[0] = GM;
double mass = 1.; // Mass of comet
double F = .01*4*pi*pi/r0/r0; // const outside force in one direction one percent of initial gravitational force
double time = 0;
//* Loop over desired number of steps using specified
// numerical method.
cout << "Enter number of steps: ";
int nStep; cin >> nStep;
cout << "Enter time step (yr): ";
double tau; cin >> tau;
cout << "Choose a numerical method:" << endl;
cout << "1) Euler, 2) Euler-Cromer, 3) Verlet " << endl;
int method; cin >> method;
double *rplot, *thplot, *tplot, *kinetic, *potential, *totalE; // Plotting variables
rplot = new double [nStep]; thplot = new double [nStep];
tplot = new double [nStep];
kinetic = new double [nStep]; potential = new double [nStep];
totalE = new double [nStep];
int iStep;
for( iStep=0; iStep < nStep; iStep++ ) {
//* Record position and energy for plotting.
double normR = sqrt( r[0]*r[0] + r[1]*r[1] );
double normV = sqrt( v[0]*v[0] + v[1]*v[1] );
rplot[iStep] = normR; // Record position for plotting
thplot[iStep] = atan2(r[1],r[0]);
tplot[iStep] = time;
kinetic[iStep] = 0.5*mass*normV*normV; // Record energies
potential[iStep] = -GM*mass/normR+F*r0*mass; // With the additional constant potential from the force
totalE[iStep] = kinetic[iStep]+potential[iStep];
//* Calculate new position and velocity using desired method.
if( method == 1 ) {
accel[0] = -GM*r[0]/(normR*normR*normR)+.01*GM*r0; // constant force effects the acceleration of the object
accel[1] = -GM*r[1]/(normR*normR*normR)+.01*GM*r0;
r[0] += tau*v[0]; // Euler step
r[1] += tau*v[1];
v[0] += tau*accel[0];
v[1] += tau*accel[1];
time += tau;
}
else if( method == 2 ) {
accel[0] =-GM*r[0]/(normR*normR*normR)+.01*GM*r0; // constant force effects the acceleration of the object
accel[1] =-GM*r[1]/(normR*normR*normR)+.01*GM*r0;
v[0] += tau*accel[0];
v[1] += tau*accel[1];
r[0] += tau*v[0]; // Euler-Cromer step
r[1] += tau*v[1];
time += tau;
}
else if( method == 3 ) {
double accel[3]; double r[3]; r[0]=r0; r[1]=0; r[2]=0; v[1]=v0; double oldr;
accel[0] =-GM*r[0]/(normR*normR*normR)+.01*GM*r0;
accel[1] =-GM*r[1]/(normR*normR*normR)+.01*GM*r0;
accel[2]=-GM*r[2]/(normR*normR*normR)+.01*GM*r0;
oldr = r[1];
r[0] = r[1] - tau*v[1]+tau*tau*accel[1];
r[1] = 2*r[1] - r[0]+tau*tau*accel[1];
r[2] = 2*r[1] - oldr+tau*tau*accel[2];
v[0] = .5*r[1]/tau +.5*r[0]/tau;
v[1] = .5*r[2]/tau +.5*oldr/tau;
time += tau;
}
}
//loop through data to find min and max for plotting purposes
double min, max;
min = kinetic[0]; max = kinetic[0];
for( iStep=0; iStep < nStep; iStep++ ) {
if(kinetic[iStep] > max)
max = kinetic[iStep];
if(potential[iStep] > max)
max = potential[iStep];
if(totalE[iStep] > max)
max = totalE[iStep];
if(kinetic[iStep] < min)
min = kinetic[iStep];
if(potential[iStep] < min)
min = potential[iStep];
if(totalE[iStep] < min)
min = totalE[iStep];
}
//Set up the canvas and divide it into two parts
TCanvas *MyC = new TCanvas("MyC", "lorenz", 1);
MyC->Divide(2,1);
MyC->cd(1);
//Graph the energy vs time
TGraph *graph = new TGraph(nStep, tplot, kinetic);
graph->SetName("gr1");
graph->Draw("A*");
graph->SetMinimum(min); //set graph's range
graph->SetMaximum(max);
TGraph *graph2 = new TGraph(nStep, tplot, potential);
graph2->SetName("gr2");
graph2->SetMarkerColor(2);
graph2->Draw("*");
TGraph *graph3 = new TGraph(nStep, tplot, totalE);
graph3->SetName("gr3");
graph3->SetMarkerColor(3);
graph3->Draw("*");
//*Provide title and axis labels
graph->SetTitle("Energies for Orbiting Comet");
graph->GetXaxis()->SetTitle("Time(yr)");
graph->GetXaxis()->CenterTitle();
graph->GetYaxis()->SetTitle("Energy (M*AU^2/yr^2");
graph->GetYaxis()->CenterTitle();
// draw the legend
TLegend *legend=new TLegend(0.6,0.65,0.88,0.85);
legend->SetTextFont(72);
legend->SetTextSize(0.04);
legend->AddEntry(gr1,"Kinetic","P");
legend->AddEntry(gr2,"Potential","P");
legend->AddEntry(gr3,"Total","P");
legend->Draw();
//switch to second half of canvas
MyC->cd(2);
// convert polar coordinates to cartesian
double *xplot, *yplot;
xplot = new double[nStep]; yplot = new double[nStep];
for( iStep=0; iStep < nStep; iStep++ ) {
xplot[iStep] = rplot[iStep]*cos(thplot[iStep]);
yplot[iStep] = rplot[iStep]*sin(thplot[iStep]);
}
//plot trajectory of comet in cartesian coordinates, with Sun at origin
TGraph *graph4 = new TGraph(nStep, xplot, yplot);
graph4->SetName("gr4");
graph4->Draw("A*");
//*Provide title and axis labels
graph4->SetTitle("Trajectory of Comet");
graph4->GetXaxis()->SetTitle("x-axis(AU)");
graph4->GetXaxis()->CenterTitle();
graph4->GetYaxis()->SetTitle("y-axis(AU)");
graph4->GetYaxis()->CenterTitle();
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment