I am trying to learn C++ but I am really struggling with it. I’d like to translate my Python code to C++. Would anyone help me?

import numpy as np
import matplotlib.pyplot as plt

# Constants
Grav = 6.6726e-11  # m^3/kg/s^2 (Universal gravitational constant)
Msun = 1.989e30    # kg
Mearth = 5.973e24  # kg
Dse = 1.496e11     # distance Sun-Earth in meters

# Initial conditions
x0 = Dse  # Initial x position of Earth
y0 = 0    # Initial y position of Earth
vx0 = 0   # Initial velocity along x-axis
vy0 = np.sqrt(Grav*Msun/Dse)  # Initial velocity along y-axis (circular orbit velocity)

# Integration parameters
steps = 48
dt = 5.1 * 365.25 * 24 * 3600 / steps  # Time step for 5 years in seconds

# Function to calculate derivatives
def derivs(x, y, vx, vy):
    r = np.sqrt(x**2 + y**2)
    dvx_dt = -Grav * Msun * x / r**3
    dvy_dt = -Grav * Msun * y / r**3
    return vx, vy, dvx_dt, dvy_dt

# 4th order Runge-Kutta integration scheme
def runge_kutta_integration(x0, y0, vx0, vy0, dt, steps):
    x = x0
    y = y0
    vx = vx0
    vy = vy0
    with open("rk4.dat", "w") as file:
        for _ in range(steps):
            k1x, k1y, k1vx, k1vy = derivs(x, y, vx, vy)
            k2x, k2y, k2vx, k2vy = derivs(x + 0.5 * k1x * dt, y + 0.5 * k1y * dt, vx + 0.5 * k1vx * dt, vy + 0.5 * k1vy * dt)
            k3x, k3y, k3vx, k3vy = derivs(x + 0.5 * k2x * dt, y + 0.5 * k2y * dt, vx + 0.5 * k2vx * dt, vy + 0.5 * k2vy * dt)
            k4x, k4y, k4vx, k4vy = derivs(x + k3x * dt, y + k3y * dt, vx + k3vx * dt, vy + k3vy * dt)
            x += (k1x + 2 * k2x + 2 * k3x + k4x) * dt / 6
            y += (k1y + 2 * k2y + 2 * k3y + k4y) * dt / 6
            vx += (k1vx + 2 * k2vx + 2 * k3vx + k4vx) * dt / 6
            vy += (k1vy + 2 * k2vy + 2 * k3vy + k4vy) * dt / 6
            file.write(f"{x} {y}n")

# Plot trajectory
def plot_trajectory():
    data = np.loadtxt("rk4.dat")
    plt.plot(data[:, 0], data[:, 1], label="4th order Runge-Kutta")
    plt.xlabel("x [m]")
    plt.ylabel("y [m]")
    plt.title("Orbit of Earth about the Sun")
    plt.legend()
    plt.grid(True)
    plt.gca().set_aspect('equal', adjustable='box')
    plt.show()

# Perform integration
runge_kutta_integration(x0, y0, vx0, vy0, dt, steps)

# Plot trajectory
plot_trajectory()

Final result of my code

I tried to go step by step in the code but I only get errors. I didn’t try to write some lines for plotting because I know how to plot in gnuplot, so that’s not a problem.

#include <iostream>
#include <fstream>
#include <cmath>
using namespace std;

// Some constants
double Grav = 6.6726e-11;        // m^3/kg/s^2 (Universal gravitational constant)
double Msun = 1.989e30;          // kg
double Mearth = 5.973e24;        // kg
double Dse = 1.496e11;           // distance Sun-Earth in meters

// Initial conditions
double x = Dse  // Initial x position of Earth
double y = 0    // Initial y position of Earth
double vx = 0   // Initial velocity along x-axis
double vy = sqrt(Grav*Msun/Dse)  // Initial velocity along y-axis (circular orbit velocity)

// Integration parameters
steps = 48
dt = 5.1 * 365.25 * 24 * 3600 / steps  // Time step for 5 years in seconds

// Function to calculate derivatives
void derivs(double x, double y, double vx, double vy) {
    double r = sqrt(x*x + y*y);
    double dvx_dt = -Grav * Msun * x / pow(r, 3);
    double dvy_dt = -Grav * Msun * y / pow(r, 3);
    return vx, vy, dvx_dt, dvy_dt;
}

// 4th order Runge-Kutta integration scheme
void runge_kutta_integration(double x, double y, double vx, double vy, double dt, int steps) {
    ofstream file("rk4.dat");
    for (int i = 0; i < steps; i++) {
        double k1x, k1y, k1vx, k1vy;
        double k2x, k2y, k2vx, k2vy;
        double k3x, k3y, k3vx, k3vy;
        double k4x, k4y, k4vx, k4vy;

        derivs(x, y, vx, vy, &k1x, &k1y, &k1vx, &k1vy);
        derivs(x + 0.5 * k1x * dt, y + 0.5 * k1y * dt, vx + 0.5 * k1vx * dt, vy + 0.5 * k1vy * dt, &k2x, &k2y, &k2vx, &k2vy);
        derivs(x + 0.5 * k2x * dt, y + 0.5 * k2y * dt, vx + 0.5 * k2vx * dt, vy + 0.5 * k2vy * dt, &k3x, &k3y, &k3vx, &k3vy);
        derivs(x + k3x * dt, y + k3y * dt, vx + k3vx * dt, vy + k3vy * dt, &k4x, &k4y, &k4vx, &k4vy);

        x += (k1x + 2 * k2x + 2 * k3x + k4x) * dt / 6;
        y += (k1y + 2 * k2y + 2 * k3y + k4y) * dt / 6;
        vx += (k1vx + 2 * k2vx + 2 * k3vx + k4vx) * dt / 6;
        vy += (k1vy + 2 * k2vy + 2 * k3vy + k4vy) * dt / 6;

        file << x << " " << y << "n";
    }
    file.close();
}

// Perform integration
int main() {
    runge_kutta_integration(x0, y0, vx0, vy0, dt, steps);
}