fork download
copy 
  1. #include <iostream>
  2. #include <fstream>
  3. #include <iomanip>
  4. #include <cmath>
  5. #include <valarray>
  6.  
  7. using namespace std;
  8.  
  9. const double G_SI = 6.67430E-11; // Newton's constant in SI units
  10. const double M_EARTH = 5.9722E24; // Earth's mass in kg
  11. const double M_MOON = 7.342E22; // Lunar mass in kg
  12. const double xL = 384400000.0; // Lunar position in meters
  13. const double yL = 0.0;
  14.  
  15. // Right-hand side function of the system of ODEs
  16. std::valarray<double> rhs(double t, std::valarray<double> y) {
  17.     double r_Earth = sqrt(y[0] * y[0] + y[1] * y[1]);
  18.     double r3_Earth = r_Earth * r_Earth * r_Earth;
  19.  
  20.     double dx_L = y[0] - xL;
  21.     double dy_L = y[1] - yL;
  22.     double r_Moon = sqrt(dx_L * dx_L + dy_L * dy_L);
  23.     double r3_Moon = r_Moon * r_Moon * r_Moon;
  24.  
  25.     std::valarray<double> dydt(4);
  26.     dydt[0] = y[2]; // dx/dt = vx
  27.     dydt[1] = y[3]; // dy/dt = vy
  28.     dydt[2] = -G_SI * M_EARTH * y[0] / r3_Earth - G_SI * M_MOON * dx_L / r3_Moon; // dvx/dt
  29.     dydt[3] = -G_SI * M_EARTH * y[1] / r3_Earth - G_SI * M_MOON * dy_L / r3_Moon; // dvy/dt
  30.  
  31.     return dydt;
  32. }
  33.  
  34. // Euler's method
  35. std::valarray<double> euler(double t0, std::valarray<double> y0, double h, double tf) {
  36.     double t = t0;
  37.     std::valarray<double> y = y0;
  38.  
  39.     // Evolution loop for Euler's method
  40.     while (t < tf) {
  41.         y += h * rhs(t, y);
  42.         t += h;
  43.     }
  44.  
  45.     return y;
  46. }
  47.  
  48. // Runge-Kutta 4th order method
  49. std::valarray<double> rungeKutta4(double t0, std::valarray<double> y0, double h, double tf) {
  50.     double t = t0;
  51.     std::valarray<double> y = y0;
  52.  
  53.     // Evolution loop for RK4
  54.     while (t < tf) {
  55.         std::valarray<double> k1 = h * rhs(t, y);
  56.         std::valarray<double> k2 = h * rhs(t + 0.5 * h, y + 0.5 * k1);
  57.         std::valarray<double> k3 = h * rhs(t + 0.5 * h, y + 0.5 * k2);
  58.         std::valarray<double> k4 = h * rhs(t + h, y + k3);
  59.  
  60.         y += (k1 + 2.0 * k2 + 2.0 * k3 + k4) / 6.0;
  61.         t += h;
  62.     }
  63.  
  64.     return y;
  65. }
  66.  
  67. // Function to solve the ODEs using specified method and output to a file
  68. void solve_and_output(double t0, double t_end, double dt_output, const std::valarray<double>& y_init, const std::string& filename, char method) {
  69.     std::ofstream data_file(filename);
  70.     data_file << std::fixed << std::setprecision(5);
  71.     data_file << "Time\tX\tY\tVX\tVY" << std::endl;
  72.  
  73.     std::valarray<double> result = y_init;
  74.     double t = t0;
  75.     double dt = 10.0; // Integration step size
  76.  
  77.     data_file << t << "\t" << result[0] << "\t" << result[1] << "\t" << result[2] << "\t" << result[3] << std::endl;
  78.  
  79.     for (int i = 0; t < t_end; ++i) {
  80.         if (method == 'E') {
  81.             result = euler(t, result, dt, t + dt);
  82.         } else if (method == 'R') {
  83.             result = rungeKutta4(t, result, dt, t + dt);
  84.         } else {
  85.             std::cerr << "Invalid method choice." << std::endl;
  86.             return;
  87.         }
  88.         t += dt;
  89.  
  90.         // Output every dt_output seconds
  91.         if (i % static_cast<int>(dt_output / dt) == 0) {
  92.             data_file << t << "\t" << result[0] << "\t" << result[1] << "\t" << result[2] << "\t" << result[3] << std::endl;
  93.         }
  94.     }
  95.  
  96.     data_file.close();
  97. }
  98.  
  99. int main() {
  100.     // Initial conditions
  101.     double t0 = 0.0;
  102.     double x0 = 0.0;
  103.     double y0 = 26378100.0; // Initial position (m)
  104.     double vx0 = 3887.3; // Initial velocity (m/s)
  105.     double vy0 = 0.0;
  106.  
  107.     // Time parameters
  108.     double t_end = 86400.0; // Final time (one day)
  109.     double dt_output = 60.0; // Output data every 60 seconds
  110.  
  111.     // Solution vector [x, y, vx, vy]
  112.     std::valarray<double> y_init = {x0, y0, vx0, vy0};
  113.  
  114.     // Solve using Euler's method
  115.     solve_and_output(t0, t_end, dt_output, y_init, "orbit_data_euler.dat", 'E');
  116.  
  117.     // Solve using RK4 method
  118.     solve_and_output(t0, t_end, dt_output, y_init, "orbit_data_rk4_moon.dat", 'R');
  119.  
  120.     std::cout << "Data has been written to 'orbit_data_euler.dat' and 'orbit_data_rk4_moon.dat'." << std::endl;
  121.  
  122.     return 0;
  123. }
  124.  
Success #stdin #stdout 0.01s 5280KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Data has been written to 'orbit_data_euler.dat' and 'orbit_data_rk4_moon.dat'.
https://ideone.com/UN9hfX
language:
C++ (gcc 8.3)
created:
11 days ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!