Coverage for pygeodesy/resections.py: 97%
365 statements
« prev ^ index » next coverage.py v7.6.1, created at 2024-11-12 16:17 -0500
« prev ^ index » next coverage.py v7.6.1, created at 2024-11-12 16:17 -0500
2# -*- coding: utf-8 -*-
4u'''3-Point resection functions L{cassini}, L{collins5}, L{pierlot}, L{pierlotx} and
5L{tienstra7}, survey functions L{snellius3} and L{wildberger3} and triangle functions
6L{triAngle}, L{triAngle5}, L{triSide}, L{triSide2} and L{triSide4}.
8@note: Functions L{pierlot} and L{pierlotx} are transcoded to Python with permission from
9 U{triangulationPierlot<http://www.Telecom.ULg.ac.BE/triangulation/doc/total_8c.html>} and
10 U{Pierlot<http://www.Telecom.ULg.ac.BE/publi/publications/pierlot/Pierlot2014ANewThree>}.
11'''
12# make sure int/int division yields float quotient
13from __future__ import division as _; del _ # PYCHOK semicolon
15from pygeodesy.basics import map1, map2, _zip, _ALL_LAZY
16from pygeodesy.constants import EPS, EPS0, EPS02, INT0, PI, PI2, PI_2, PI_4, \
17 _0_0, _0_5, _1_0, _N_1_0, _2_0, _N_2_0, _4_0, \
18 _16_0, _180_0, _360_0, _copysign_0_0, isnear0, \
19 _over, _umod_360
20from pygeodesy.errors import _and, _or, TriangleError, _ValueError, _xcallable, \
21 _xkwds, _xkwds_pop2
22from pygeodesy.fmath import favg, Fdot, fidw, fmean, hypot, hypot2_
23from pygeodesy.fsums import _Fsumf_, fsumf_, fsum1, fsum1f_
24from pygeodesy.interns import _a_, _A_, _area_, _b_, _B_, _c_, _C_, _coincident_, \
25 _colinear_, _d_, _eps_, _invalid_, _negative_, _not_, \
26 _rIn_, _SPACE_
27# from pygeodesy.lazily import _ALL_LAZY # from .basics
28from pygeodesy.named import _NamedTuple, _Pass, Fmt
29# from pygeodesy.streprs import Fmt # from .named
30from pygeodesy.units import Degrees, Distance, Radians
31from pygeodesy.utily import acos1, asin1, sincos2, sincos2_, sincos2d, sincos2d_
32from pygeodesy.vector3d import _otherV3d, Vector3d
34from math import cos, atan2, degrees, fabs, radians, sin, sqrt
36__all__ = _ALL_LAZY.resections
37__version__ = '24.11.04'
39_concyclic_ = 'concyclic'
40_PA_ = 'PA'
41_PB_ = 'PB'
42_PC_ = 'PC'
43_pointH_ = 'pointH'
44_pointP_ = 'pointP'
45_positive_ = 'positive'
46_radA_ = 'radA'
47_radB_ = 'radB'
48_radC_ = 'radC'
51class Collins5Tuple(_NamedTuple):
52 '''5-Tuple C{(pointP, pointH, a, b, c)} with survey C{pointP}, auxiliary
53 C{pointH}, each an instance of B{C{pointA}}'s (sub-)class and triangle
54 sides C{a}, C{b} and C{c} in C{meter}, conventionally.
55 '''
56 _Names_ = (_pointP_, _pointH_, _a_, _b_, _c_)
57 _Units_ = (_Pass, _Pass, Distance, Distance, Distance)
60class ResectionError(_ValueError):
61 '''Error raised for issues in L{pygeodesy.resections}.
62 '''
63 pass
66class Survey3Tuple(_NamedTuple):
67 '''3-Tuple C{(PA, PB, PC)} with distance from survey point C{P} to each of
68 the triangle corners C{A}, C{B} and C{C} in C{meter}, conventionally.
69 '''
70 _Names_ = (_PA_, _PB_, _PC_)
71 _Units_ = ( Distance, Distance, Distance)
74class Tienstra7Tuple(_NamedTuple):
75 '''7-Tuple C{(pointP, A, B, C, a, b, c)} with survey C{pointP}, interior
76 triangle angles C{A}, C{B} and C{C} in C{degrees} and triangle sides
77 C{a}, C{b} and C{c} in C{meter}, conventionally.
78 '''
79 _Names_ = (_pointP_, _A_, _B_, _C_, _a_, _b_, _c_)
80 _Units_ = (_Pass, Degrees, Degrees, Degrees, Distance, Distance, Distance)
83class TriAngle5Tuple(_NamedTuple):
84 '''5-Tuple C{(radA, radB, radC, rIn, area)} with the interior angles at
85 triangle corners C{A}, C{B} and C{C} in C{radians}, the C{InCircle}
86 radius C{rIn} aka C{inradius} in C{meter} and the triangle C{area}
87 in C{meter} I{squared}, conventionally.
88 '''
89 _Names_ = (_radA_, _radB_, _radC_, _rIn_, _area_)
90 _Units_ = ( Radians, Radians, Radians, Distance, _Pass)
93class TriSide2Tuple(_NamedTuple):
94 '''2-Tuple C{(a, radA)} with triangle side C{a} in C{meter}, conventionally
95 and angle C{radA} at the opposite triangle corner in C{radians}.
96 '''
97 _Names_ = (_a_, _radA_)
98 _Units_ = ( Distance, Radians)
101class TriSide4Tuple(_NamedTuple):
102 '''4-Tuple C{(a, b, radC, d)} with interior angle C{radC} at triangle corner
103 C{C} in C{radians}with the length of triangle sides C{a} and C{b} and
104 with triangle height C{d} perpendicular to triangle side C{c}, in the
105 same units as triangle sides C{a} and C{b}.
106 '''
107 _Names_ = (_a_, _b_, _radC_, _d_)
108 _Units_ = ( Distance, Distance, Radians, Distance)
111def _ABC3(useZ, pointA, pointB, pointC):
112 '''(INTERNAL) Helper for L{cassini} and L{tienstra7}.
113 '''
114 return (_otherV3d(useZ=useZ, pointA=pointA),
115 _otherV3d(useZ=useZ, pointB=pointB),
116 _otherV3d(useZ=useZ, pointC=pointC))
119def _B3(useZ, point1, point2, point3):
120 '''(INTERNAL) Helper for L{pierlot} and L{pierlotx}.
121 '''
122 return (_otherV3d(useZ=useZ, point1=point1),
123 _otherV3d(useZ=useZ, point2=point2),
124 _otherV3d(useZ=useZ, point3=point3))
127def cassini(pointA, pointB, pointC, alpha, beta, useZ=False, **Clas_and_kwds):
128 '''3-Point resection using U{Cassini<https://NL.WikiPedia.org/wiki/Achterwaartse_insnijding>}'s method.
130 @arg pointA: First point (C{Cartesian}, L{Vector3d}, C{Vector3Tuple},
131 C{Vector4Tuple} or C{Vector2Tuple} if C{B{useZ}=False}).
132 @arg pointB: Second point (C{Cartesian}, L{Vector3d}, C{Vector3Tuple},
133 C{Vector4Tuple} or C{Vector2Tuple} if C{B{useZ}=False}).
134 @arg pointC: Center point (C{Cartesian}, L{Vector3d}, C{Vector3Tuple},
135 C{Vector4Tuple} or C{Vector2Tuple} if C{B{useZ}=False}).
136 @arg alpha: Angle subtended by triangle side B{C{pointA}} to B{C{pointC}}
137 (C{degrees}, non-negative).
138 @arg beta: Angle subtended by triangle side B{C{pointB}} to B{C{pointC}}
139 (C{degrees}, non-negative).
140 @kwarg useZ: If C{True}, use and interpolate the Z component, otherwise
141 force C{z=INT0} (C{bool}).
142 @kwarg Clas_and_kwds: Optional class C{B{Clas}=B{pointA}.classof} to
143 return the survey point with optionally other B{C{Clas}}
144 keyword arguments to instantiate the survey point.
146 @note: Typically, B{C{pointC}} is between B{C{pointA}} and B{C{pointB}}.
148 @return: The survey point, an instance of B{C{Clas}} or B{C{pointA}}'s
149 (sub-)class.
151 @raise ResectionError: Near-coincident, -colinear or -concyclic points
152 or negative or invalid B{C{alpha}} or B{C{beta}}.
154 @raise TypeError: Invalid B{C{pointA}}, B{C{pointB}} or B{C{pointM}}.
156 @see: U{Three Point Resection Problem<https://Dokumen.tips/documents/
157 three-point-resection-problem-introduction-kaestner-burkhardt-method.html>}
158 and functions L{collins5}, L{pierlot}, L{pierlotx} and L{tienstra7}.
159 '''
161 def _H(A, C, sa):
162 s, c = sincos2d(sa)
163 if isnear0(s):
164 raise ValueError(_or(_coincident_, _colinear_))
165 t = s, c, c
166 x = Fdot(t, A.x, C.y, -A.y).fover(s)
167 y = Fdot(t, A.y, -C.x, A.x).fover(s)
168 return x, y
170 A, B, C = _ABC3(useZ, pointA, pointB, pointC)
171 try:
172 sa, sb = map1(float, alpha, beta)
173 if min(sa, sb) < 0:
174 raise ValueError(_negative_)
175 if fsumf_(_360_0, -sa, -sb) < EPS0:
176 raise ValueError()
178 x1, y1 = _H(A, C, sa)
179 x2, y2 = _H(B, C, -sb)
181 x = x1 - x2
182 y = y1 - y2
183 if isnear0(x) or isnear0(y):
184 raise ValueError(_SPACE_(_concyclic_, (x, y)))
186 m = y / x
187 n = x / y
188 N = n + m
189 if isnear0(N):
190 raise ValueError(_SPACE_(_concyclic_, (m, n, N)))
192 t = n, m, _1_0, _N_1_0
193 x = Fdot(t, C.x, x1, C.y, y1).fover(N)
194 y = Fdot(t, y1, C.y, C.x, x1).fover(N)
195 z = _zidw(x, y, useZ, A, B, C)
196 return _Clas(cassini, pointA, Clas_and_kwds, x, y, z)
198 except (TypeError, ValueError) as x:
199 raise ResectionError(pointA=pointA, pointB=pointB, pointC=pointC,
200 alpha=alpha, beta=beta, cause=x)
203def _Clas(which, point, Clas_and_kwds, *args):
204 '''(INTERNAL) Return a C{B{Clas}=point.classof} survey point.
205 '''
206 Clas, kwds = _xkwds_pop2(Clas_and_kwds, Clas=point.classof)
207 return Clas(*args, **_xkwds(kwds, name=which.__name__))
210def collins5(pointA, pointB, pointC, alpha, beta, useZ=False, **Clas_and_kwds):
211 '''3-Point resection using U{Collins<https://Dokumen.tips/documents/
212 three-point-resection-problem-introduction-kaestner-burkhardt-method.html>}' method.
214 @arg pointA: First point (C{Cartesian}, L{Vector3d}, C{Vector3Tuple},
215 C{Vector4Tuple} or C{Vector2Tuple} if C{B{useZ}=False}).
216 @arg pointB: Second point (C{Cartesian}, L{Vector3d}, C{Vector3Tuple},
217 C{Vector4Tuple} or C{Vector2Tuple} if C{B{useZ}=False}).
218 @arg pointC: Center point (C{Cartesian}, L{Vector3d}, C{Vector3Tuple},
219 C{Vector4Tuple} or C{Vector2Tuple} if C{B{useZ}=False}).
220 @arg alpha: Angle subtended by triangle side C{b} from B{C{pointA}} to
221 B{C{pointC}} (C{degrees}, non-negative).
222 @arg beta: Angle subtended by triangle side C{a} from B{C{pointB}} to
223 B{C{pointC}} (C{degrees}, non-negative).
224 @kwarg useZ: If C{True}, use and interpolate the Z component, otherwise
225 force C{z=INT0} (C{bool}).
226 @kwarg Clas_and_kwds: Optional class C{B{Clas}=B{pointA}.classof} to
227 return the survey point with optionally other B{C{Clas}}
228 keyword arguments to instantiate the survey point.
230 @note: Typically, B{C{pointC}} is between B{C{pointA}} and B{C{pointB}}.
232 @return: L{Collins5Tuple}C{(pointP, pointH, a, b, c)} with survey C{pointP},
233 auxiliary C{pointH}, each an instance of B{C{Clas}} or B{C{pointA}}'s
234 (sub-)class and triangle sides C{a}, C{b} and C{c} in C{meter},
235 conventionally.
237 @raise ResectionError: Near-coincident, -colinear or -concyclic points
238 or negative or invalid B{C{alpha}} or B{C{beta}}.
240 @raise TypeError: Invalid B{C{pointA}}, B{C{pointB}} or B{C{pointM}}.
242 @see: U{Collins' methode<https://NL.WikiPedia.org/wiki/Achterwaartse_insnijding>}
243 and functions L{cassini}, L{pierlot}, L{pierlotx} and L{tienstra7}.
244 '''
246 def _azi_len2(A, B, pi2=PI2):
247 v = B.minus(A)
248 r = atan2(v.x, v.y)
249 if r < 0 and pi2:
250 r += pi2
251 return r, v.length
253 def _xyz(d, r, A, B, C, useZ):
254 s, c = sincos2(r)
255 x = d * s + A.x # fma(d, s, A.x)
256 y = d * c + A.y # fma(d, c, A.y)
257 z = _zidw(x, y, useZ, A, B, C)
258 return x, y, z
260 A, B, C = _ABC3(useZ, pointA, pointB, pointC)
261 try:
262 ra, rb = radians(alpha), radians(beta)
263 if min(ra, rb) < 0:
264 raise ValueError(_negative_)
266 sra, srH = sin(ra), sin(ra + rb - PI) # rH = PI - ((PI - ra) + (PI - rb))
267 if isnear0(sra) or isnear0(srH):
268 raise ValueError(_or(_coincident_, _colinear_, _concyclic_))
270# za, a = _azi_len2(C, B)
271 zb, b = _azi_len2(C, A)
272 zc, c = _azi_len2(A, B, 0)
274# d = c * sin(PI - rb) / srH # B.minus(H).length
275 d = c * sin(PI - ra) / srH # A.minus(H).length
276 r = zc + PI - rb # zh = zc + (PI - rb)
277 H = _xyz(d, r, A, B, C, useZ)
279 zh, _ = _azi_len2(C, Vector3d(*H))
281# d = a * sin(za - zh) / sin(rb) # B.minus(P).length
282 d = b * sin(zb - zh) / sra # A.minus(P).length
283 r = zh - ra # zb - PI + (PI - ra - (zb - zh))
284 P = _xyz(d, r, A, B, C, useZ)
286 P = _Clas(collins5, pointA, Clas_and_kwds, *P)
287 H = _Clas(collins5, pointA, Clas_and_kwds, *H)
288 a = B.minus(C).length
289 return Collins5Tuple(P, H, a, b, c, name=collins5.__name__)
291 except (TypeError, ValueError) as x:
292 raise ResectionError(pointA=pointA, pointB=pointB, pointC=pointC,
293 alpha=alpha, beta=beta, cause=x)
296def pierlot(point1, point2, point3, alpha12, alpha23, useZ=False, eps=EPS,
297 **Clas_and_kwds):
298 '''3-Point resection using U{Pierlot<http://www.Telecom.ULg.ac.BE/publi/publications/
299 pierlot/Pierlot2014ANewThree>}'s method C{ToTal} with I{approximate} limits for
300 the (pseudo-)singularities.
302 @arg point1: First point (C{Cartesian}, L{Vector3d}, C{Vector3Tuple},
303 C{Vector4Tuple} or C{Vector2Tuple} if C{B{useZ}=False}).
304 @arg point2: Second point (C{Cartesian}, L{Vector3d}, C{Vector3Tuple},
305 C{Vector4Tuple} or C{Vector2Tuple} if C{B{useZ}=False}).
306 @arg point3: Third point (C{Cartesian}, L{Vector3d}, C{Vector3Tuple},
307 C{Vector4Tuple} or C{Vector2Tuple} if C{B{useZ}=False}).
308 @arg alpha12: Angle subtended from B{C{point1}} to B{C{point2}} or
309 B{C{alpha2 - alpha1}} (C{degrees}).
310 @arg alpha23: Angle subtended from B{C{point2}} to B{C{point3}} or
311 B{C{alpha3 - alpha2}}(C{degrees}).
312 @kwarg useZ: If C{True}, interpolate the survey point's Z component,
313 otherwise use C{z=INT0} (C{bool}).
314 @kwarg eps: Tolerance for C{cot}angent (pseudo-)singularities (C{float}).
315 @kwarg Clas_and_kwds: Optional class C{B{Clas}=B{point1}.classof} to
316 return the survey point with optionally other B{C{Clas}}
317 keyword arguments to instantiate the survey point.
319 @note: Typically, B{C{point1}}, B{C{point2}} and B{C{point3}} are ordered
320 by angle, modulo 360, counter-clockwise.
322 @return: The survey (or robot) point, an instance of B{C{Clas}} or B{C{point1}}'s
323 (sub-)class.
325 @raise ResectionError: Near-coincident, -colinear or -concyclic points
326 or invalid B{C{alpha12}} or B{C{alpha23}} or
327 non-positive B{C{eps}}.
329 @raise TypeError: Invalid B{C{point1}}, B{C{point2}} or B{C{point3}}.
331 @see: I{Pierlot}'s C function U{triangulationPierlot<http://www.Telecom.ULg.ac.BE/
332 triangulation/doc/total_8c_source.html>}, U{V. Pierlot, M. Van Droogenbroeck,
333 "A New Three Object Triangulation Algorithm for Mobile Robot Positioning"
334 <https://ORBi.ULiege.BE/bitstream/2268/157469/1/Pierlot2014ANewThree.pdf>},
335 U{Vincent Pierlot, Marc Van Droogenbroeck, "18 Triangulation Algorithms for 2D
336 Positioning (also known as the Resection Problem)"<http://www.Telecom.ULg.ac.BE/
337 triangulation>} and functions L{pierlotx}, L{cassini}, L{collins5} and L{tienstra7}.
338 '''
340 def _cot(s, c): # -eps < I{approximate} cotangent < eps
341 if eps > 0:
342 return c / (min(s, -eps) if s < 0 else max(s, eps))
343 raise ValueError(_SPACE_(_eps_, _not_, _positive_))
345 B1, B2, B3 = _B3(useZ, point1, point2, point3)
346 try:
347 xyz = _pierlot3(B1, B2, B3, alpha12, alpha23, useZ, _cot)
348 return _Clas(pierlot, point1, Clas_and_kwds, *xyz)
350 except (TypeError, ValueError) as x:
351 raise ResectionError(point1=point1, point2=point2, point3=point3,
352 alpha12=alpha12, alpha23=alpha23, eps=eps, cause=x)
355def _pierlot3(B1, B2, B3, a12, a23, useZ, cot):
356 '''(INTERNAL) Shared L{pierlot} and L{pierlotx}.
357 '''
358 x1_, y1_, _ = B1.minus(B2).xyz3
359 x3_, y3_, _ = B3.minus(B2).xyz3
361 s12, c12, s23, c23 = sincos2d_(a12, a23)
362 # cot31 = (1 - cot12 * cot23) / (cot12 + cot32)
363 # = (1 - c12 / s12 * c23 / s23) / (c12 / s12 + c23 / s23)
364 # = (1 - (c12 * c23) / (s12 * s23)) / (c12 * s23 + s12 * c23) / (s12 * s23)
365 # = (s12 * s23 - c12 * c23) / (c12 * s23 + s12 * c23)
366 # = c31 / s31
367 cot31 = cot(fsum1f_(c12 * s23, s12 * c23), # s31
368 fsum1f_(s12 * s23, -c12 * c23)) # c31
370 K = _Fsumf_(x3_ * x1_, cot31 * (y3_ * x1_),
371 y3_ * y1_, -cot31 * (x3_ * y1_))
372 if K:
373 cot12 = cot(s12, c12)
374 cot23 = cot(s23, c23)
376 # x12 = x1_ + cot12 * y1_
377 # y12 = y1_ - cot12 * x1_
379 # x23 = x3_ - cot23 * y3_
380 # y23 = y3_ + cot23 * x3_
382 # x31 = x3_ + x1_ + cot31 * (y3_ - y1_)
383 # y31 = y3_ + y1_ - cot31 * (x3_ - x1_)
385 # x12 - x23 = x1_ + cot12 * y1_ - x3_ + cot23 * y3_
386 X12_23 = _Fsumf_(x1_, cot12 * y1_, -x3_, cot23 * y3_)
387 # y12 - y23 = y1_ - cot12 * x1_ - y3_ - cot23 * x3_
388 Y12_23 = _Fsumf_(y1_, -cot12 * x1_, -y3_, -cot23 * x3_)
390 # x31 - x23 = x3_ + x1_ + cot31 * (y3_ - y1_) - x3_ + cot23 * y3_
391 # = x1_ + cot31 * y3_ - cot31 * y1_ + cot23 * y3_
392 X31_23 = _Fsumf_(x1_, -cot31 * y1_, cot31 * y3_, cot23 * y3_)
393 # y31 - y23 = y3_ + y1_ - cot31 * (x3_ - x1_) - y3_ - cot23 * x3_
394 # = y1_ - cot31 * x3_ + cot31 * x1_ - cot23 * x3_
395 Y31_23 = _Fsumf_(y1_, cot31 * x1_, -cot31 * x3_, -cot23 * x3_)
397 # d = (x12 - x23) * (y23 - y31) + (x31 - x23) * (y12 - y23)
398 # = (x31 - x23) * (y12 - y23) - (x12 - x23) * (y31 - y23)
399 # x = (d * B2.x + K * Y12_23).fover(d)
400 # y = (d * B2.y - K * X12_23).fover(d)
401 x, y = _pierlotxy2(B2, -K, Y12_23, X12_23, (X31_23 * Y12_23 -
402 X12_23 * Y31_23))
403 else:
404 x, y, _ = B2.xyz3
405 return x, y, _zidw(x, y, useZ, B1, B2, B3)
408def pierlotx(point1, point2, point3, alpha1, alpha2, alpha3, useZ=False,
409 **Clas_and_kwds):
410 '''3-Point resection using U{Pierlot<http://www.Telecom.ULg.ac.BE/publi/
411 publications/pierlot/Pierlot2014ANewThree>}'s method C{ToTal} with
412 I{exact} limits for the (pseudo-)singularities.
414 @arg point1: First point (C{Cartesian}, L{Vector3d}, C{Vector3Tuple},
415 C{Vector4Tuple} or C{Vector2Tuple} if C{B{useZ}=False}).
416 @arg point2: Second point (C{Cartesian}, L{Vector3d}, C{Vector3Tuple},
417 C{Vector4Tuple} or C{Vector2Tuple} if C{B{useZ}=False}).
418 @arg point3: Third point (C{Cartesian}, L{Vector3d}, C{Vector3Tuple},
419 C{Vector4Tuple} or C{Vector2Tuple} if C{B{useZ}=False}).
420 @arg alpha1: Angle at B{C{point1}} (C{degrees}, counter-clockwise).
421 @arg alpha2: Angle at B{C{point2}} (C{degrees}, counter-clockwise).
422 @arg alpha3: Angle at B{C{point3}} (C{degrees}, counter-clockwise).
423 @kwarg useZ: If C{True}, interpolate the survey point's Z component,
424 otherwise use C{z=INT0} (C{bool}).
425 @kwarg Clas_and_kwds: Optional class C{B{Clas}=B{point1}.classof} to
426 return the survey point with optionally other B{C{Clas}}
427 keyword arguments to instantiate the survey point.
429 @return: The survey (or robot) point, an instance of B{C{Clas}} or
430 B{C{point1}}'s (sub-)class.
432 @raise ResectionError: Near-coincident, -colinear or -concyclic points or
433 invalid B{C{alpha1}}, B{C{alpha2}} or B{C{alpha3}}.
435 @raise TypeError: Invalid B{C{point1}}, B{C{point2}} or B{C{point3}}.
437 @see: I{Pierlot}'s C function U{triangulationPierlot2<http://www.Telecom.ULg.ac.BE/
438 triangulation/doc/total_8c_source.html>} and function L{pierlot}, L{cassini},
439 L{collins5} and L{tienstra7}.
440 '''
442 def _a_z_Bs(Bs, *alphas):
443 ds = map2(_umod_360, alphas) # 0 <= alphas < 360
444 ds, Bs = zip(*sorted(_zip(ds, Bs))) # unzip
445 for p, d, B in _zip(ds, _rotate(ds), Bs):
446 d -= p # a12 = a2 - a1, ...
447 z = isnear0(fabs(d) % _180_0)
448 yield d, z, B
450 def _cot(s, c): # I{exact} cotangent
451 try:
452 return (c / s) if c else _copysign_0_0(s)
453 except ZeroDivisionError:
454 raise ValueError(_or(_coincident_, _colinear_))
456 Bs = _B3(useZ, point1, point2, point3)
457 try:
458 Cs = [0] # pseudo-global, passing the exception Case
459 xyz = _pierlotx3(_a_z_Bs(Bs, alpha1, alpha2, alpha3),
460 useZ, _cot, Cs.append)
461 return _Clas(pierlotx, point1, Clas_and_kwds, *xyz)
463 except (TypeError, ValueError) as x:
464 raise ResectionError(point1=point1, point2=point2, point3=point3, C=Cs.pop(),
465 alpha1=alpha1, alpha2=alpha2, alpha3=alpha3, cause=x)
468def _pierlotx3(a_z_Bs, useZ, cot, Cs):
469 '''(INTERNAL) Core of L{pierlotx}.
470 '''
471 (a12, z12, B1), \
472 (a23, z23, B2), \
473 (a31, z31, B3) = a_z_Bs
474 if z12 and not z23:
475 Cs(1)
476 elif z23 and not z31:
477 Cs(2)
478 a23, B1, B2, B3 = a31, B2, B3, B1
479 elif z31 and not z12:
480 Cs(3)
481 a23, B2, B3 = a12, B3, B2
482 else:
483 Cs(4)
484 return _pierlot3(B1, B2, B3, a12, a23, useZ, cot)
486 x1_, y1_, _ = B1.minus(B3).xyz3
487 x2_, y2_, _ = B2.minus(B3).xyz3
489 K = _Fsumf_(y1_ * x2_, -x1_ * y2_)
490 if K:
491 cot23 = cot(*sincos2d(a23))
493 # x23 = x2_ + cot23 * y2_
494 # y23 = y2_ - cot23 * x2_
496 # x31 = x1_ + cot23 * y1_
497 # y31 = y1_ - cot23 * x1_
499 # x31 - x23 = x1_ + cot23 * y1_ - x2_ - cot23 * y2_
500 X31_23 = _Fsumf_(x1_, cot23 * y1_, -x2_, -cot23 * y2_)
501 # y31 - y23 = y1_ - cot23 * x1_ - y2_ + cot23 * x2_
502 Y31_23 = _Fsumf_(y1_, -cot23 * x1_, -y2_, cot23 * x2_)
504 # d = (x31 - x23) * (x2_ - x1_) + (y31 - y23) * (y2_ - y1_)
505 # x = (D * B3.x - K * Y31_23).fover(d)
506 # y = (D * B3.y + K * X31_23).fover(d)
507 x, y = _pierlotxy2(B3, K, Y31_23, X31_23, (X31_23 * _Fsumf_(x2_, -x1_) +
508 Y31_23 * _Fsumf_(y2_, -y1_)))
509 else:
510 x, y, _ = B3.xyz3
511 return x, y, _zidw(x, y, useZ, B1, B2, B3)
514def _pierlotxy2(B, K, X, Y, D):
515 '''(INTERNAL) Helper for C{_pierlot3} and C{_pierlotx3}.
516 '''
517 d = float(D)
518 if isnear0(d):
519 raise ValueError(_or(_coincident_, _colinear_, _concyclic_))
520 x = (D * B.x - K * X).fover(d)
521 y = (D * B.y + K * Y).fover(d)
522 return x, y
525def _rotate(xs, n=1):
526 '''Rotate list or tuple C{xs} by C{n} items, right if C{n > 0} else left.
527 '''
528 return xs[n:] + xs[:n]
531def snellius3(a, b, degC, alpha, beta):
532 '''Snellius' surveying using U{Snellius Pothenot<https://WikiPedia.org/wiki/Snellius–Pothenot_problem>}.
534 @arg a: Length of the triangle side between corners C{B} and C{C} and opposite of
535 triangle corner C{A} (C{scalar}, non-negative C{meter}, conventionally).
536 @arg b: Length of the triangle side between corners C{C} and C{A} and opposite of
537 triangle corner C{B} (C{scalar}, non-negative C{meter}, conventionally).
538 @arg degC: Angle at triangle corner C{C}, opposite triangle side C{c} (non-negative C{degrees}).
539 @arg alpha: Angle subtended by triangle side B{C{b}} (non-negative C{degrees}).
540 @arg beta: Angle subtended by triangle side B{C{a}} (non-negative C{degrees}).
542 @return: L{Survey3Tuple}C{(PA, PB, PC)} with distance from survey point C{P} to
543 each of the triangle corners C{A}, C{B} and C{C}, same units as triangle
544 sides B{C{a}}, B{C{b}} and B{C{c}}.
546 @raise TriangleError: Invalid B{C{a}}, B{C{b}} or B{C{degC}} or negative B{C{alpha}}
547 or B{C{beta}}.
549 @see: Function L{wildberger3}.
550 '''
551 try:
552 a, b, degC, alpha, beta = t = map1(float, a, b, degC, alpha, beta)
553 if min(t) < 0:
554 raise ValueError(_negative_)
555 ra, rb, rC = map1(radians, alpha, beta, degC)
557 r = fsum1f_(ra, rb, rC) * _0_5
558 k = PI - r
559 if min(k, r) < 0:
560 raise ValueError(_or(_coincident_, _colinear_))
562 sa, sb = map1(sin, ra, rb)
563 p = atan2(sa * a, sb * b)
564 sp, cp, sr, cr = sincos2_(PI_4 - p, r)
565 p = atan2(sp * sr, cp * cr)
566 pa = k + p
567 pb = k - p
569 if fabs(sb) > fabs(sa):
570 pc = fabs(a * sin(pb) / sb)
571 elif sa:
572 pc = fabs(b * sin(pa) / sa)
573 else:
574 raise ValueError(_or(_colinear_, _coincident_))
576 pa = _triSide(b, pc, fsumf_(PI, -ra, -pa))
577 pb = _triSide(a, pc, fsumf_(PI, -rb, -pb))
578 return Survey3Tuple(pa, pb, pc, name=snellius3.__name__)
580 except (TypeError, ValueError) as x:
581 raise TriangleError(a=a, b=b, degC=degC, alpha=alpha, beta=beta, cause=x)
584def tienstra7(pointA, pointB, pointC, alpha, beta=None, gamma=None,
585 useZ=False, **Clas_and_kwds):
586 '''3-Point resection using U{Tienstra<https://WikiPedia.org/wiki/Tienstra_formula>}'s formula.
588 @arg pointA: First point (C{Cartesian}, L{Vector3d}, C{Vector3Tuple}, C{Vector4Tuple} or
589 C{Vector2Tuple} if C{B{useZ}=False}).
590 @arg pointB: Second point (C{Cartesian}, L{Vector3d}, C{Vector3Tuple}, C{Vector4Tuple} or
591 C{Vector2Tuple} if C{B{useZ}=False}).
592 @arg pointC: Third point (C{Cartesian}, L{Vector3d}, C{Vector3Tuple}, C{Vector4Tuple} or
593 C{Vector2Tuple} if C{B{useZ}=False}).
594 @arg alpha: Angle subtended by triangle side C{a} from B{C{pointB}} to B{C{pointC}}
595 (C{degrees}, non-negative).
596 @kwarg beta: Angle subtended by triangle side C{b} from B{C{pointA}} to B{C{pointC}}
597 (C{degrees}, non-negative) or C{None} if C{B{gamma} is not None}.
598 @kwarg gamma: Angle subtended by triangle side C{c} from B{C{pointA}} to B{C{pointB}}
599 (C{degrees}, non-negative) or C{None} if C{B{beta} is not None}.
600 @kwarg useZ: If C{True}, use and interpolate the Z component, otherwise force C{z=INT0}
601 (C{bool}).
602 @kwarg Clas_and_kwds: Optional class C{B{Clas}=B{pointA}.classof} to return the survey
603 point with optionally other B{C{Clas}} keyword arguments to instantiate
604 the survey point.
606 @note: Points B{C{pointA}}, B{C{pointB}} and B{C{pointC}} are ordered clockwise.
608 @return: L{Tienstra7Tuple}C{(pointP, A, B, C, a, b, c)} with survey C{pointP}, an
609 instance of B{C{Clas}} or B{C{pointA}}'s (sub-)class, with triangle angles C{A}
610 at B{C{pointA}}, C{B} at B{C{pointB}} and C{C} at B{C{pointC}} in C{degrees}
611 and with triangle sides C{a}, C{b} and C{c} in C{meter}, conventionally.
613 @raise ResectionError: Near-coincident, -colinear or -concyclic points or sum of
614 B{C{alpha}}, B{C{beta}} and B{C{gamma}} not C{360} or negative
615 B{C{alpha}}, B{C{beta}} or B{C{gamma}}.
617 @raise TypeError: Invalid B{C{pointA}}, B{C{pointB}} or B{C{pointC}}.
619 @see: U{3-Point Resection Solver<http://MesaMike.org/geocache/GC1B0Q9/tienstra/>},
620 U{V. Pierlot, M. Van Droogenbroeck, "A New Three Object Triangulation..."
621 <http://www.Telecom.ULg.ac.BE/publi/publications/pierlot/Pierlot2014ANewThree/>},
622 U{18 Triangulation Algorithms...<http://www.Telecom.ULg.ac.BE/triangulation/>} and
623 functions L{cassini}, L{collins5}, L{pierlot} and L{pierlotx}.
624 '''
626 def _deg_ks(r, s, ks, N):
627 if isnear0(fsumf_(PI, r, -s)): # r + (PI2 - s) == PI
628 raise ValueError(Fmt.PARENSPACED(concyclic=N))
629 # k = 1 / (cot(r) - cot(s))
630 # = 1 / (cos(r) / sin(r) - cos(s) / sin(s))
631 # = 1 / (cos(r) * sin(s) - cos(s) * sin(r)) / (sin(r) * sin(s))
632 # = sin(r) * sin(s) / (cos(r) * sin(s) - cos(s) * sin(r))
633 sr, cr, ss, cs = sincos2_(r, s)
634 c = fsum1f_(cr * ss, -cs * sr)
635 if isnear0(c):
636 raise ValueError(Fmt.PARENSPACED(cotan=N))
637 ks.append(sr * ss / c)
638 return Degrees(degrees(r), name=N) # C degrees
640 A, B, C = _ABC3(useZ, pointA, pointB, pointC)
641 try:
642 sa, sb, sc = map1(radians, alpha, (beta or 0), (gamma or 0))
643 if beta is None:
644 if gamma is None:
645 raise ValueError(_and(Fmt.EQUAL(beta=beta), Fmt.EQUAL(gamma=gamma)))
646 sb = fsumf_(PI2, -sa, -sc)
647 elif gamma is None:
648 sc = fsumf_(PI2, -sa, -sb)
649 else: # subtended angles must add to 360 degrees
650 r = fsum1f_(sa, sb, sc)
651 if fabs(r - PI2) > EPS:
652 raise ValueError(Fmt.EQUAL(sum=degrees(r)))
653 if min(sa, sb, sc) < 0:
654 raise ValueError(_negative_)
656 # triangle sides
657 a = B.minus(C).length
658 b = A.minus(C).length
659 c = A.minus(B).length
661 ks = [] # 3 Ks and triangle angles
662 dA = _deg_ks(_triAngle(b, c, a), sa, ks, _A_)
663 dB = _deg_ks(_triAngle(a, c, b), sb, ks, _B_)
664 dC = _deg_ks(_triAngle(a, b, c), sc, ks, _C_)
666 k = fsum1(ks)
667 if isnear0(k):
668 raise ValueError(Fmt.EQUAL(K=k))
669 x = Fdot(ks, A.x, B.x, C.x).fover(k)
670 y = Fdot(ks, A.y, B.y, C.y).fover(k)
671 z = _zidw(x, y, useZ, A, B, C)
673 P = _Clas(tienstra7, pointA, Clas_and_kwds, x, y, z)
674 return Tienstra7Tuple(P, dA, dB, dC, a, b, c, name=tienstra7.__name__)
676 except (TypeError, ValueError) as x:
677 raise ResectionError(pointA=pointA, pointB=pointB, pointC=pointC,
678 alpha=alpha, beta=beta, gamma=gamma, cause=x)
681def triAngle(a, b, c):
682 '''Compute one angle of a triangle.
684 @arg a: Adjacent triangle side length (C{scalar}, non-negative
685 C{meter}, conventionally).
686 @arg b: Adjacent triangle side length (C{scalar}, non-negative
687 C{meter}, conventionally).
688 @arg c: Opposite triangle side length (C{scalar}, non-negative
689 C{meter}, conventionally).
691 @return: Angle in C{radians} at triangle corner C{C}, opposite
692 triangle side B{C{c}}.
694 @raise TriangleError: Invalid or negative B{C{a}}, B{C{b}} or B{C{c}}.
696 @see: Functions L{triAngle5} and L{triSide}.
697 '''
698 try:
699 return _triAngle(a, b, c)
700 except (TypeError, ValueError) as x:
701 raise TriangleError(a=a, b=b, c=c, cause=x)
704def _triAngle(a, b, c):
705 # (INTERNAL) To allow callers to embellish errors
706 a, b, c = map1(float, a, b, c)
707 if a < b:
708 a, b = b, a
709 if b < 0 or c < 0:
710 raise ValueError(_negative_)
711 if a < EPS0:
712 raise ValueError(_coincident_)
713 b_a = b / a
714 if b_a < EPS0:
715 raise ValueError(_coincident_)
716 t = fsumf_(_1_0, b_a**2, -(c / a)**2) / (b_a * _2_0)
717 return acos1(t)
720def triAngle5(a, b, c):
721 '''Compute the angles of a triangle.
723 @arg a: Length of the triangle side opposite of triangle corner C{A}
724 (C{scalar}, non-negative C{meter}, conventionally).
725 @arg b: Length of the triangle side opposite of triangle corner C{B}
726 (C{scalar}, non-negative C{meter}, conventionally).
727 @arg c: Length of the triangle side opposite of triangle corner C{C}
728 (C{scalar}, non-negative C{meter}, conventionally).
730 @return: L{TriAngle5Tuple}C{(radA, radB, radC, rIn, area)} with angles
731 C{radA}, C{radB} and C{radC} at triangle corners C{A}, C{B}
732 and C{C}, all in C{radians}, the C{InCircle} radius C{rIn}
733 aka C{inradius}, same units as triangle sides B{C{a}},
734 B{C{b}} and B{C{c}} and the triangle C{area} in those same
735 units I{squared}.
737 @raise TriangleError: Invalid or negative B{C{a}}, B{C{b}} or B{C{c}}.
739 @see: Functions L{triAngle} and L{triArea}.
740 '''
741 try:
742 x, y, z = map1(float, a, b, c)
743 ab = x < y
744 if ab:
745 x, y = y, x
746 bc = y < z
747 if bc:
748 y, z = z, y
750 if z > EPS0: # z = min(a, b, c)
751 s = fsum1f_(z, y, x) * _0_5
752 sa, sb, r = (s - x), (s - y), (s - z)
753 r *= _over(sa * sb, s)
754 if r < EPS02:
755 raise ValueError(_coincident_)
756 r = sqrt(r)
757 rA = atan2(r, sa) * _2_0
758 rB = atan2(r, sb) * _2_0
759 rC = fsumf_(PI, -rA, -rB)
760 if min(rA, rB, rC) < 0:
761 raise ValueError(_colinear_)
762 s *= r # Heron's area
763 elif z < 0:
764 raise ValueError(_negative_)
765 else: # 0 <= c <= EPS0
766 rA = rB = PI_2
767 rC = r = s = _0_0
769 if bc:
770 rB, rC = rC, rB
771 if ab:
772 rA, rB = rB, rA
773 return TriAngle5Tuple(rA, rB, rC, r, s, name=triAngle5.__name__)
775 except (TypeError, ValueError) as x:
776 raise TriangleError(a=a, b=b, c=c, cause=x)
779def triArea(a, b, c):
780 '''Compute the area of a triangle using U{Heron's<https://
781 WikiPedia.org/wiki/Heron%27s_formula>} C{stable} formula.
783 @arg a: Length of the triangle side opposite of triangle corner C{A}
784 (C{scalar}, non-negative C{meter}, conventionally).
785 @arg b: Length of the triangle side opposite of triangle corner C{B}
786 (C{scalar}, non-negative C{meter}, conventionally).
787 @arg c: Length of the triangle side opposite of triangle corner C{C}
788 (C{scalar}, non-negative C{meter}, conventionally).
790 @return: The triangle area (C{float}, conventionally C{meter} or
791 same units as B{C{a}}, B{C{b}} and B{C{c}} I{squared}).
793 @raise TriangleError: Invalid or negative B{C{a}}, B{C{b}} or B{C{c}}.
794 '''
795 try:
796 r, y, x = sorted(map1(float, a, b, c))
797 if r > 0: # r = min(a, b, c)
798 ab = x - y
799 bc = y - r
800 y += r
801 r = (x + y) * (r - ab) * (r + ab) * (x + bc)
802 if r:
803 r = sqrt(r / _16_0)
804 elif r < 0:
805 raise ValueError(_negative_)
806 return r
808 except (TypeError, ValueError) as x:
809 raise TriangleError(a=a, b=b, c=c, cause=x)
812def triSide(a, b, radC):
813 '''Compute one side of a triangle.
815 @arg a: Adjacent triangle side length (C{scalar},
816 non-negative C{meter}, conventionally).
817 @arg b: Adjacent triangle side length (C{scalar},
818 non-negative C{meter}, conventionally).
819 @arg radC: Angle included by sides B{C{a}} and B{C{b}},
820 opposite triangle side C{c} (C{radians}).
822 @return: Length of triangle side C{c}, opposite triangle
823 corner C{C} and angle B{C{radC}}, same units as
824 B{C{a}} and B{C{b}}.
826 @raise TriangleError: Invalid B{C{a}}, B{C{b}} or B{C{radC}}.
828 @see: Functions L{sqrt_a}, L{triAngle}, L{triSide2} and L{triSide4}.
829 '''
830 try:
831 return _triSide(a, b, radC)
832 except (TypeError, ValueError) as x:
833 raise TriangleError(a=a, b=b, radC=radC, cause=x)
836def _triSide(a, b, radC):
837 # (INTERNAL) To allow callers to embellish errors
838 a, b, r = t = map1(float, a, b, radC)
839 if min(t) < 0:
840 raise ValueError(_negative_)
842 if a < b:
843 a, b = b, a
844 if a > EPS0:
845 ba = b / a
846 c2 = fsumf_(_1_0, ba**2, _N_2_0 * ba * cos(r))
847 if c2 > EPS02:
848 return a * sqrt(c2)
849 elif c2 < 0:
850 raise ValueError(_invalid_)
851 return hypot(a, b)
854def triSide2(b, c, radB):
855 '''Compute a side and its opposite angle of a triangle.
857 @arg b: Adjacent triangle side length (C{scalar},
858 non-negative C{meter}, conventionally).
859 @arg c: Adjacent triangle side length (C{scalar},
860 non-negative C{meter}, conventionally).
861 @arg radB: Angle included by sides B{C{a}} and B{C{c}},
862 opposite triangle side C{b} (C{radians}).
864 @return: L{TriSide2Tuple}C{(a, radA)} with triangle angle
865 C{radA} in C{radians} and length of the opposite
866 triangle side C{a}, same units as B{C{b}} and B{C{c}}.
868 @raise TriangleError: Invalid B{C{b}} or B{C{c}} or either
869 B{C{b}} or B{C{radB}} near zero.
871 @see: Functions L{sqrt_a}, L{triSide} and L{triSide4}.
872 '''
873 try:
874 return _triSide2(b, c, radB)
875 except (TypeError, ValueError) as x:
876 raise TriangleError(b=b, c=c, radB=radB, cause=x)
879def _triSide2(b, c, radB):
880 # (INTERNAL) To allow callers to embellish errors
881 b, c, rB = map1(float, b, c, radB)
882 if min(b, c, rB) < 0:
883 raise ValueError(_negative_)
884 sB, cB = sincos2(rB)
885 if isnear0(sB):
886 if not isnear0(b):
887 raise ValueError(_invalid_)
888 a, rA = ((b + c), PI) if cB < 0 else (fabs(b - c), _0_0)
889 elif isnear0(b):
890 raise ValueError(_invalid_)
891 else:
892 rA = fsumf_(PI, -rB, -asin1(c * sB / b))
893 a = sin(rA) * b / sB
894 return TriSide2Tuple(a, rA, name=triSide2.__name__)
897def triSide4(radA, radB, c):
898 '''Compute two sides and the height of a triangle.
900 @arg radA: Angle at triangle corner C{A}, opposite triangle side C{a}
901 (non-negative C{radians}).
902 @arg radB: Angle at triangle corner C{B}, opposite triangle side C{b}
903 (non-negative C{radians}).
904 @arg c: Length of triangle side between triangle corners C{A} and C{B},
905 (C{scalar}, non-negative C{meter}, conventionally).
907 @return: L{TriSide4Tuple}C{(a, b, radC, d)} with triangle sides C{a} and
908 C{b} and triangle height C{d} perpendicular to triangle side
909 B{C{c}}, all in the same units as B{C{c}} and interior angle
910 C{radC} in C{radians} at triangle corner C{C}, opposite
911 triangle side B{C{c}}.
913 @raise TriangleError: Invalid or negative B{C{radA}}, B{C{radB}} or B{C{c}}.
915 @see: U{Triangulation, Surveying<https://WikiPedia.org/wiki/Triangulation_(surveying)>}
916 and functions L{sqrt_a}, L{triSide} and L{triSide2}.
917 '''
918 try:
919 rA, rB, c = map1(float, radA, radB, c)
920 rC = fsumf_(PI, -rA, -rB)
921 if min(rC, rA, rB, c) < 0:
922 raise ValueError(_negative_)
923 sa, ca, sb, cb = sincos2_(rA, rB)
924 sc = fsum1f_(sa * cb, sb * ca)
925 if sc < EPS0 or min(sa, sb) < 0:
926 raise ValueError(_invalid_)
927 sc = c / sc
928 return TriSide4Tuple((sa * sc), (sb * sc), rC, (sa * sb * sc),
929 name=triSide4.__name__)
931 except (TypeError, ValueError) as x:
932 raise TriangleError(radA=radA, radB=radB, c=c, cause=x)
935def wildberger3(a, b, c, alpha, beta, R3=min):
936 '''Snellius' surveying using U{Rational Trigonometry
937 <https://WikiPedia.org/wiki/Snellius–Pothenot_problem>}.
939 @arg a: Length of the triangle side between corners C{B} and C{C} and opposite of
940 triangle corner C{A} (C{scalar}, non-negative C{meter}, conventionally).
941 @arg b: Length of the triangle side between corners C{C} and C{A} and opposite of
942 triangle corner C{B} (C{scalar}, non-negative C{meter}, conventionally).
943 @arg c: Length of the triangle side between corners C{A} and C{B} and opposite of
944 triangle corner C{C} (C{scalar}, non-negative C{meter}, conventionally).
945 @arg alpha: Angle subtended by triangle side B{C{b}} (C{degrees}, non-negative).
946 @arg beta: Angle subtended by triangle side B{C{a}} (C{degrees}, non-negative).
947 @kwarg R3: Callable to determine C{R3} from C{(R3 - C)**2 = D}, typically standard
948 Python function C{min} or C{max}, invoked with 2 arguments.
950 @return: L{Survey3Tuple}C{(PA, PB, PC)} with distance from survey point C{P} to
951 each of the triangle corners C{A}, C{B} and C{C}, same units as B{C{a}},
952 B{C{b}} and B{C{c}}.
954 @raise TriangleError: Invalid B{C{a}}, B{C{b}} or B{C{c}} or negative B{C{alpha}} or
955 B{C{beta}} or B{C{R3}} not C{callable}.
957 @see: U{Wildberger, Norman J.<https://Math.Sc.Chula.ac.TH/cjm/content/
958 survey-article-greek-geometry-rational-trigonometry-and-snellius-–-pothenot-surveying>},
959 U{Devine Proportions, page 252<http://www.MS.LT/derlius/WildbergerDivineProportions.pdf>}
960 and function L{snellius3}.
961 '''
962 def _s(x):
963 return sin(x)**2
965 def _vpa(r3, q2, q3, s2, s3):
966 r1 = s2 * q3 / s3
967 r = r1 * r3 * _4_0
968 n = (r - _Fsumf_(r1, r3, -q2)**2).fover(s3)
969 if n < 0 or r < EPS0:
970 raise ValueError(_coincident_)
971 return sqrt((n / r) * q3) if n else _0_0
973 try:
974 a, b, c, da, db = q = map1(float, a, b, c, alpha, beta)
975 if min(q) < 0:
976 raise ValueError(_negative_)
978 q1, q2, q3 = q = a**2, b**2, c**2
979 if min(q) < EPS02:
980 raise ValueError(_coincident_)
982 ra, rb = map1(radians, da, db)
983 s1, s2, s3 = s = map1(_s, rb, ra, ra + rb) # rb, ra!
984 if min(s) < EPS02:
985 raise ValueError(_or(_coincident_, _colinear_))
987 q4 = hypot2_(*q) * _2_0 # a**4 + ...
988 Qs = _Fsumf_(*q) # == hypot2_(a, b, c)
989 d0 = (Qs**2 - q4).fmul(s1 * s2).fover(s3)
990 if d0 < 0:
991 raise ValueError(_negative_)
992 s += _Fsumf_(*s), # == fsum1(s),
993 C0 = Fdot(s, q1, q2, q3, -Qs * _0_5)
994 r3 = C0.fover(-s3) # C0 /= -s3
995 if d0 > EPS02: # > c0
996 _xcallable(R3=R3)
997 d0 = sqrt(d0)
998 r3 = R3(float(C0 + d0), float(C0 - d0)) # XXX min or max
1000 pa = _vpa(r3, q2, q3, s2, s3)
1001 pb = _vpa(r3, q1, q3, s1, s3)
1002 pc = favg(_triSide2(b, pa, ra).a,
1003 _triSide2(a, pb, rb).a)
1004 return Survey3Tuple(pa, pb, pc, name=wildberger3.__name__)
1006 except (TypeError, ValueError) as x:
1007 raise TriangleError(a=a, b=b, c=c, alpha=alpha, beta=beta, R3=R3, cause=x)
1010def _zidw(x, y, useZ, *ABC):
1011 if useZ: # interpolate z or coplanar with A, B and C?
1012 t = tuple(_.z for _ in ABC)
1013 v = Vector3d(x, y, fmean(t))
1014 z = fidw(t, (v.minus(T).length for T in ABC))
1015 else:
1016 z = INT0
1017 return z
1019# **) MIT License
1020#
1021# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.
1022#
1023# Permission is hereby granted, free of charge, to any person obtaining a
1024# copy of this software and associated documentation files (the "Software"),
1025# to deal in the Software without restriction, including without limitation
1026# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1027# and/or sell copies of the Software, and to permit persons to whom the
1028# Software is furnished to do so, subject to the following conditions:
1029#
1030# The above copyright notice and this permission notice shall be included
1031# in all copies or substantial portions of the Software.
1032#
1033# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1034# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1035# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
1036# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1037# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1038# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1039# OTHER DEALINGS IN THE SOFTWARE.