/******************************************************************** * Two bodies problem - with sfml * *********************************************************************/ #include #include #include #include #include #include using namespace std; const double G = 6.67428e-11; class Vector2d { public: double x, y; Vector2d(double xin = 0, double yin = 0) : x(xin), y(yin) {} void operator +=(Vector2d& v) { //defenition of new operators += and *= x += v.x; y += v.y; } Vector2d operator*(double scalar) { return Vector2d(scalar * x, scalar * y); } }; class Body { private: string name; Vector2d position; Vector2d velocity; double mass; //sf::CircleShape; public: Body(string n, Vector2d p = { 0,0 }, Vector2d v = { 0,0 }, double m = 0) : name(n), position(p), velocity(v), mass(m) {} void setPosition(Vector2d& p) { position += p; } void setVelocity(Vector2d& v) { velocity += v; } Vector2d getPosition() { return position; } Vector2d getVelocity() { return velocity; } string getName() { return name; } double getMass() { return mass; } }; ostream& operator<<(ostream& os, Body& b) { return os << b.getPosition().x << "," << b.getPosition().y << "\n"; } /******************************************************************* function that calculates r *******************************************************************/ double Distance(Body a, Body b) { return sqrt((a.getPosition().x - b.getPosition().x) * (a.getPosition().x - b.getPosition().x) //sqrt( (x2-x1)^2 + (y2-y1)^2 + (a.getPosition().y - b.getPosition().y) * (a.getPosition().y - b.getPosition().y)); } /****************************************************************** Function that updates the acceleration, velocity and position ********************************************************************/ void update(std::vector& bodies, double dT) { // vector that holds acceleration vector accelerations(bodies.size()); //cout << "bodies.size() = " << bodies.size() << endl; // loop that computes total acceleration of each body. for (int i = 0; i < int(bodies.size()); i++) { for (int j = 0; j < int(bodies.size()); j++) { if (i != j) { double r = Distance(bodies[i], bodies[j]); accelerations[i] = Vector2d(((G * bodies[j].getMass()) //(G*mass)/R * (position 1 - positon2) * (bodies[j].getPosition().x - bodies[i].getPosition().x)) / (r * r * r), //r^3 ((G * bodies[j].getMass()) * (bodies[j].getPosition().y - bodies[i].getPosition().y)) / (r * r * r)); } } } // loop to make sure acceleration is correct /** for (int i = 0; i < int(bodies.size()); ++i) { cout << "acc[" << i << "]=(" << accelerations[i].x << ", " << accelerations[i].y << ")\n"; } cin.get(); */ // vectors that hold the velocities and positions of the bodies vector changeInVelocities(bodies.size()); vector changeInPositions(bodies.size()); // loop for updating the new velocities and positions for (int i = 0; i < int(bodies.size()); ++i) { changeInVelocities[i] = accelerations[i] * dT; bodies[i].setVelocity(changeInVelocities[i]); changeInPositions[i] = bodies[i].getVelocity() * dT; bodies[i].setPosition(changeInPositions[i]); } } //************************************************************* /* Cesar, That's just great! I added a few file streams to put the data into comma separated lists so I could plot them out in Excel and they look convincingly orbital...I added a second month to see how stable the orbits are. The moon looks pretty stable, but the Earth is a little wobbly. Also, the sidereal lunar month is 27 d 7 h 43 min 11.6 s 27 d * 24 hr/d * 3600 sec/hr = 2,332,800 sec 7 hr * 3600 sec/hr = 25,200 sec 43 min * 60 sec/min = 2,580 sec + 11.6 sec = 2,360,591.6 sec and your moon is completing the orbit in only about 22 days. The center of mass is 7.35e22 kg *3.85e8 m/(7.35e22 kg + 5.98e24)~4.76e6 m from the center of the earth. As for the speed of the Earth, it's traveling a distance 2pi*4.76e6 m /2,360,591.6 sec~12.6m/sec */ int main() { //GH: Added ofstreams for gathering data: ofstream moon("moon.csv"), earth("earth.csv"); moon << "Moon\n x,y\n"; earth << "Earth\n x,y\n"; // name ,position,velocity,mass Body Earth("Earth", Vector2d(-4.672e6, 0), Vector2d(0, -12.435), 5.972e24); Body Moon("Moon", Vector2d(3.7973e8, 0), Vector2d(0, 1010.7), 7.34767e22); vector bodies{ Earth,Moon }; double dT = 1.0; //DeltaT for functions int dT2 = 1; //for the loop that prints out the output while (dT2 <= 4721184) //GH: doubled the total time 2360592) { update(bodies, dT); dT2++; //loop counter if (dT2 % 86400 == 0) { //PRINTING OUT EVERY DAY cout << "Time Elapsed: " << dT2 / 86400 << " Days" << endl; for (int i = 0; i < int(bodies.size()); ++i) { cout << bodies[i]; // .getName() //<< " will be at the new points: ( " //<< bodies[i].getPosition().x << " , " //<< bodies[i].getPosition().y << " )\n" << endl; } //GH: added these two lines for data files moon << bodies[1];// .getPosition().x << ',' << bodies[1].getPosition().y << endl; earth << bodies[0];// .getPosition().x << ',' << bodies[0].getPosition().y << endl; } } }