1. Limited time only! Sign up for a free 30min personal tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

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
    2016 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.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted



Similar Discussions: Velocity Verlet applied to solar system with C++
  1. C++ Free fall velocity (Replies: 9)

Loading...