Physical Moon Earth System, problems with Moon's orbit
I am trying to obtain the trajectories of the orbit of earth and moon around the Sun using LeapFrog. I have managed to obtain the correct orbit of the Earth but in the moon case I haven't been able to make it orbit the Earth while both orbit the Sun, I think it might be a problem of initial conditions but I am not sure, this is the code that I have written:
#include <iostream>
#include <fstream>
#include <cmath>
using namespace std;
double Grav = 6.6726e-11; // m^3/kg/s^2 (Universal gravitational constant)
double Msun = 1.989e30; // kg
double Mearth = 5.973e24; // kg
double Mmoon = 7.349e22; // kg
double Dse = 1.496e11; // distance Sun-Earth in meters
double Dem = 3.844e8; // distance Earth-Moon in meters
int nsteps = 100; // Timesteps per orbit
void LeapFrog(int n)
{
// Initial Conditions
double x0 = Dse;
double y0 = 0;
double vx0 = 0;
double vy0 = sqrt(Grav * Msun / Dse);
double px0 = Mearth * vx0;
double py0 = Mearth * vy0;
// Initial Conditions for Moon
double xm0 = Dse + Dem;
double ym0 = 0;
double vxm0 = sqrt(Grav * Mearth / Dem); // La velocidad en x es opuesta a la velocidad orbital de la Tierra
double vym0 = vy0; // La velocidad en y es igual a la velocidad orbital de la Tierra
double pxm0 = Mmoon * vxm0;
double pym0 = Mmoon * vym0;
// Timestep
double dt = (2 * 3.141592654 * Dse / vy0) / nsteps; // fraction of a year
// Create a file for saving the coordinates of Earth
ofstream earthFile("LFEarth.dat");
// Create a file for saving the coordinates of Moon
ofstream moonFile("LFMoon.dat");
earthFile << x0 << " " << y0 << '\n';
moonFile << xm0 << " " << ym0 << '\n';
for (int i = 0; i < n; i++)
{
// Update Earth position
double x = x0 + dt * px0 / (2 * Mearth);
double y = y0 + dt * py0 / (2 * Mearth);
// Update Moon position
double xm = xm0 + dt * pxm0 / (2 * Mmoon);
double ym = ym0 + dt * pym0 / (2 * Mmoon);
// Calculate distances to the Sun
double r_earth_sun = sqrt(x * x + y * y);
double r_moon_sun = sqrt(xm * xm + ym * ym);
double r_moon_earth = sqrt((xm-x)*(xm-x) + (ym-y)*(ym-y));
// Update Earth momentum
double px = px0 - dt * Grav * Msun * Mearth * x / pow(r_earth_sun, 3) - dt * Grav * Mmoon * Mearth * (x - xm) / pow(r_moon_earth, 3);
double py = py0 - dt * Grav * Msun * Mearth * y / pow(r_earth_sun, 3) - dt * Grav * Mmoon * Mearth * y / pow(r_moon_earth, 3);
// Update Moon momentum
double pxm = pxm0 - dt * Grav * (Msun * Mmoon) * (xm) / pow(r_moon_sun, 3) - dt * Grav * (Mearth * Mmoon) * (xm - x) / pow(r_moon_earth, 3);
double pym = pym0 - dt * Grav * (Msun * Mmoon) * (ym) / pow(r_moon_sun, 3) - dt * Grav * (Mearth * Mmoon) * (ym - y) / pow(r_moon_earth, 3);
// Another Step for Earth position
double x2 = x + dt * px / (2 * Mearth);
double y2 = y + dt * py / (2 * Mearth);
// Another Step for Moon position
double xm2 = xm + dt * pxm / (2 * Mmoon);
double ym2 = ym + dt * pym / (2 * Mmoon);
earthFile << x2 << " " << y2 << std::endl;
moonFile << xm2 << " " << ym2 << std::endl;
// Initialise the values for the following step
x0 = x2;
y0 = y2;
xm0 = xm2;
ym0 = ym2;
px0 = px;
py0 = py;
pxm0 = pxm;
pym0 = pym;
}
// Close the files
earthFile.close();
moonFile.close();
}
int main()
{
LeapFrog(nsteps * 10); // Ten orbits
return 0;
}
And this is what I have obtained:
I expected to get the moon to orbit the Earth while they both orbit the Sun.
Solution
There are a number of things wrong with your code ... and your expectations.
The following have all been alluded to by other people in comments:
(1) You need more steps per (earth) orbit. After all, there are 12 lunar months in a year (roughly).
(2) The expression for py is wrong, where it needs relative displacement of moon and earth.
(3) Your initial conditions are slightly wrong (alternative version sketched as a comment in the code below.)
But not mentioned yet are your expectations in the end.
(4) If you plot BOTH earth and moon paths on the same graph, then, since the distance between earth and sun is about 500 times that between earth and moon, you would need both a monitor and eyesight better than I possess to see a difference in orbits.
In the code below I add another output file "LFrelative.dat" (sorry!) which shows just the moon's position relative to the earth. Plotting this gives the orbit of our nearest celestial body - note that the scale is completely different.
#include <iostream>
#include <fstream>
#include <cmath>
using namespace std;
const double Grav = 6.6726e-11; // m^3/kg/s^2 (Universal gravitational constant)
const double Msun = 1.989e30; // kg
const double Mearth = 5.973e24; // kg
const double Mmoon = 7.349e22; // kg
const double Dse = 1.496e11; // distance Sun-Earth in meters
const double Dem = 3.844e8; // distance Earth-Moon in meters
const int nsteps = 1000; // Timesteps per orbit
const double PI = 4 * atan( 1.0 );
void LeapFrog(int n)
{
// Initial Conditions
double x0 = Dse;
double y0 = 0;
double vx0 = 0;
double vy0 = sqrt(Grav * Msun / Dse);
double px0 = Mearth * vx0;
double py0 = Mearth * vy0;
// Initial Conditions for Moon
double xm0 = x0; // <--M
double ym0 = Dem; // |
double vxm0 = sqrt(Grav * Mearth / Dem); // S---------------------------E
double vym0 = vy0;
double pxm0 = Mmoon * vxm0;
double pym0 = Mmoon * vym0;
// Timestep
double dt = (2 * PI * Dse / vy0) / nsteps; // fraction of a year
// Create a file for saving the coordinates of Earth
ofstream earthFile("LFEarth.dat");
// Create a file for saving the coordinates of Moon
ofstream moonFile("LFMoon.dat");
ofstream relativeFile("LFRelative.dat");
earthFile << x0 << " " << y0 << '\n';
moonFile << xm0 << " " << ym0 << '\n';
relativeFile << xm0 - x0 << " " << ym0 - y0 << '\n';
for (int i = 0; i < n; i++)
{
// Update Earth position
double x = x0 + dt * px0 / (2 * Mearth);
double y = y0 + dt * py0 / (2 * Mearth);
// Update Moon position
double xm = xm0 + dt * pxm0 / (2 * Mmoon);
double ym = ym0 + dt * pym0 / (2 * Mmoon);
// Calculate distances to the Sun
double r_earth_sun = sqrt(x * x + y * y);
double r_moon_sun = sqrt(xm * xm + ym * ym);
double r_moon_earth = sqrt((xm-x)*(xm-x) + (ym-y)*(ym-y));
// Update Earth momentum
double px = px0 - dt * Grav * Msun * Mearth * x / pow(r_earth_sun, 3) - dt * Grav * Mmoon * Mearth * (x - xm) / pow(r_moon_earth, 3);
double py = py0 - dt * Grav * Msun * Mearth * y / pow(r_earth_sun, 3) - dt * Grav * Mmoon * Mearth * (y - ym) / pow(r_moon_earth, 3);
// Update Moon momentum
double pxm = pxm0 - dt * Grav * Msun * Mmoon * xm / pow(r_moon_sun, 3) - dt * Grav * Mearth * Mmoon * (xm - x) / pow(r_moon_earth, 3);
double pym = pym0 - dt * Grav * Msun * Mmoon * ym / pow(r_moon_sun, 3) - dt * Grav * Mearth * Mmoon * (ym - y) / pow(r_moon_earth, 3);
// Another Step for Earth position
double x2 = x + dt * px / (2 * Mearth);
double y2 = y + dt * py / (2 * Mearth);
// Another Step for Moon position
double xm2 = xm + dt * pxm / (2 * Mmoon);
double ym2 = ym + dt * pym / (2 * Mmoon);
earthFile << x2 << " " << y2 << std::endl;
moonFile << xm2 << " " << ym2 << std::endl;
relativeFile << xm2 - x2 << " " << ym2 - y2 << std::endl;
// Initialise the values for the following step
x0 = x2;
y0 = y2;
xm0 = xm2;
ym0 = ym2;
px0 = px;
py0 = py;
pxm0 = pxm;
pym0 = pym;
}
// Close the files
earthFile.close();
moonFile.close();
relativeFile.close();
}
int main()
{
LeapFrog(nsteps * 10); // Ten orbits
return 0;
}
Earth and moon on SOLAR ORBIT scale:
Moon ORBIT RELATIVE TO THE EARTH - note the difference in scale:
- Storing value infinity to a variable
- C compilers default behavior with constants
- Reading a 4 bytes date (big endian) from binary file
- -Werror=format: how can the compiler know
- check unused variable in C project
- When using clang on OSX, what controls the automatically added command line options?
- How to avoid using printf in a signal handler?
- Remove console window with SDL application
- gcc-arm-none-eabi 11.3 "is not implemented and will always fail"
- is this the right way to use Includes in C?
- Unexpected output from right rotate function
- How can I convert virtual address from user space to physical address?(in linux device driver)
- What is the difference between open and creat system call in c?
- Are all elements of the dirent structure open to use?
- Is there a way to utilize all XMM registers?
- In the C if-else statement, should the condition which is more likely to be true come first?
- Optimization for global and static variables
- Having trouble with char and scanf_s()
- Allegro error in DevC++
- Compress string using bits Operation in C
- What is the difference in atomic_load() and assignment?
- bind: invalid arguments AF_UNIX
- Explicit cast of nested pointers and const qualification
- Ignoring or redefining GCC Standard Predefined Macros
- Format strings safely when vsnprintf is not available
- How can I print to standard error in C with 'printf'?
- __STDC_LIB_EXT1__ availability in gcc and clang
- How to reassign void pointer pointing to a string to a new string in C
- memset() or value initialization to zero out a struct?
- Can quicksort be implemented in C without stack and recursion?