Problem with Leap Frog Integration simulating the Earth around the Sun

I have been trying to implement Leap Frog integration for the system of the Earth orbiting the Sun:

#include <iostream>
#include <fstream>
#include <cmath>

using namespace std;


#define Grav    6.6726e-11             // m^3/kg/s^2
#define Msun    1.989e30               // kg
#define Mearth  5.973e24               // kg
#define Dse     1.496e11               // distance Sun-Earth in meters



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;

    // Delta time
    double dt = 0.01;

    //we create a file for saving the coordinates
    std::ofstream outputFile("LF.dat");
    outputFile << x0 << " " << y0 << std::endl;


    for (int i = 0; i<n; i++) {

        double x = x0 + dt*px0/(2*Mearth);
        double y = y0 + dt*py0/(2*Mearth);

        outputFile << x << " " << y << std::endl;

        double px = px0 - dt* Grav* Msun * x/pow(sqrt(x*x + y*y), 3);
        double py = py0 - dt* Grav* Msun* y/pow(sqrt(x*x + y*y), 3);

        double x2 = x + dt*px/(2*Mearth);
        double y2 = y + dt*py/(2*Mearth);

        
        outputFile << x2 << " " << y2 << std::endl;

        //initialize the values for the following step
        x0 = x;
        y0 = y;

        px0 = px;
        py0 = py;

    }

    outputFile.close();
}



int main() {

    LeapFrog(48);

    return 0;

}

...and I keep getting the same error where I only get a straight vertical line instead of an orbit.

I have tried almost everything and different implementations of the algorithm, the problem is also that I need to do a kick-drift-kick integration. I was expecting getting a correct orbit of the Earth around the Sun.

OK, you have a few problems - mostly physical.

Your code is going to simulate 48 * 0.01 = 0.48s, but there are 3.16e7 s in a year, so it won't be much of an orbit. Increase the number of steps and make the timestep a realistic fraction of the orbital time (circumference/velocity).

px and py refer to momentum, so you need total force, not total acceleration in their change: you forgot Mearth! Change to

double px = px0 - dt* Grav* Msun * Mearth * x/pow(sqrt(x*x + y*y), 3);

and similarly for py.

You should update x0 and y0 to x2 and y2 at the end of the step, not x and y. Also, you only need to output once a timestep.

The following plots a circular orbit. Change the initial velocity a bit if you want an ellipse.

#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 Dse = 1.496e11;           // distance Sun-Earth in meters
int nsteps = 1000;               // 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;

    // Timestep
    double dt = ( 2 * 3.141592654 * Dse / vy0 ) / nsteps;      // fraction of a year

    // Create a file for saving the coordinates
    ofstream outputFile("LF.dat");
    outputFile << x0 << " " << y0 << '\n';

    for (int i = 0; i<n; i++) {
        double x = x0 + dt*px0/(2*Mearth);
        double y = y0 + dt*py0/(2*Mearth);

        double px = px0 - dt * Grav* Msun * Mearth * x / pow(sqrt(x*x + y*y), 3);        // Don't forget the earth!
        double py = py0 - dt * Grav* Msun * Mearth * y / pow(sqrt(x*x + y*y), 3);

        double x2 = x + dt*px/(2*Mearth);
        double y2 = y + dt*py/(2*Mearth);

        outputFile << x2 << " " << y2 << std::endl;

        // Initialise the values for the following step
        x0 = x2;
        y0 = y2;

        px0 = px;
        py0 = py;
    }
}


int main()
{
    LeapFrog( nsteps * 1 );                   // One orbit
}