Matlab twobody algorithm is not working precisely

8 次查看（过去 30 天）
Mehmet Mahmudoglu
Mehmet Mahmudoglu 2025-7-25
评论： Mehmet Mahmudoglu 2025-7-25
Ran in:
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');
●

请先登录，再回答此问题。

回答（1 个）

Torsten
Torsten 2025-7-25
移动：Torsten 2025-7-25
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)
  1 个评论
Mehmet Mahmudoglu
Mehmet Mahmudoglu 2025-7-25
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

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Aerospace Applications 的更多信息

标签

产品


版本

R2025a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by