物理月球地球系统，月球轨道问题

问题描述 投票：0回答：1

我正在尝试使用 LeapFrog 获取地球和月球绕太阳运行的轨道。我已经设法获得了地球的正确轨道，但在月球的情况下，我无法使其绕地球运行，同时绕太阳运行，我认为这可能是初始条件的问题，但我不确定，这是我写的代码：

#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;
}

这就是我所获得的： Orbits Obtained

我希望让月球绕地球运行，同时它们都绕太阳运行。

c++ leapfrog-integration
1个回答
0
投票

您的代码有很多问题......或者您的期望。

以下内容均已被其他人在评论中提及：

(1) 每个（地球）轨道需要更多步数。毕竟，一年有 12 个阴历月（大约）。

(2) py 的表达式是错误的，需要月球和地球的相对位移。

(3) 您的初始条件略有错误（下面的代码中列出了替代版本。

但最后还没有提到你的期望。

(4) 如果你在同一张图表上绘制地球和月球的路径，那么，由于地球和太阳之间的距离大约是地球和月球之间距离的 500 倍，那么你需要一个显示器和比我更好的视力才能看到轨道差异。

在下面的代码中，我添加了另一个输出文件“LFrelative.dat”（抱歉！），它仅显示月球相对于地球的位置。绘制此图给出了距离我们最近的天体的轨道 - 请注意，比例尺完全不同。

#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;
}

太阳轨道尺度上的地球和月球：

月球相对于地球的轨道 - 注意比例尺上的距离：

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.