Velocity Verlet applied to solar system with C++

[Rebecca1990]

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

[/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::count << "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::count << "Final values at t = " << dt*t << " seconds:\n";
   std::count << "Sun: " << sun << "\n";
   std::count << "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::count << "Final values at t = " << dt*t << " seconds:\n";
   std::count << "Earth: " << Earth << "\n";
   std::count << "Mass: "<< mass << "\n";
*/
   return 0;
}

[/B]

 

[mfb]

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.