Principia icon indicating copy to clipboard operation
Principia copied to clipboard

Numerical stability and angles

Open eggrobin opened this issue 10 years ago • 0 comments

The codebase has several angle computations whose numerical stability has not been investigated. In particular, https://github.com/mockingbirdnest/Principia/blob/69cf1a4/geometry/r3_element_body.hpp#L108 computes the latitude as arc sin(z / sqrt(x^2 + y^2 + z^2)), which is likely to cause trouble because of arc sin's singularity. arc tg(z / sqrt(x^2 + y^2)), computed with the 2-parameter arctangent, should be stabler, although a detailed analysis would be desirable.

In addition, https://github.com/mockingbirdnest/Principia/blob/69cf1a4/physics/ephemeris_test.cpp#L710 computes an angle between two vector as the 2-parameter arctangent of the norm of the cross product. While this does not present the singularity-induced instability of the arccosine method, this is not the method recommended by Kahan in https://www.cs.berkeley.edu/~wkahan/Mindless.pdf (p. 47). In any case, the angle of two vectors is a useful enough operation that it deserves a dedicated function in geometry: astronomical tests, e.g., @pdn4kd's eclipses, will need this.

eggrobin avatar Nov 24 '15 17:11 eggrobin