projectile physics – How to calculate trajectory on a planet with drag factor

I’m curious what is the valid calculation of a projectile trajectory in case when:

  • surface is not flat but spherical like a planet
  • drag is enabled (for example 0.1)
  • and we are in 3D

I’ve found calculations where either the drag or the spherical surface factor was missing. I assume in each simulation step I have to recalculate the direction of the gravity which is the direction between the current position of the projectile and center of the planet.

enter image description here

I’ve tried to reuse the SimulateArc method in one of these answers: https://stackoverflow.com/questions/61125224/2d-projectile-trajectory-predictionunity-2d

I changed Vector2 to Vector3 and added a gravity direction recalculation on the loop after the new position is calculated. But Drag factor is missing yet. Is it a good approach to solve this issue? If not how should I calculate the trajectory also with the Drag factor?

private void CalculateTrajectory()
{
    float simulateForDuration = 5f;//simulate for 5 secs in the furture
    float simulationStep = 0.1f;//Will add a point every 0.1 secs.

    int steps = (int)(simulateForDuration/simulationStep);//50 in this example
    Vector3 calculatedPosition;
    Vector3 directionVector = unit.CannonDirection;
    Vector3 launchPosition = unit.StartPosition;
    float launchSpeed = unit.ShootPower;

    for(int i = 0; i < steps; ++i)
    {
        calculatedPosition = launchPosition + ( directionVector * (launchSpeed * i * simulationStep));
        //Calculate gravity
        Vector3 gravity = (planet.transform.position - transform.position).normalized * gravityConst;
        calculatedPosition += gravity * ( i * simulationStep) *  ( i * simulationStep);
            
        Gizmos.DrawSphere(calculatedPosition, 0.5f);
    }

}

EDIT 1:
I tried to code the formula provided by Sacha but it seems it is not perfect yet (only on flat surface and without drag for now, if it works then I will add the drag factor and a sphere surface):

private void CalculateTrajectory()
{
    float simulateForDuration = 5f;
    float simulationStep = 0.1f;

    int steps = (int)(simulateForDuration / simulationStep);
    float mass = 0.5f;
    Vector3 previousPosition = transform.position;
    Vector3 previousVelocity = transform.forward.normalized * 5;
    Vector3 gravity = Physics.gravity;

    for (int i = 0; i < steps; ++i)
    {
        Vector3 acceleration = (previousVelocity + gravity) / mass;
        Vector3 newVelocity = previousVelocity + acceleration * simulationStep;
        Vector3 newPosition = previousPosition + newVelocity * simulationStep;

        Gizmos.DrawSphere(newPosition, 0.2f);

        previousPosition = newPosition;
        previousVelocity = newVelocity;
    }

}

enter image description here

Here is my projectile creation code:

GameObject clone = GameObject.Instantiate(sphere);
clone.transform.position = transform.position;
clone.transform.forward = transform.forward;

Rigidbody rb = clone.GetComponent<Rigidbody>();
rb.velocity = clone.transform.forward.normalized * 5f;

What did I miss?

EDIT 2:
I tried another approach:

float simulateForDuration = 5f;
float simulationStep = 0.02f;

int steps = (int)(simulateForDuration / simulationStep);
Vector3 initialVelocity = transform.forward.normalized * 5f;

for (int i = 0; i < steps; ++i)
{
    Vector3 pos = new Vector3(
        transform.position.x + initialVelocity.x * (i * simulationStep),
        transform.position.y + initialVelocity.y * (i * simulationStep) - 0.5f * -Physics.gravity.y * Mathf.Pow(i * simulationStep, 2),
        transform.position.z + initialVelocity.z * (i * simulationStep));
        
    Gizmos.DrawSphere(pos, 0.05f);
}

Result:
enter image description here

This seems fine when I have only gravity on axis Y. We can say that the first version is almost fine, but I cannot figure out what is missing.

EDIT Final:
Based on Sacha updated answer regarding P = m * g the calculation is fine. I also had to update the AddForce method on the projectile itself because I added 9.81 * normalizedDirectionToTheCenter for the projectile as gravity. Now it is 9.81 * projectileMass * normalizedDirectionToTheCenter .

Here is the full code:

private void CalculateTrajectory()
{
    float simulateForDuration = 10f;
    float simulationStep = 0.02f;

    int steps = (int)(simulateForDuration / simulationStep);
    float mass = r.mass; // Mass of the projectile.

    Vector3 previousPosition = transform.position;
    Vector3 previousVelocity = transform.forward.normalized * power;

    for (int i = 1; i < steps; ++i)
    {
        Vector3 directionToCenter = (planet.position - previousPosition).normalized;
        Vector3 gravity = directionToCenter * 9.81f;

        // Separated vectors to understand better.
        Vector3 force = gravity * mass;
        Vector3 acceleration = force / mass;
        Vector3 newVelocity = previousVelocity + acceleration * simulationStep;
        Vector3 newPosition = previousPosition + newVelocity * simulationStep;

        Gizmos.DrawSphere(newPosition, 0.05f);

        previousPosition = newPosition;
        previousVelocity = newVelocity;
    }
}

Thanks Sacha!

enter image description here