search

pbrod / nvector

Nvector is a suite of tools written in Python to solve geographical position calculations.
Other
57 stars 7 forks source link

Cross track distance ill conditioned #9

Closed oysstu closed 3 years ago

oysstu commented 3 years ago

Hi,

I was doing some quick experimentation with nvector, and found this calculation of the cross-track distance on the great circle. https://github.com/pbrod/nvector/blob/bf1cf5e1e210b74a57ea4bb2c277b388308bdba9/src/nvector/_core.py#L1140-L1143

The arccos function is inaccurate for small angles, as mentioned in the comment. When this function is called from cross_track_distance however, it is taken between the cross product of the n-vectors describing the path (c^E) and the position B. For small cross-track errors, shouldn't the expression with arccos be the most accurate (since arcsin is inaccurate at -pi/2 and pi/2)?

pbrod commented 3 years ago

No, the following example shows that np.arcsin delivers double precision both for small and large angles, but np.arcos solution gives less than double precision for most of the angles (even though the precision is ok for most applications) :

>>> import numpy as np
>>> theta = np.logspace(-50, 0)*np.pi/2
>>> sin_theta = np.sin(theta)
>>> t1 = np.arcsin(sin_theta)
>>> t2 = np.arcos(-sin_theta)-np.pi/2
>>> np.abs(t1-theta)/theta
array([0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,
   0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,
   0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,
   0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,
   0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,
   0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,
   0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,
   0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,
   0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,
   0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,
   0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,
   0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 1.8090627e-16,
   0.0000000e+00, 0.0000000e+00])

>>> np.abs(t2-theta)/theta
array([3.17768732e-07, 8.41285987e-08, 2.99824097e-08, 2.51113553e-07,
   3.65447084e-08, 6.71601921e-08, 1.13440618e-08, 6.73537798e-09,
   2.06761534e-08, 3.76403350e-09, 2.41976409e-09, 9.60252863e-10,
   1.34828226e-09, 3.70225909e-10, 4.83429429e-10, 1.75196678e-10,
   8.52621110e-11, 6.57949754e-11, 1.75711340e-10, 4.99285763e-12,
   8.42401296e-11, 4.23587289e-11, 3.14708133e-11, 1.77221187e-11,
   4.08716375e-12, 8.06297736e-12, 3.40705458e-12, 4.44254594e-13,
   1.95394667e-13, 1.00059067e-12, 7.88776657e-14, 4.48770241e-13,
   1.34252409e-13, 1.45198665e-13, 1.68480928e-14, 1.21205710e-14,
   4.80642950e-14, 2.70152645e-14, 1.37805543e-14, 1.01907463e-14,
   5.30814909e-15, 2.65430892e-15, 8.29544787e-16, 1.62960687e-15,
   1.85198952e-16, 3.47278700e-16, 4.34136709e-16, 1.80906270e-16,
   2.26152806e-16, 0.00000000e+00])
oysstu commented 3 years ago

Aha, I see. I assumed that the comment was related to footnote 2 in Gade (2010) regarding the numerical precision of the inverse trigonometric functions near the endpoints of their domains due to the steep slope.

I can confirm your findings on Linux, where presumably numpy uses the trigonometric functions of gcc/glibc. It looks like the asin and acos functions in glibc are implemented in terms of the atan machine instruction. The division by x in acos may be the cause for the numerical issues, but then again using long double does not help so perhaps not. asin = atan (x / sqrt(1 - x^2)) acos = atan (sqrt(1 - x^2) / x)

import numpy as np

theta = np.logspace(-50, 0, dtype=np.float64)*np.pi/2
cos_theta = np.cos(theta)
sin_theta = np.sin(theta)

t1 = np.arcsin(sin_theta)
t2 = np.arccos(-sin_theta) - np.pi/2
t3 = np.arctan2(sin_theta, cos_theta)

print('theta\n', theta)
print('t1\n', t1)
print('t2\n', t2)
print('t3\n', t3)

Output:

theta
 [1.57079633e-50 1.64637226e-49 1.72558439e-48 1.80860766e-47
 1.89562545e-46 1.98682993e-45 2.08242254e-44 2.18261442e-43
 2.28762684e-42 2.39769174e-41 2.51305220e-40 2.63396302e-39
 2.76069123e-38 2.89351674e-37 3.03273290e-36 3.17864719e-35
 3.33158186e-34 3.49187471e-33 3.65987975e-32 3.83596803e-31
 4.02052848e-30 4.21396870e-29 4.41671594e-28 4.62921799e-27
 4.85194418e-26 5.08538642e-25 5.33006030e-24 5.58650620e-23
 5.85529052e-22 6.13700690e-21 6.43227754e-20 6.74175457e-19
 7.06612151e-18 7.40609477e-17 7.76242520e-16 8.13589980e-15
 8.52734344e-14 8.93762066e-13 9.36763760e-12 9.81834400e-11
 1.02907353e-09 1.07858548e-08 1.13047961e-07 1.18487053e-06
 1.24187836e-05 1.30162902e-04 1.36425448e-03 1.42989303e-02
 1.49868967e-01 1.57079633e+00]
t1
 [1.57079633e-50 1.64637226e-49 1.72558439e-48 1.80860766e-47
 1.89562545e-46 1.98682993e-45 2.08242254e-44 2.18261442e-43
 2.28762684e-42 2.39769174e-41 2.51305220e-40 2.63396302e-39
 2.76069123e-38 2.89351674e-37 3.03273290e-36 3.17864719e-35
 3.33158186e-34 3.49187471e-33 3.65987975e-32 3.83596803e-31
 4.02052848e-30 4.21396870e-29 4.41671594e-28 4.62921799e-27
 4.85194418e-26 5.08538642e-25 5.33006030e-24 5.58650620e-23
 5.85529052e-22 6.13700690e-21 6.43227754e-20 6.74175457e-19
 7.06612151e-18 7.40609477e-17 7.76242520e-16 8.13589980e-15
 8.52734344e-14 8.93762066e-13 9.36763760e-12 9.81834400e-11
 1.02907353e-09 1.07858548e-08 1.13047961e-07 1.18487053e-06
 1.24187836e-05 1.30162902e-04 1.36425448e-03 1.42989303e-02
 1.49868967e-01 1.57079633e+00]
t2
 [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
 0.00000000e+00 2.22044605e-16 8.88178420e-16 8.21565038e-15
 8.52651283e-14 8.93729535e-13 9.36761779e-12 9.81834614e-11
 1.02907349e-09 1.07858549e-08 1.13047961e-07 1.18487053e-06
 1.24187836e-05 1.30162902e-04 1.36425448e-03 1.42989303e-02
 1.49868967e-01 1.57079633e+00]
t3
 [1.57079633e-50 1.64637226e-49 1.72558439e-48 1.80860766e-47
 1.89562545e-46 1.98682993e-45 2.08242254e-44 2.18261442e-43
 2.28762684e-42 2.39769174e-41 2.51305220e-40 2.63396302e-39
 2.76069123e-38 2.89351674e-37 3.03273290e-36 3.17864719e-35
 3.33158186e-34 3.49187471e-33 3.65987975e-32 3.83596803e-31
 4.02052848e-30 4.21396870e-29 4.41671594e-28 4.62921799e-27
 4.85194418e-26 5.08538642e-25 5.33006030e-24 5.58650620e-23
 5.85529052e-22 6.13700690e-21 6.43227754e-20 6.74175457e-19
 7.06612151e-18 7.40609477e-17 7.76242520e-16 8.13589980e-15
 8.52734344e-14 8.93762066e-13 9.36763760e-12 9.81834400e-11
 1.02907353e-09 1.07858548e-08 1.13047961e-07 1.18487053e-06
 1.24187836e-05 1.30162902e-04 1.36425448e-03 1.42989303e-02
 1.49868967e-01 1.57079633e+00]