How to calculate the precise ignition duration during orbit changes

Hello everyone!
In my project, there is a sun-synchronous low-earth orbit satellite with a mass of approximately 50–150 kg. The initial orbit parameters of the satellite are known (Keplerian elements; the average altitude ranges from 350 km to 550 km). A Hall electric thruster is used for propulsion, providing a constant thrust of around 9–15 millinewtons. Fuel consumption can be neglected. The goal is to determine the required ignition duration to reach a target altitude (for example, 1500 meters above the current orbit). The maximum ignition time is 3 hours; if this limit is exceeded, multiple ignition sessions will be necessary.
I have developed a calculation algorithm, but the results differ by 10–200 seconds compared to the actual orbital changes.

public static double calculateDuration(
        double initialHeight,
        double targetHeight,
        double mass,
        double thrust) {

   
    double a1 = OrekitUtil.EARTH_RADIUS + initialHeight;
    double a2 = OrekitUtil.EARTH_RADIUS + targetHeight;

    // F = [m * √μ * (1/√a1 - 1/√a2)] / t
    double term1 = 1.0 / Math.sqrt(a1);
    double term2 = 1.0 / Math.sqrt(a2);
    double duration = mass * Math.sqrt(OrekitUtil.MU) * (term1 - term2) / thrust;

    return Math.abs(duration);
}

I’ve also considered extrapolation methods, but the semi-major axis varies significantly, making it impossible to accurately determine whether the target altitude is reached.
Is there any more accurate calculation method? Thank you all!

As you are using low thrust and long continuous burns, I doubt you can use simple Keplerian kinematics formulas. The propulsion formula (thrust, flow rate, duration…) does hold but it is the velocity increment that needs to be properly integrated. Is your formula already an analytical model of the integrated continuous thrust? What is the total duration of the maneuver? If you get 10 seconds difference on a 3 hours burn, the approximation is already good.

Hi @pig357, @luc,

Your formula looks good to me. @luc, yes, the formula is the integration of the first Gauss equation for a pure constant tangential thrust with a constant mass.
Are you using perturbations like the J2 ? I guess so since you say that the semi-major axis varies greatly. What may be happening is that you’re looking at the osculating elements while the formulas you’re using is, as Luc said, based on pure Keplerian (i.e. two-body) motion.
Can you try the same calculations with a pure Keplerian motion (i.e. no perturbations) and see if you get the same issue ?

Cheers,
Maxime