Why acos fails?

Question

OldCar 2015-10-25

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/250971-why-acos-fails

回答： Star Strider 2015-10-25

I have a Space Mechanics problem: I have a position vector and a velocity vector of a circular orbit with i=0. I have a problem with the function acos. In fact acos(something) is not 0 when something is 1. the code is

mi=398600.4;
x=7000;
y=0;
z=0;
u=0;
v=sqrt(mi/7000) ;
w=0;
r=sqrt(x^2+y^2+z^2);
V=sqrt((u^2+v^2+w^2));
a=mi*r/(2*mi-r*V^2);
e=sqrt((((V^2-mi/r)*x-(x*u+y*v+z*w)*u)/mi)^2+...
       (((V^2-mi/r)*y-(x*u+y*v+z*w)*v)/mi)^2+...
       (((V^2-mi/r)*z-(x*u+y*v+z*w)*w)/mi)^2);
if e>=1
    error('l''orbita non è una ellisse')
end
E=atan2(e*sqrt(1/mi*a)*(x*u+y*v+z*w),e*(1-r/a));
%tau=-(a^3/mi)^(1/3)*(E-e*sin(E));  %tempo dal passaggio al perigeo
ni=2*atan(sqrt((1+e)/(1-e))*tan(E/2));
H=sqrt(mi*a*(1-e^2));
i=acos((x*v-y*u)/H);
%if i<10^-5
 %   error('l''orbita è poco inclinata e bisogna corregere con un trick')
%end
raan=atan2((y*w-z*v)/(H*sin(i)),(x*w-z*u)/(H*sin(i)));
ot=atan2(z/(r*sin(i)),y/r*sin(raan)+x/r*cos(raan));    %omega+teta
omegap=ot-ni;