# Homework Help: Velocity Verlet applied to solar system with C++

Tags:
1. Nov 22, 2017

### Rebecca1990

1. The problem statement, all variables and given/known data
Hello, I am attempting to use the velocity Verlet algorithm of integration applied to the solar system in c++. My goal is be able to use the outputted position components in a plot to see if the trajectory of my object is elliptical/parabolic/hyperbolic resulting from the gravitational interaction with the Sun. I have coded it using c++. I then load the output file into matlab and plot. I did an example that mimicked a falling object that interacted with Earth, with a dt = 0.1. It worked fine when plotting. I attempted using x and y velocity components I found for Halley's comet using dt = 3.154e+7s (one year). When I load the data, and plot the x and y components of the position, I thought I would get an elliptical shaped object.Instead, I am getting a straight line that is decreasing as x increases. Can someone help me?

2. Relevant equations

3. The attempt at a solution
Code (C):
[/B]
#include <string>
#include <iostream>
#include <vector>
#include <cmath>
#include <fstream>
#include <cassert>

const double G = 6.67e-11;  // N*(m/kg)^2
//double dt = 0.1; // seconds for falling
double dt = 3.154e+7; // comet
class CelestialObject
{
private:
std::string m_name;
double m_mass;
double m_positionx;
double m_positiony;
double m_velocityx;
double m_velocityy;

public:
CelestialObject(std::string name, double mass, double* position, double* velocity) // constructor
{
m_name = name;
m_mass = mass;
m_positionx = position[0];
m_positiony = position[1];
m_velocityx = velocity[0];
m_velocityy = velocity[1];
}

// setters and getters

void setPosition(double* position)
{
m_positionx = position[0];
m_positiony = position[1];
}
void setVelocity(double* velocity)
{
m_velocityx = velocity[0];
m_velocityy = velocity[1];
}
double getMass() const
{
return m_mass;
}
double getPositionX() const
{
return m_positionx;
}
double getPositionY() const
{
return m_positiony;
}
double getVelocityX() const
{
return m_velocityx;
}
double getVelocityY() const
{
return m_velocityy;
}

friend std::eek:stream& operator<<(std::eek:stream& out, const CelestialObject& object)
{
return out << object.getPositionX() << '\t' << object.getPositionY() << '\t' << object.getVelocityX() << '\t' << object.getVelocityY() << std::endl;
}
};

// returns square of distance between objects
double distanceSquared(const CelestialObject &a, const CelestialObject &b)
{
// distance squared is (dy^2 + dx^2)
double rSquared = pow(a.getPositionY()-b.getPositionY(),2) + pow(a.getPositionX()-b.getPositionX(),2);
return rSquared;
}
// returns distance between objects
double distance(const CelestialObject &a, const CelestialObject &b)
{
// distance is (dy^2 + dx^2)^(1/2)
double r = sqrt(pow(a.getPositionY()-b.getPositionY(),2) + pow(a.getPositionX()-b.getPositionX(),2));
return r;
}

// returns magnitude of the force between the objects
double force(const CelestialObject &a, const CelestialObject &b)
{
//  F=(G * m1 * m1)/(r^2) in the direction a->b and b->a
double forceGrav = G*a.getMass()*b.getMass()/distanceSquared(a, b);
return forceGrav;
}

// returns the angle from a to b
double angle(const CelestialObject &a, const CelestialObject &b)
{
double angle = atan2f(b.getPositionY()-a.getPositionY(),b.getPositionX()-a.getPositionX());
return angle;
}

// Velocity Verlet algorithm

void updatePosition(CelestialObject &a, CelestialObject &b )
{
double F = force(a,b);
double theta = angle(a,b);
double accelerationA = F/a.getMass();
double accelerationB = -F/b.getMass();

// now that we have the acceleration of both objects, update positions
// x = x +v *dt + a*dt*dt/2 = x + dt * (v + a*dt/2)
double newPositionA[2] = {a.getPositionX() + dt * (a.getVelocityX() + accelerationA*((b.getPositionX() - a.getPositionX())/distance(a,b))*dt/2), a.getPositionY() + dt * (a.getVelocityY() + accelerationA*((b.getPositionY() - a.getPositionY())/distance(a,b))*dt/2)};
a.setPosition(newPositionA);
double newPositionB[2] = {b.getPositionX() + dt * (b.getVelocityX() + accelerationB*((b.getPositionX() - a.getPositionX())/distance(a,b))*dt/2), b.getPositionY() + dt * (b.getVelocityY() + accelerationB*((b.getPositionY() - a.getPositionY())/distance(a,b))*dt/2)};
b.setPosition(newPositionB);

// get new acceleration a'
F = force(a,b);
double thetaNew = angle(a,b);
double newAccelerationA = F/a.getMass();
double newAccelerationB = -F/b.getMass();

// update velocities
// v = v + (a + a')*dt/2
double newVelocityA[2] = {a.getVelocityX() + (accelerationA*((b.getPositionX() - a.getPositionX())/distance(a,b)) + newAccelerationA*((b.getPositionX() - a.getPositionX())/distance(a,b)))*dt/2, a.getVelocityY() + (accelerationA*((b.getPositionY() - a.getPositionY())/distance(a,b)) + newAccelerationA*((b.getPositionY() - a.getPositionY())/distance(a,b)))*dt/2};
a.setVelocity(newVelocityA);
double newVelocityB[2] = {b.getVelocityX() + (accelerationB*((b.getPositionX() - a.getPositionX())/distance(a,b)) + newAccelerationB*((b.getPositionX() - a.getPositionX())/distance(a,b)))*dt/2, b.getVelocityY() + (accelerationB*((b.getPositionY() - a.getPositionY())/distance(a,b)) + newAccelerationB*((b.getPositionY() - a.getPositionY())/distance(a,b)))*dt/2};
b.setVelocity(newVelocityB);
}

int main()
{
// example of halley's comet
double positionHalley[2] = {48696201283, -6.8734238e+10};
double positionSun[2] = {5.6378e7,5.6378e7}; // in metres, take origin at centre so equal component on both axes
double velocityHalley[2] = {-9.096111, -6.916686};
double velocitySun[2] = {0,0}; // may have to switch
CelestialObject halley("halley", 2.2e+14, positionHalley, velocityHalley);
CelestialObject sun("sun", 1.99e+30, positionSun, velocitySun);
std::cout << "Initial values:\n" << sun << halley;

std::eek:fstream write_output("halley2.txt");
assert (write_output.is_open());

int t;
for (t=0; t < 80; t++)
{
write_output << dt*t << '\t' << halley;
updatePosition(sun, halley);
}
std::cout << "Final values at t = " << dt*t << " seconds:\n";
std::cout << "Sun: " << sun << "\n";
std::cout << "Halley: "<< halley << "\n";
write_output.close();
/*
// example of mass falling off building file output
double positionMass[2] = {0, 370};
double positionEarth[2] = {0, -6.378e6};
double velocityMass[2] = {0, 0};
double velocityEarth[2] = {0,0};
CelestialObject mass("mass", 70, positionMass, velocityMass);
CelestialObject earth("earth", 5.97e+24, positionEarth, velocityEarth);

std::eek:fstream write_output("falling.txt");
assert (write_output.is_open());

int t;
for (t=0; mass.getPositionY() > 0; t++)
{
write_output << dt*t << '\t' << mass;
updatePosition(mass, earth);
}
std::cout << "Final values at t = " << dt*t << " seconds:\n";
std::cout << "Earth: " << earth << "\n";
std::cout << "Mass: "<< mass << "\n";
*/

return 0;
}

Last edited by a moderator: Nov 22, 2017
2. Nov 23, 2017

### Staff: Mentor

Did you try smaller steps?
Did you check the size of all the relevant numbers? A straight line could mean you underestimate the strength of gravity.

Share this great discussion with others via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted