- #1
Rebecca1990
- 1
- 0
Homework Statement
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?
Homework Equations
The Attempt at a Solution
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: