I am trying to calculate radius of created lens by two overlapped spheres. In this regard, I tried both trigonometric method and another method based on just algebraic. I compared the results by these two methods with various data sets and find a small number of contradictions on just some of those data sets. The results are the same in most cases. The problem can be reproduced by the following example (on 3-5 indices):
poss = np.array([[[-0.884, -3.45, -0.99 ], [-0.901, -3.43, -0.995]], [[-0.993, -3.44, -0.97 ], [-1.01, -3.46, -1. ]],
[[-0.993, -3.44, -0.97 ], [-0.998, -3.45, -1. ]], [[0.885 , 0.967, -1.02 ], [0.885, 0.964, -1.02] ],
[[-0.252, -3.3 , -0.777], [-0.197, -3.3 , -0.777]], [[0.26 , -1.68, -0.803], [0.288, -1.67, -0.799]],
[[0.599 , 2.04 , -0.857], [0.607 , 2.04 , -0.84 ]], [[0.615 , 2. , -0.833], [0.633, 2. , -0.855]],
[[0.698 , 2.06 , -0.921], [0.679 , 2.06 , -0.914]]])
rad = np.array([[0.0108, 0.0205], [0.0231, 0.0259], [0.0231 , 0.0304], [0.0154, 0.0124], [0.0137, 0.0413],
[0.027 , 0.003 ], [0.0102, 0.022 ], [0.00221, 0.0268], [0.0147, 0.0124]])
# The length of the overlaps; lenses' heights
gap = np.array([-4.57922157e-03, -9.13773714e-03, -2.14843788e-02, -2.48000000e-02, -1.38777878e-17, -2.42861287e-17,
-1.34117058e-02, -5.84659193e-04, -6.85154327e-03])
The functions are:
def trigonometric(r_active, gap):
r_add = np.add.reduce(r_active, axis=1)
paired_cent_dis = np.sum((r_add, gap), axis=0)
intersect_angle_0 = np.arccos(np.clip((r_active[:, 0] ** 2 +
paired_cent_dis ** 2 - r_active[:, 1] ** 2) /
(2 * r_active[:, 0] * paired_cent_dis), -1, 1))
intersect_plane_rad = r_active[:, 0] * np.sin(intersect_angle_0)
return intersect_plane_rad
def algebraic(r, gap):
items_ = np.empty((len(gap), 1), dtype=np.float64)
for i in range(len(gap)):
r0, r1 = r[i]
cur_gap = gap[i]
paired_cent_dis = r0 + r1 + cur_gap
intersect_plane_rad = 0.5 * abs((-paired_cent_dis + r0 + r1) *
( paired_cent_dis + r0 + r1) * (-paired_cent_dis - r0 + r1) *
(-paired_cent_dis + r0 - r1)) ** 0.5 / paired_cent_dis
items_[i] = intersect_plane_rad
return items_.ravel()
trigonometric(rad, gap)
algebraic(rad, gap)
The results:
# repr trigonometric:
array([7.59403901e-03, 1.42126146e-02, 2.08670250e-02, 0.00000000e+00,
4.56484128e-10, 0.00000000e+00, 1.01747354e-02, 1.45347671e-03,
8.94740633e-03])
# repr algebraic:
array([7.59403901e-03, 1.42126146e-02, 2.08670250e-02, 4.69938148e-10,
5.34354024e-10, 3.68549655e-10, 1.01747354e-02, 1.45347671e-03,
8.94740633e-03])
As it can be seen by the results, there are some different resulted values on indices 3, 4, and 5. AFAIK, the two methods do the same job; It is proved by various data volumes. But some such differences may be happened on some indices in rare cases. In this example, just the 3rd index is affected by np.clip (this index in this small example gets 0 by trigonometric method, but it gets a nonzero value in my main code!? That nonzero value, too, was different from the same index resulted value by algebraic method i.e. 4.69938148e-10). As it is obvious in the images, and by focusing on the gap values (that are very small or near the diameter size of the smaller sphere), it seems the problem (differences between the results on some contacts) will be due to calculation precisions or something like that.
The final algebraic result shows the number of decimals for suspected indices will be in a reasonable range (here it is 10) and it seems trigonometric method is misled during the process.
I would be grateful to find
- where is the problem source,
- why the 4th index of trigonometric result gets a nonzero value but by a different magnitude although the 4th and 5th gap values are near the same,
- and how could cure trigonometric method if it could.