Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

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

  1. Nov 22, 2017 #1
    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. jcsd
  3. Nov 23, 2017 #2

    mfb

    User Avatar
    2017 Award

    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