Matlab twobody algorithm is not working precisely

I checked matlab twobody algorithm with accoring to the AGI STK with Howard Curtis and David Vallado books and Matlab propagateOrbit function . AGI STK Curtis and Vallado book are compatible with each other. Matlab function inject an error 1m for every day totally 60 m for 60 day propagation. This Bug effects also HPOP algorithm

mu = 398600.4415;
load('STKtwoBody.mat')
r0 = [-6786.08143734318400674965  -1210.16092926545320551668 12.95464131748462754956];% in km
v0 = [   -0.15900031559336824660      0.97232513551400689966  7.54021643425097654045];% in km
startTime=datetime(2020,1,1,6,0,0); 
stopTime =datetime(2020,3,1,6,0,0); 
timeStep=60;
time = startTime:seconds(timeStep):stopTime;
Nsample=length(time);
rSTK=rSTK(:,1:Nsample);
vSTK=vSTK(:,1:Nsample);
km2m=1e3;
[rMatlab,vMatlab]= propagateOrbit(time,r0'*km2m,v0'*km2m,PropModel="two-body-keplerian");
rMatlab=rMatlab/km2m;
rVallado = zeros(3, Nsample); 
vVallado = zeros(3, Nsample);
rVallado(:,1)=r0;
vVallado(:,1)=v0;

rCurtis=rVallado;
vCurtis=vVallado;
for i=1:Nsample-1
    [rVallado(:,i+1),vVallado(:,i+1)] = kepler(r0,v0,timeStep*i,mu);
    [rCurtis(:,i+1),vCurtis(:,i+1)]   = rv_from_r0v0(r0,v0,timeStep*i,mu);
end
r=rSTK;
figure('color',[1 1 1]);
plot(time,vecnorm(r-rCurtis)*km2m,'b','linewidth',2); hold on;    
plot(time,vecnorm(r-rVallado)*km2m,'m','linewidth',2); hold on;
plot(time,vecnorm(r-rMatlab)*km2m,'k','linewidth',2); hold on;
xlabel('Time (day)'); ylabel('Position Magnitude Absolute Error (meters)');
legend('Curtis vs STK', 'Vallado vs STK', 'MATLAB vs STK', 'Location', 'best');
title('Comparison of Two-Body Propagators against STK');
grid on;
datetick('x','dd-mmm-yy','keeplimits');

●

Answer 1

[Torsten]

If you think you found a bug, contact MATLAB Support (Product Usage):

https://uk.mathworks.com/support/contact_us.html

Did you try the numerical method ?

numericalPropOpts = Aero.spacecraft.NumericalPropagatorOptions( ...
      ODESet=odeset(RelTol=1e-8,AbsTol=1e-8));
[rMatlab,vMatlab]= propagateOrbit(time,r0'*km2m,v0'*km2m,PropModel="numerical",NumericalPropagatorOptions=numericalPropOpts)

[Mehmet Mahmudoglu]

numericalPropOpts = Aero.spacecraft.NumericalPropagatorOptions( ...

ODESolver="ode89",...

ODESet=odeset(RelTol=1e-13,AbsTol=1e-15,MaxStep=1000), ...

IncludeThirdBodyGravity=false, ...

IncludeSRP=false,...

IncludeAtmosDrag=false,...

UseSpaceWeatherDataFile=false,...

IncludeAnomalousOxygen=false,...

GravitationalPotentialModel="point-mass");

[rMatlab,vMatlab]= propagateOrbit(time,r0'*km2m, v0'*km2m, ...

NumericalPropagatorOptions=numericalPropOpts);

Okey i will report Product usage i tried numeric method same result. I think twobody is core function executed both . Thanks