Coverage for pygeodesy/resections.py: 97%

360 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-08-06 15:27 -0400

1 

2# -*- coding: utf-8 -*- 

3 

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{triAngle4}, L{triSide}, L{triSide2} and L{triSide4}. 

7 

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 

14 

15from pygeodesy.basics import map1, _ALL_LAZY 

16from pygeodesy.constants import EPS, EPS0, EPS02, INT0, NEG0, PI, PI2, PI_2, PI_4, \ 

17 isnear0, _umod_360, _0_0, _0_5, _1_0, _N_1_0, _2_0, \ 

18 _N_2_0, _4_0, _180_0, _360_0 

19from pygeodesy.errors import _and, _or, TriangleError, _ValueError, _xkwds 

20from pygeodesy.fmath import favg, Fdot, fidw, fmean, hypot, hypot2_ 

21from pygeodesy.fsums import Fsum, fsumf_, fsum1, fsum1f_ 

22from pygeodesy.interns import _a_, _A_, _b_, _B_, _c_, _C_, _coincident_, _colinear_, \ 

23 _d_, _eps_, _invalid_, _negative_, _not_, _rIn_, _SPACE_ 

24# from pygeodesy.lazily import _ALL_LAZY # from .basics 

25from pygeodesy.named import Fmt, _NamedTuple, _Pass 

26# from pygeodesy.streprs import Fmt # from .named 

27from pygeodesy.units import Degrees, Distance, Radians 

28from pygeodesy.utily import acos1, asin1, sincos2, sincos2_, sincos2d, sincos2d_ 

29from pygeodesy.vector3d import _otherV3d, Vector3d 

30 

31from math import cos, atan2, degrees, fabs, radians, sin, sqrt 

32 

33__all__ = _ALL_LAZY.resections 

34__version__ = '23.06.07' 

35 

36_concyclic_ = 'concyclic' 

37_PA_ = 'PA' 

38_PB_ = 'PB' 

39_PC_ = 'PC' 

40_pointH_ = 'pointH' 

41_pointP_ = 'pointP' 

42_positive_ = 'positive' 

43_R3__ = 'R3 ' 

44_radA_ = 'radA' 

45_radB_ = 'radB' 

46_radC_ = 'radC' 

47 

48 

49class Collins5Tuple(_NamedTuple): 

50 '''5-Tuple C{(pointP, pointH, a, b, c)} with survey C{pointP}, auxiliary 

51 C{pointH}, each an instance of B{C{pointA}}'s (sub-)class and triangle 

52 sides C{a}, C{b} and C{c} in C{meter}, conventionally. 

53 ''' 

54 _Names_ = (_pointP_, _pointH_, _a_, _b_, _c_) 

55 _Units_ = (_Pass, _Pass, Distance, Distance, Distance) 

56 

57 

58class ResectionError(_ValueError): 

59 '''Error raised for issues in L{pygeodesy.resections}. 

60 ''' 

61 pass 

62 

63 

64class Survey3Tuple(_NamedTuple): 

65 '''3-Tuple C{(PA, PB, PC)} with distance from survey point C{P} to each of 

66 the triangle corners C{A}, C{B} and C{C} in C{meter}, conventionally. 

67 ''' 

68 _Names_ = (_PA_, _PB_, _PC_) 

69 _Units_ = ( Distance, Distance, Distance) 

70 

71 

72class Tienstra7Tuple(_NamedTuple): 

73 '''7-Tuple C{(pointP, A, B, C, a, b, c)} with survey C{pointP}, interior 

74 triangle angles C{A}, C{B} and C{C} in C{degrees} and triangle sides 

75 C{a}, C{b} and C{c} in C{meter}, conventionally. 

76 ''' 

77 _Names_ = (_pointP_, _A_, _B_, _C_, _a_, _b_, _c_) 

78 _Units_ = (_Pass, Degrees, Degrees, Degrees, Distance, Distance, Distance) 

79 

80 

81class TriAngle4Tuple(_NamedTuple): 

82 '''4-Tuple C{(radA, radB, radC, rIn)} with the interior angles at triangle 

83 corners C{A}, C{B} and C{C} in C{radians} and the C{InCircle} radius 

84 C{rIn} aka C{inradius} in C{meter}, conventionally. 

85 ''' 

86 _Names_ = (_radA_, _radB_, _radC_, _rIn_) 

87 _Units_ = ( Radians, Radians, Radians, Distance) 

88 

89 

90class TriSide2Tuple(_NamedTuple): 

91 '''2-Tuple C{(a, radA)} with triangle side C{a} in C{meter}, conventionally 

92 and angle C{radA} at the opposite triangle corner in C{radians}. 

93 ''' 

94 _Names_ = (_a_, _radA_) 

95 _Units_ = ( Distance, Radians) 

96 

97 

98class TriSide4Tuple(_NamedTuple): 

99 '''4-Tuple C{(a, b, radC, d)} with interior angle C{radC} at triangle corner 

100 C{C} in C{radians}with the length of triangle sides C{a} and C{b} and 

101 with triangle height C{d} perpendicular to triangle side C{c}, in the 

102 same units as triangle sides C{a} and C{b}. 

103 ''' 

104 _Names_ = (_a_, _b_, _radC_, _d_) 

105 _Units_ = ( Distance, Distance, Radians, Distance) 

106 

107 

108def _ABC3(useZ, pointA, pointB, pointC): 

109 '''(INTERNAL) Helper for L{cassini} and L{tienstra7}. 

110 ''' 

111 return (_otherV3d(useZ=useZ, pointA=pointA), 

112 _otherV3d(useZ=useZ, pointB=pointB), 

113 _otherV3d(useZ=useZ, pointC=pointC)) 

114 

115 

116def _B3(useZ, point1, point2, point3): 

117 '''(INTERNAL) Helper for L{pierlot} and L{pierlotx}. 

118 ''' 

119 return (_otherV3d(useZ=useZ, point1=point1), 

120 _otherV3d(useZ=useZ, point2=point2), 

121 _otherV3d(useZ=useZ, point3=point3)) 

122 

123 

124def cassini(pointA, pointB, pointC, alpha, beta, useZ=False, Clas=None, **Clas_kwds): 

125 '''3-Point resection using U{Cassini<https://NL.WikiPedia.org/wiki/Achterwaartse_insnijding>}'s method. 

126 

127 @arg pointA: First point (C{Cartesian}, L{Vector3d}, C{Vector3Tuple}, 

128 C{Vector4Tuple} or C{Vector2Tuple} if C{B{useZ}=False}). 

129 @arg pointB: Second point (C{Cartesian}, L{Vector3d}, C{Vector3Tuple}, 

130 C{Vector4Tuple} or C{Vector2Tuple} if C{B{useZ}=False}). 

131 @arg pointC: Center point (C{Cartesian}, L{Vector3d}, C{Vector3Tuple}, 

132 C{Vector4Tuple} or C{Vector2Tuple} if C{B{useZ}=False}). 

133 @arg alpha: Angle subtended by triangle side B{C{pointA}} to B{C{pointC}} 

134 (C{degrees}, non-negative). 

135 @arg beta: Angle subtended by triangle side B{C{pointB}} to B{C{pointC}} 

136 (C{degrees}, non-negative). 

137 @kwarg useZ: If C{True}, use and interpolate the Z component, otherwise 

138 force C{z=INT0} (C{bool}). 

139 @kwarg Clas: Optional class to return the survey point or C{None} for 

140 B{C{pointA}}'s (sub-)class. 

141 @kwarg Clas_kwds: Optional, additional keyword argument for the survey 

142 point instance. 

143 

144 @note: Typically, B{C{pointC}} is between B{C{pointA}} and B{C{pointB}}. 

145 

146 @return: The survey point, an instance of B{C{Clas}} or if C{B{Clas} is 

147 None} an instance of B{C{pointA}}'s (sub-)class. 

148 

149 @raise ResectionError: Near-coincident, -colinear or -concyclic points 

150 or negative or invalid B{C{alpha}} or B{C{beta}}. 

151 

152 @raise TypeError: Invalid B{C{pointA}}, B{C{pointB}} or B{C{pointM}}. 

153 

154 @see: U{Three Point Resection Problem<https://Dokumen.tips/documents/ 

155 three-point-resection-problem-introduction-kaestner-burkhardt-method.html>} 

156 and functions L{pygeodesy.collins5}, L{pygeodesy.pierlot} and 

157 L{pygeodesy.tienstra7}. 

158 ''' 

159 

160 def _H(A, C, sa): 

161 s, c = sincos2d(sa) 

162 if isnear0(s): 

163 raise ValueError(_or(_coincident_, _colinear_)) 

164 t = s, c, c 

165 x = Fdot(t, A.x, C.y, -A.y).fover(s) 

166 y = Fdot(t, A.y, -C.x, A.x).fover(s) 

167 return x, y 

168 

169 A, B, C = _ABC3(useZ, pointA, pointB, pointC) 

170 try: 

171 sa, sb = map1(float, alpha, beta) 

172 if min(sa, sb) < 0: 

173 raise ValueError(_negative_) 

174 if fsumf_(_360_0, -sa, -sb) < EPS0: 

175 raise ValueError() 

176 

177 x1, y1 = _H(A, C, sa) 

178 x2, y2 = _H(B, C, -sb) 

179 

180 x = x1 - x2 

181 y = y1 - y2 

182 if isnear0(x) or isnear0(y): 

183 raise ValueError(_SPACE_(_concyclic_, (x, y))) 

184 

185 m = y / x 

186 n = x / y 

187 N = n + m 

188 if isnear0(N): 

189 raise ValueError(_SPACE_(_concyclic_, (m, n, N))) 

190 

191 t = n, m, _1_0, _N_1_0 

192 x = Fdot(t, C.x, x1, C.y, y1).fover(N) 

193 y = Fdot(t, y1, C.y, C.x, x1).fover(N) 

194 z = _zidw(x, y, useZ, A, B, C) 

195 

196 clas = Clas or pointA.classof 

197 return clas(x, y, z, **_xkwds(Clas_kwds, name=cassini.__name__)) 

198 

199 except (TypeError, ValueError) as x: 

200 raise ResectionError(pointA=pointA, pointB=pointB, pointC=pointC, 

201 alpha=alpha, beta=beta, cause=x) 

202 

203 

204def collins5(pointA, pointB, pointC, alpha, beta, useZ=False, Clas=None, **Clas_kwds): 

205 '''3-Point resection using U{Collins<https://Dokumen.tips/documents/ 

206 three-point-resection-problem-introduction-kaestner-burkhardt-method.html>}' method. 

207 

208 @arg pointA: First point (C{Cartesian}, L{Vector3d}, C{Vector3Tuple}, 

209 C{Vector4Tuple} or C{Vector2Tuple} if C{B{useZ}=False}). 

210 @arg pointB: Second point (C{Cartesian}, L{Vector3d}, C{Vector3Tuple}, 

211 C{Vector4Tuple} or C{Vector2Tuple} if C{B{useZ}=False}). 

212 @arg pointC: Center point (C{Cartesian}, L{Vector3d}, C{Vector3Tuple}, 

213 C{Vector4Tuple} or C{Vector2Tuple} if C{B{useZ}=False}). 

214 @arg alpha: Angle subtended by triangle side C{b} from B{C{pointA}} to 

215 B{C{pointC}} (C{degrees}, non-negative). 

216 @arg beta: Angle subtended by triangle side C{a} from B{C{pointB}} to 

217 B{C{pointC}} (C{degrees}, non-negative). 

218 @kwarg useZ: If C{True}, use and interpolate the Z component, otherwise 

219 force C{z=INT0} (C{bool}). 

220 @kwarg Clas: Optional class to return the survey and auxiliary point 

221 or C{None} for B{C{pointA}}'s (sub-)class. 

222 @kwarg Clas_kwds: Optional, additional keyword argument for the survey 

223 and auxiliary point instance. 

224 

225 @note: Typically, B{C{pointC}} is between B{C{pointA}} and B{C{pointB}}. 

226 

227 @return: L{Collins5Tuple}C{(pointP, pointH, a, b, c)} with survey C{pointP}, 

228 auxiliary C{pointH}, each an instance of B{C{Clas}} or if C{B{Clas} 

229 is None} an instance of B{C{pointA}}'s (sub-)class and triangle 

230 sides C{a}, C{b} and C{c} in C{meter}, conventionally. 

231 

232 @raise ResectionError: Near-coincident, -colinear or -concyclic points 

233 or negative or invalid B{C{alpha}} or B{C{beta}}. 

234 

235 @raise TypeError: Invalid B{C{pointA}}, B{C{pointB}} or B{C{pointM}}. 

236 

237 @see: U{Collins' methode<https://NL.WikiPedia.org/wiki/Achterwaartse_insnijding>} 

238 and functions L{pygeodesy.cassini}, L{pygeodesy.pierlot} and 

239 L{pygeodesy.tienstra7}. 

240 ''' 

241 

242 def _azi_len2(A, B, pi2): 

243 v = B.minus(A) 

244 r = atan2(v.x, v.y) 

245 if pi2 and r < 0: 

246 r += pi2 

247 return r, v.length 

248 

249 def _cV3(d, r, A, B, C, useZ, V3, **kwds): 

250 s, c = sincos2(r) 

251 x = A.x + d * s 

252 y = A.y + d * c 

253 z = _zidw(x, y, useZ, A, B, C) 

254 return V3(x, y, z, **kwds) 

255 

256 A, B, C = _ABC3(useZ, pointA, pointB, pointC) 

257 try: 

258 ra, rb = radians(alpha), radians(beta) 

259 if min(ra, rb) < 0: 

260 raise ValueError(_negative_) 

261 

262 sra, srH = sin(ra), sin(ra + rb - PI) # rH = PI - ((PI - ra) + (PI - rb)) 

263 if isnear0(sra) or isnear0(srH): 

264 raise ValueError(_or(_coincident_, _colinear_, _concyclic_)) 

265 

266 clas = Clas or pointA.classof 

267 kwds = _xkwds(Clas_kwds, name=collins5.__name__) 

268 

269# za, a = _azi_len2(C, B, PI2) 

270 zb, b = _azi_len2(C, A, PI2) 

271 zc, c = _azi_len2(A, B, 0) 

272 

273# d = c * sin(PI - rb) / srH # B.minus(H).length 

274 d = c * sin(PI - ra) / srH # A.minus(H).length 

275 r = zc + PI - rb # zh = zc + (PI - rb) 

276 H = _cV3(d, r, A, B, C, useZ, Vector3d) 

277 

278 zh, _ = _azi_len2(C, H, PI2) 

279 

280# d = a * sin(za - zh) / sin(rb) # B.minus(P).length 

281 d = b * sin(zb - zh) / sra # A.minus(P).length 

282 r = zh - ra # zb - PI + (PI - ra - (zb - zh)) 

283 P = _cV3(d, r, A, B, C, useZ, clas, **kwds) 

284 

285 H = clas(H.x, H.y, H.z, **kwds) 

286 a = B.minus(C).length 

287 

288 return Collins5Tuple(P, H, a, b, c, name=collins5.__name__) 

289 

290 except (TypeError, ValueError) as x: 

291 raise ResectionError(pointA=pointA, pointB=pointB, pointC=pointC, 

292 alpha=alpha, beta=beta, cause=x) 

293 

294 

295def pierlot(point1, point2, point3, alpha12, alpha23, useZ=False, eps=EPS, 

296 Clas=None, **Clas_kwds): 

297 '''3-Point resection using U{Pierlot<http://www.Telecom.ULg.ac.BE/publi/publications/ 

298 pierlot/Pierlot2014ANewThree>}'s method C{ToTal} with I{approximate} limits for 

299 the (pseudo-)singularities. 

300 

301 @arg point1: First point (C{Cartesian}, L{Vector3d}, C{Vector3Tuple}, 

302 C{Vector4Tuple} or C{Vector2Tuple} if C{B{useZ}=False}). 

303 @arg point2: Second point (C{Cartesian}, L{Vector3d}, C{Vector3Tuple}, 

304 C{Vector4Tuple} or C{Vector2Tuple} if C{B{useZ}=False}). 

305 @arg point3: Third point (C{Cartesian}, L{Vector3d}, C{Vector3Tuple}, 

306 C{Vector4Tuple} or C{Vector2Tuple} if C{B{useZ}=False}). 

307 @arg alpha12: Angle subtended from B{C{point1}} to B{C{point2}} or 

308 B{C{alpha2 - alpha1}} (C{degrees}). 

309 @arg alpha23: Angle subtended from B{C{point2}} to B{C{point3}} or 

310 B{C{alpha3 - alpha2}}(C{degrees}). 

311 @kwarg useZ: If C{True}, interpolate the survey point's Z component, 

312 otherwise use C{z=INT0} (C{bool}). 

313 @kwarg eps: Tolerance for C{cot} (pseudo-)singularities (C{float}). 

314 @kwarg Clas: Optional class to return the survey point, if C{None} use 

315 B{C{point1}}'s (sub-)class. 

316 @kwarg Clas_kwds: Optional, additional keyword arguments for the survey 

317 point instance. 

318 

319 @note: Typically, B{C{point1}}, B{C{point2}} and B{C{point3}} are ordered 

320 by angle, modulo 360, counter-clockwise. 

321 

322 @return: The survey (or robot) point, an instance of B{C{Clas}} or if 

323 C{B{Clas} is None} an instance of B{C{point1}}'s (sub-)class. 

324 

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}}. 

328 

329 @raise TypeError: Invalid B{C{point1}}, B{C{point2}} or B{C{point3}}. 

330 

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{pygeodesy.pierlotx}, L{pygeodesy.cassini}, 

338 L{pygeodesy.collins5} and L{pygeodesy.tienstra7}. 

339 ''' 

340 

341 def _cot(s, c): # -eps < I{approximate} cotangent < eps 

342 if eps > 0: 

343 return c / (min(s, -eps) if s < 0 else max(s, eps)) 

344 raise ValueError(_SPACE_(_eps_, _not_, _positive_)) 

345 

346 B1, B2, B3 = _B3(useZ, point1, point2, point3) 

347 try: 

348 x, y, z = _pierlot3(B1, B2, B3, alpha12, alpha23, useZ, _cot) 

349 clas = Clas or point1.classof 

350 return clas(x, y, z, **_xkwds(Clas_kwds, name=pierlot.__name__)) 

351 

352 except (TypeError, ValueError) as x: 

353 raise ResectionError(point1=point1, point2=point2, point3=point3, 

354 alpha12=alpha12, alpha23=alpha23, eps=eps, cause=x) 

355 

356 

357def _pierlot3(B1, B2, B3, a12, a23, useZ, cot): 

358 '''(INTERNAL) Shared L{pierlot} and L{pierlotx}. 

359 ''' 

360 x1_, y1_, _ = B1.minus(B2).xyz 

361 x3_, y3_, _ = B3.minus(B2).xyz 

362 

363 s12, c12, s23, c23 = sincos2d_(a12, a23) 

364 # cot31 = (1 - cot12 * cot23) / (cot12 + cot32) 

365 # = (1 - c12 / s12 * c23 / s23) / (c12 / s12 + c23 / s23) 

366 # = (1 - (c12 * c23) / (s12 * s23)) / (c12 * s23 + s12 * c23) / (s12 * s23) 

367 # = (s12 * s23 - c12 * c23) / (c12 * s23 + s12 * c23) 

368 # = c31 / s31 

369 cot31 = cot(fsum1f_(c12 * s23, s12 * c23), # s31 

370 fsum1f_(s12 * s23, -c12 * c23)) # c31 

371 

372 K = Fsum(x3_ * x1_, cot31 * (y3_ * x1_), 

373 y3_ * y1_, -cot31 * (x3_ * y1_)) 

374 if K: 

375 cot12 = cot(s12, c12) 

376 cot23 = cot(s23, c23) 

377 

378 # x12 = x1_ + cot12 * y1_ 

379 # y12 = y1_ - cot12 * x1_ 

380 

381 # x23 = x3_ - cot23 * y3_ 

382 # y23 = y3_ + cot23 * x3_ 

383 

384 # x31 = x3_ + x1_ + cot31 * (y3_ - y1_) 

385 # y31 = y3_ + y1_ - cot31 * (x3_ - x1_) 

386 

387 # x12 - x23 = x1_ + cot12 * y1_ - x3_ + cot23 * y3_ 

388 X12_23 = Fsum(x1_, cot12 * y1_, -x3_, cot23 * y3_) 

389 # y12 - y23 = y1_ - cot12 * x1_ - y3_ - cot23 * x3_ 

390 Y12_23 = Fsum(y1_, -cot12 * x1_, -y3_, -cot23 * x3_) 

391 

392 # x31 - x23 = x3_ + x1_ + cot31 * (y3_ - y1_) - x3_ + cot23 * y3_ 

393 # = x1_ + cot31 * y3_ - cot31 * y1_ + cot23 * y3_ 

394 X31_23 = Fsum(x1_, -cot31 * y1_, cot31 * y3_, cot23 * y3_) 

395 # y31 - y23 = y3_ + y1_ - cot31 * (x3_ - x1_) - y3_ - cot23 * x3_ 

396 # = y1_ - cot31 * x3_ + cot31 * x1_ - cot23 * x3_ 

397 Y31_23 = Fsum(y1_, cot31 * x1_, -cot31 * x3_, -cot23 * x3_) 

398 

399 # d = (x12 - x23) * (y23 - y31) + (x31 - x23) * (y12 - y23) 

400 # = (x31 - x23) * (y12 - y23) - (x12 - x23) * (y12 - y23) 

401 d = float(X31_23 * Y12_23 - X12_23 * Y31_23) 

402 if isnear0(d): 

403 raise ValueError(_or(_coincident_, _colinear_, _concyclic_)) 

404 

405 x = (B2.x * d + K * Y12_23).fover(d) 

406 y = (B2.y * d - K * X12_23).fover(d) 

407 else: 

408 x, y, _ = B2.xyz 

409 return x, y, _zidw(x, y, useZ, B1, B2, B3) 

410 

411 

412def _pierlotx3(a_z_Bs, useZ, cot, C): 

413 '''(INTERNAL) Core of L{pierlotx}. 

414 ''' 

415 (a12, z12, B1), \ 

416 (a23, z23, B2), \ 

417 (a31, z31, B3) = a_z_Bs 

418 if z12 and not z23: 

419 C(1) 

420 elif z23 and not z31: 

421 C(2) 

422 a23, B1, B2, B3 = a31, B2, B3, B1 

423 elif z31 and not z12: 

424 C(3) 

425 a23, B2, B3 = a12, B3, B2 

426 else: 

427 C(4) 

428 return _pierlot3(B1, B2, B3, a12, a23, useZ, cot) 

429 

430 x1_, y1_, _ = B1.minus(B3).xyz 

431 x2_, y2_, _ = B2.minus(B3).xyz 

432 

433 K = Fsum(_1_0, y1_ * x2_, -x1_ * y2_, _N_1_0) # 1-primed 

434 if K: 

435 cot23 = cot(*sincos2d(a23)) 

436 

437 # x23 = x2_ + cot23 * y2_ 

438 # y23 = y2_ - cot23 * x2_ 

439 

440 # x31 = x1_ + cot23 * y1_ 

441 # y31 = y1_ - cot23 * x1_ 

442 

443 # x31 - x23 = x1_ + cot23 * y1_ - x2_ - cot23 * y2_ 

444 X31_23 = Fsum(x1_, cot23 * y1_, -x2_, -cot23 * y2_) 

445 # y31 - y23 = y1_ - cot23 * x1_ - y2_ + cot23 * x2_ 

446 Y31_23 = Fsum(y1_, -cot23 * x1_, -y2_, cot23 * x2_) 

447 

448 # d = (x31 - x23) * (x2_ - x1_) - (y31 - y23) * (y1_ - y2_) 

449 d = float(X31_23 * x2_ - X31_23 * x1_ + 

450 Y31_23 * y2_ - Y31_23 * y1_) 

451 if isnear0(d): 

452 raise ValueError(_or(_coincident_, _colinear_, _concyclic_)) 

453 

454 x = (B3.x * d - K * Y31_23).fover(d) 

455 y = (B3.y * d + K * X31_23).fover(d) 

456 else: 

457 x, y, _ = B3.xyz 

458 return x, y, _zidw(x, y, useZ, B1, B2, B3) 

459 

460 

461def pierlotx(point1, point2, point3, alpha1, alpha2, alpha3, useZ=False, 

462 Clas=None, **Clas_kwds): 

463 '''3-Point resection using U{Pierlot<http://www.Telecom.ULg.ac.BE/publi/ 

464 publications/pierlot/Pierlot2014ANewThree>}'s method C{ToTal} with 

465 I{exact} limits for the (pseudo-)singularities. 

466 

467 @arg point1: First point (C{Cartesian}, L{Vector3d}, C{Vector3Tuple}, 

468 C{Vector4Tuple} or C{Vector2Tuple} if C{B{useZ}=False}). 

469 @arg point2: Second point (C{Cartesian}, L{Vector3d}, C{Vector3Tuple}, 

470 C{Vector4Tuple} or C{Vector2Tuple} if C{B{useZ}=False}). 

471 @arg point3: Third point (C{Cartesian}, L{Vector3d}, C{Vector3Tuple}, 

472 C{Vector4Tuple} or C{Vector2Tuple} if C{B{useZ}=False}). 

473 @arg alpha1: Angle at B{C{point1}} (C{degrees}, counter-clockwise). 

474 @arg alpha2: Angle at B{C{point2}} (C{degrees}, counter-clockwise). 

475 @arg alpha3: Angle at B{C{point3}} (C{degrees}, counter-clockwise). 

476 @kwarg useZ: If C{True}, interpolate the survey point's Z component, 

477 otherwise use C{z=INT0} (C{bool}). 

478 @kwarg Clas: Optional class to return the survey point, if C{None} use 

479 B{C{point1}}'s (sub-)class. 

480 @kwarg Clas_kwds: Optional, additional keyword arguments for the survey 

481 point instance. 

482 

483 @return: The survey (or robot) point, an instance of B{C{Clas}} or if 

484 C{B{Clas} is None} an instance of B{C{point1}}'s (sub-)class. 

485 

486 @raise ResectionError: Near-coincident, -colinear or -concyclic points or 

487 invalid B{C{alpha1}}, B{C{alpha2}} or B{C{alpha3}}. 

488 

489 @raise TypeError: Invalid B{C{point1}}, B{C{point2}} or B{C{point3}}. 

490 

491 @see: I{Pierlot's} C function U{triangulationPierlot2<http://www.Telecom.ULg.ac.BE/ 

492 triangulation/doc/total_8c_source.html>} and function L{pygeodesy.pierlot}. 

493 ''' 

494 

495 def _a_z_Bs(Bs, *alphas): 

496 a3 = map(_umod_360, alphas) # 0 <= alphas < 360 

497 (a1, a2, a3), (B1, B2, B3) = zip(*sorted(zip(a3, Bs))) 

498 for a, d, B in ((a1, a2, B1), (a2, a3, B2), (a3, a1, B3)): 

499 d -= a # a12 = a2 - a1, ... 

500 z = isnear0(fabs(d) % _180_0) 

501 yield d, z, B 

502 

503 def _cot(s, c): # I{exact} cotangent 

504 try: 

505 return (c / s) if c else (NEG0 if s < 0 else _0_0) 

506 except ZeroDivisionError: 

507 raise ValueError(_or(_coincident_, _colinear_)) 

508 

509 Bs = _B3(useZ, point1, point2, point3) 

510 try: 

511 C = [0] # pseudo-global, passing the exception Case 

512 x, y, z = _pierlotx3(_a_z_Bs(Bs, alpha1, alpha2, alpha3), 

513 useZ, _cot, C.append) 

514 clas = Clas or point1.classof 

515 return clas(x, y, z, **_xkwds(Clas_kwds, name=pierlotx.__name__)) 

516 

517 except (TypeError, ValueError) as x: 

518 raise ResectionError(point1=point1, point2=point2, point3=point3, C=C.pop(), 

519 alpha1=alpha1, alpha2=alpha2, alpha3=alpha3, cause=x) 

520 

521 

522def snellius3(a, b, degC, alpha, beta): 

523 '''Snellius' surveying using U{Snellius Pothenot<https://WikiPedia.org/wiki/Snellius–Pothenot_problem>}. 

524 

525 @arg a: Length of the triangle side between corners C{B} and C{C} and opposite of 

526 triangle corner C{A} (C{scalar}, non-negative C{meter}, conventionally). 

527 @arg b: Length of the triangle side between corners C{C} and C{A} and opposite of 

528 triangle corner C{B} (C{scalar}, non-negative C{meter}, conventionally). 

529 @arg degC: Angle at triangle corner C{C}, opposite triangle side C{c} (non-negative C{degrees}). 

530 @arg alpha: Angle subtended by triangle side B{C{b}} (non-negative C{degrees}). 

531 @arg beta: Angle subtended by triangle side B{C{a}} (non-negative C{degrees}). 

532 

533 @return: L{Survey3Tuple}C{(PA, PB, PC)} with distance from survey point C{P} to 

534 each of the triangle corners C{A}, C{B} and C{C}, same units as triangle 

535 sides B{C{a}}, B{C{b}} and B{C{c}}. 

536 

537 @raise TriangleError: Invalid B{C{a}}, B{C{b}} or B{C{degC}} or negative B{C{alpha}} 

538 or B{C{beta}}. 

539 

540 @see: Function L{pygeodesy.wildberger3}. 

541 ''' 

542 try: 

543 a, b, degC, alpha, beta = t = map1(float, a, b, degC, alpha, beta) 

544 if min(t) < 0: 

545 raise ValueError(_negative_) 

546 ra, rb, rC = map1(radians, alpha, beta, degC) 

547 

548 r = fsumf_(ra, rb, rC) * _0_5 

549 k = PI - r 

550 if min(k, r) < 0: 

551 raise ValueError(_or(_coincident_, _colinear_)) 

552 

553 sa, sb = sin(ra), sin(rb) 

554 p = atan2(a * sa, b * sb) 

555 sp, cp, sr, cr = sincos2_(PI_4 - p, r) 

556 w = atan2(sp * sr, cp * cr) 

557 x = k + w 

558 y = k - w 

559 

560 s = fabs(sa) 

561 if fabs(sb) > s: 

562 pc = fabs(a * sin(y) / sb) 

563 elif s: 

564 pc = fabs(b * sin(x) / sa) 

565 else: 

566 raise ValueError(_or(_colinear_, _coincident_)) 

567 

568 pa = _triSide(b, pc, fsumf_(PI, -ra, -x)) 

569 pb = _triSide(a, pc, fsumf_(PI, -rb, -y)) 

570 return Survey3Tuple(pa, pb, pc, name=snellius3.__name__) 

571 

572 except (TypeError, ValueError) as x: 

573 raise TriangleError(a=a, b=b, degC=degC, alpha=alpha, beta=beta, cause=x) 

574 

575 

576def tienstra7(pointA, pointB, pointC, alpha, beta=None, gamma=None, 

577 useZ=False, Clas=None, **Clas_kwds): 

578 '''3-Point resection using U{Tienstra<https://WikiPedia.org/wiki/Tienstra_formula>}'s formula. 

579 

580 @arg pointA: First point (C{Cartesian}, L{Vector3d}, C{Vector3Tuple}, C{Vector4Tuple} or 

581 C{Vector2Tuple} if C{B{useZ}=False}). 

582 @arg pointB: Second point (C{Cartesian}, L{Vector3d}, C{Vector3Tuple}, C{Vector4Tuple} or 

583 C{Vector2Tuple} if C{B{useZ}=False}). 

584 @arg pointC: Third point (C{Cartesian}, L{Vector3d}, C{Vector3Tuple}, C{Vector4Tuple} or 

585 C{Vector2Tuple} if C{B{useZ}=False}). 

586 @arg alpha: Angle subtended by triangle side C{a} from B{C{pointB}} to B{C{pointC}} 

587 (C{degrees}, non-negative). 

588 @kwarg beta: Angle subtended by triangle side C{b} from B{C{pointA}} to B{C{pointC}} 

589 (C{degrees}, non-negative) or C{None} if C{B{gamma} is not None}. 

590 @kwarg gamma: Angle subtended by triangle side C{c} from B{C{pointA}} to B{C{pointB}} 

591 (C{degrees}, non-negative) or C{None} if C{B{beta} is not None}. 

592 @kwarg useZ: If C{True}, use and interpolate the Z component, otherwise force C{z=INT0} 

593 (C{bool}). 

594 @kwarg Clas: Optional class to return the survey point or C{None} for B{C{pointA}}'s 

595 (sub-)class. 

596 @kwarg Clas_kwds: Optional, additional keyword arguments for the survey point instance. 

597 

598 @note: Points B{C{pointA}}, B{C{pointB}} and B{C{pointC}} are ordered clockwise. 

599 

600 @return: L{Tienstra7Tuple}C{(pointP, A, B, C, a, b, c)} with survey C{pointP}, an 

601 instance of B{C{Clas}} or if C{B{Clas} is None} an instance of B{C{pointA}}'s 

602 (sub-)class, with triangle angles C{A} at B{C{pointA}}, C{B} at B{C{pointB}} 

603 and C{C} at B{C{pointC}} in C{degrees} and with triangle sides C{a}, C{b} and 

604 C{c} in C{meter}, conventionally. 

605 

606 @raise ResectionError: Near-coincident, -colinear or -concyclic points or sum of 

607 B{C{alpha}}, B{C{beta}} and B{C{gamma}} not C{360} or negative 

608 B{C{alpha}}, B{C{beta}} or B{C{gamma}}. 

609 

610 @raise TypeError: Invalid B{C{pointA}}, B{C{pointB}} or B{C{pointC}}. 

611 

612 @see: U{3-Point Resection Solver<http://MesaMike.org/geocache/GC1B0Q9/tienstra/>}, 

613 U{V. Pierlot, M. Van Droogenbroeck, "A New Three Object Triangulation..." 

614 <http://www.Telecom.ULg.ac.BE/publi/publications/pierlot/Pierlot2014ANewThree/>}, 

615 U{18 Triangulation Algorithms...<http://www.Telecom.ULg.ac.BE/triangulation/>} and 

616 functions L{pygeodesy.cassini}, L{pygeodesy.collins5} and L{pygeodesy.pierlot}. 

617 ''' 

618 

619 def _deg_ks(r, s, ks, N): 

620 if isnear0(fsumf_(PI, r, -s)): # r + (PI2 - s) == PI 

621 raise ValueError(Fmt.PARENSPACED(concyclic=N)) 

622 # k = 1 / (cot(r) - cot(s)) 

623 # = 1 / (cos(r) / sin(r) - cos(s) / sin(s)) 

624 # = 1 / (cos(r) * sin(s) - cos(s) * sin(r)) / (sin(r) * sin(s)) 

625 # = sin(r) * sin(s) / (cos(r) * sin(s) - cos(s) * sin(r)) 

626 sr, cr, ss, cs = sincos2_(r, s) 

627 c = fsum1f_(cr * ss, -cs * sr) 

628 if isnear0(c): 

629 raise ValueError(Fmt.PARENSPACED(cotan=N)) 

630 ks.append(sr * ss / c) 

631 return Degrees(degrees(r), name=N) # C degrees 

632 

633 A, B, C = _ABC3(useZ, pointA, pointB, pointC) 

634 try: 

635 sa, sb, sc = map1(radians, alpha, (beta or 0), (gamma or 0)) 

636 if beta is None: 

637 if gamma is None: 

638 raise ValueError(_and(Fmt.EQUAL(beta=beta), Fmt.EQUAL(gamma=gamma))) 

639 sb = fsumf_(PI2, -sa, -sc) 

640 elif gamma is None: 

641 sc = fsumf_(PI2, -sa, -sb) 

642 else: # subtended angles must add to 360 degrees 

643 r = fsum1f_(sa, sb, sc) 

644 if fabs(r - PI2) > EPS: 

645 raise ValueError(Fmt.EQUAL(sum=degrees(r))) 

646 if min(sa, sb, sc) < 0: 

647 raise ValueError(_negative_) 

648 

649 # triangle sides 

650 a = B.minus(C).length 

651 b = A.minus(C).length 

652 c = A.minus(B).length 

653 

654 ks = [] # 3 Ks and triangle angles 

655 dA = _deg_ks(_triAngle(b, c, a), sa, ks, _A_) 

656 dB = _deg_ks(_triAngle(a, c, b), sb, ks, _B_) 

657 dC = _deg_ks(_triAngle(a, b, c), sc, ks, _C_) 

658 

659 k = fsum1(ks, floats=True) 

660 if isnear0(k): 

661 raise ValueError(Fmt.EQUAL(K=k)) 

662 x = Fdot(ks, A.x, B.x, C.x).fover(k) 

663 y = Fdot(ks, A.y, B.y, C.y).fover(k) 

664 z = _zidw(x, y, useZ, A, B, C) 

665 

666 n = tienstra7.__name__ 

667 clas = Clas or pointA.classof 

668 P = clas(x, y, z, **_xkwds(Clas_kwds, name=n)) 

669 return Tienstra7Tuple(P, dA, dB, dC, a, b, c, name=n) 

670 

671 except (TypeError, ValueError) as x: 

672 raise ResectionError(pointA=pointA, pointB=pointB, pointC=pointC, 

673 alpha=alpha, beta=beta, gamma=gamma, cause=x) 

674 

675 

676def triAngle(a, b, c): 

677 '''Compute one angle of a triangle. 

678 

679 @arg a: Adjacent triangle side length (C{scalar}, non-negative 

680 C{meter}, conventionally). 

681 @arg b: Adjacent triangle side length (C{scalar}, non-negative 

682 C{meter}, conventionally). 

683 @arg c: Opposite triangle side length (C{scalar}, non-negative 

684 C{meter}, conventionally). 

685 

686 @return: Angle in C{radians} at triangle corner C{C}, opposite 

687 triangle side B{C{c}}. 

688 

689 @raise TriangleError: Invalid or negative B{C{a}}, B{C{b}} or B{C{c}}. 

690 

691 @see: Functions L{pygeodesy.triAngle4} and L{pygeodesy.triSide}. 

692 ''' 

693 try: 

694 return _triAngle(a, b, c) 

695 except (TypeError, ValueError) as x: 

696 raise TriangleError(a=a, b=b, c=c, cause=x) 

697 

698 

699def _triAngle(a, b, c): 

700 # (INTERNAL) To allow callers to embellish errors 

701 a, b, c = map1(float, a, b, c) 

702 if a < b: 

703 a, b = b, a 

704 if b < 0 or c < 0: 

705 raise ValueError(_negative_) 

706 if a < EPS0: 

707 raise ValueError(_coincident_) 

708 b_a = b / a 

709 if b_a < EPS0: 

710 raise ValueError(_coincident_) 

711 t = fsumf_(_1_0, b_a**2, -(c / a)**2) / (b_a * _2_0) 

712 return acos1(t) 

713 

714 

715def triAngle4(a, b, c): 

716 '''Compute the angles of a triangle. 

717 

718 @arg a: Length of the triangle side opposite of triangle corner C{A} 

719 (C{scalar}, non-negative C{meter}, conventionally). 

720 @arg b: Length of the triangle side opposite of triangle corner C{B} 

721 (C{scalar}, non-negative C{meter}, conventionally). 

722 @arg c: Length of the triangle side opposite of triangle corner C{C} 

723 (C{scalar}, non-negative C{meter}, conventionally). 

724 

725 @return: L{TriAngle4Tuple}C{(radA, radB, radC, rIn)} with angles C{radA}, 

726 C{radB} and C{radC} at triangle corners C{A}, C{B} and C{C}, all 

727 in C{radians} and the C{InCircle} radius C{rIn} aka C{inradius}, 

728 same units as triangle sides B{C{a}}, B{C{b}} and B{C{c}}. 

729 

730 @raise TriangleError: Invalid or negative B{C{a}}, B{C{b}} or B{C{c}}. 

731 

732 @see: Function L{pygeodesy.triAngle}. 

733 ''' 

734 try: 

735 a, b, c = map1(float, a, b, c) 

736 ab = a < b 

737 if ab: 

738 a, b = b, a 

739 bc = b < c 

740 if bc: 

741 b, c = c, b 

742 

743 if c > EPS0: # c = min(a, b, c) 

744 s = fsumf_(a, b, c) * _0_5 

745 if s < EPS0: 

746 raise ValueError(_negative_) 

747 sa, sb, sc = (s - a), (s - b), (s - c) 

748 r = sa * sb * sc / s 

749 if r < EPS02: 

750 raise ValueError(_coincident_) 

751 r = sqrt(r) 

752 rA = atan2(r, sa) * _2_0 

753 rB = atan2(r, sb) * _2_0 

754 rC = fsumf_(PI, -rA, -rB) 

755 if min(rA, rB, rC) < 0: 

756 raise ValueError(_colinear_) 

757 elif c < 0: 

758 raise ValueError(_negative_) 

759 else: # 0 <= c <= EPS0 

760 rA = rB = PI_2 

761 rC = r = _0_0 

762 

763 if bc: 

764 rB, rC = rC, rB 

765 if ab: 

766 rA, rB = rB, rA 

767 return TriAngle4Tuple(rA, rB, rC, r, name=triAngle4.__name__) 

768 

769 except (TypeError, ValueError) as x: 

770 raise TriangleError(a=a, b=b, c=c, cause=x) 

771 

772 

773def triSide(a, b, radC): 

774 '''Compute one side of a triangle. 

775 

776 @arg a: Adjacent triangle side length (C{scalar}, 

777 non-negative C{meter}, conventionally). 

778 @arg b: Adjacent triangle side length (C{scalar}, 

779 non-negative C{meter}, conventionally). 

780 @arg radC: Angle included by sides B{C{a}} and B{C{b}}, 

781 opposite triangle side C{c} (C{radians}). 

782 

783 @return: Length of triangle side C{c}, opposite triangle 

784 corner C{C} and angle B{C{radC}}, same units as 

785 B{C{a}} and B{C{b}}. 

786 

787 @raise TriangleError: Invalid B{C{a}}, B{C{b}} or B{C{radC}}. 

788 

789 @see: Functions L{pygeodesy.sqrt_a}, L{pygeodesy.triAngle}, 

790 L{pygeodesy.triSide2} and L{pygeodesy.triSide4}. 

791 ''' 

792 try: 

793 return _triSide(a, b, radC) 

794 except (TypeError, ValueError) as x: 

795 raise TriangleError(a=a, b=b, radC=radC, cause=x) 

796 

797 

798def _triSide(a, b, radC): 

799 # (INTERNAL) To allow callers to embellish errors 

800 a, b, r = t = map1(float, a, b, radC) 

801 if min(t) < 0: 

802 raise ValueError(_negative_) 

803 

804 if a < b: 

805 a, b = b, a 

806 if a > EPS0: 

807 ba = b / a 

808 c2 = fsumf_(_1_0, ba**2, _N_2_0 * ba * cos(r)) 

809 if c2 > EPS02: 

810 return a * sqrt(c2) 

811 elif c2 < 0: 

812 raise ValueError(_invalid_) 

813 return hypot(a, b) 

814 

815 

816def triSide2(b, c, radB): 

817 '''Compute one side and the opposite angle of a triangle. 

818 

819 @arg b: Adjacent triangle side length (C{scalar}, 

820 non-negative C{meter}, conventionally). 

821 @arg c: Adjacent triangle side length (C{scalar}, 

822 non-negative C{meter}, conventionally). 

823 @arg radB: Angle included by sides B{C{a}} and B{C{c}}, 

824 opposite triangle side C{b} (C{radians}). 

825 

826 @return: L{TriSide2Tuple}C{(a, radA)} with triangle angle 

827 C{radA} in C{radians} and length of the opposite 

828 triangle side C{a}, same units as B{C{b}} and B{C{c}}. 

829 

830 @raise TriangleError: Invalid B{C{b}} or B{C{c}} or either 

831 B{C{b}} or B{C{radB}} near zero. 

832 

833 @see: Functions L{pygeodesy.sqrt_a}, L{pygeodesy.triSide} 

834 and L{pygeodesy.triSide4}. 

835 ''' 

836 try: 

837 return _triSide2(b, c, radB) 

838 except (TypeError, ValueError) as x: 

839 raise TriangleError(b=b, c=c, radB=radB, cause=x) 

840 

841 

842def _triSide2(b, c, radB): 

843 # (INTERNAL) To allow callers to embellish errors 

844 b, c, rB = map1(float, b, c, radB) 

845 if min(b, c, rB) < 0: 

846 raise ValueError(_negative_) 

847 sB, cB = sincos2(rB) 

848 if isnear0(sB): 

849 if not isnear0(b): 

850 raise ValueError(_invalid_) 

851 if cB < 0: 

852 a, rA = (b + c), PI 

853 else: 

854 a, rA = fabs(b - c), _0_0 

855 elif isnear0(b): 

856 raise ValueError(_invalid_) 

857 else: 

858 rA = fsumf_(PI, -rB, -asin1(c * sB / b)) 

859 a = sin(rA) * b / sB 

860 return TriSide2Tuple(a, rA, name=triSide2.__name__) 

861 

862 

863def triSide4(radA, radB, c): 

864 '''Compute two sides and the height of a triangle. 

865 

866 @arg radA: Angle at triangle corner C{A}, opposite triangle side C{a} 

867 (non-negative C{radians}). 

868 @arg radB: Angle at triangle corner C{B}, opposite triangle side C{b} 

869 (non-negative C{radians}). 

870 @arg c: Length of triangle side between triangle corners C{A} and C{B}, 

871 (C{scalar}, non-negative C{meter}, conventionally). 

872 

873 @return: L{TriSide4Tuple}C{(a, b, radC, d)} with triangle sides C{a} and 

874 C{b} and triangle height C{d} perpendicular to triangle side 

875 B{C{c}}, all in the same units as B{C{c}} and interior angle 

876 C{radC} in C{radians} at triangle corner C{C}, opposite 

877 triangle side B{C{c}}. 

878 

879 @raise TriangleError: Invalid or negative B{C{radA}}, B{C{radB}} or B{C{c}}. 

880 

881 @see: U{Triangulation, Surveying<https://WikiPedia.org/wiki/Triangulation_(surveying)>} 

882 and functions L{pygeodesy.sqrt_a}, L{pygeodesy.triSide} and L{pygeodesy.triSide2}. 

883 ''' 

884 try: 

885 rA, rB, c = map1(float, radA, radB, c) 

886 rC = fsumf_(PI, -rA, -rB) 

887 if min(rC, rA, rB, c) < 0: 

888 raise ValueError(_negative_) 

889 sa, ca, sb, cb = sincos2_(rA, rB) 

890 sc = fsum1f_(sa * cb, sb * ca) 

891 if sc < EPS0 or min(sa, sb) < 0: 

892 raise ValueError(_invalid_) 

893 sc = c / sc 

894 return TriSide4Tuple((sa * sc), (sb * sc), rC, (sa * sb * sc), 

895 name=triSide4.__name__) 

896 

897 except (TypeError, ValueError) as x: 

898 raise TriangleError(radA=radA, radB=radB, c=c, cause=x) 

899 

900 

901def wildberger3(a, b, c, alpha, beta, R3=min): 

902 '''Snellius' surveying using U{Rational Trigonometry 

903 <https://WikiPedia.org/wiki/Snellius–Pothenot_problem>}. 

904 

905 @arg a: Length of the triangle side between corners C{B} and C{C} and opposite of 

906 triangle corner C{A} (C{scalar}, non-negative C{meter}, conventionally). 

907 @arg b: Length of the triangle side between corners C{C} and C{A} and opposite of 

908 triangle corner C{B} (C{scalar}, non-negative C{meter}, conventionally). 

909 @arg c: Length of the triangle side between corners C{A} and C{B} and opposite of 

910 triangle corner C{C} (C{scalar}, non-negative C{meter}, conventionally). 

911 @arg alpha: Angle subtended by triangle side B{C{b}} (C{degrees}, non-negative). 

912 @arg beta: Angle subtended by triangle side B{C{a}} (C{degrees}, non-negative). 

913 @kwarg R3: Callable to determine C{R3} from C{(R3 - C)**2 = D}, typically standard 

914 Python function C{min} or C{max}, invoked with 2 arguments. 

915 

916 @return: L{Survey3Tuple}C{(PA, PB, PC)} with distance from survey point C{P} to 

917 each of the triangle corners C{A}, C{B} and C{C}, same units as B{C{a}}, 

918 B{C{b}} and B{C{c}}. 

919 

920 @raise TriangleError: Invalid B{C{a}}, B{C{b}} or B{C{c}} or negative B{C{alpha}} or 

921 B{C{beta}} or B{C{R3}} not C{callable}. 

922 

923 @see: U{Wildberger, Norman J.<https://Math.Sc.Chula.ac.TH/cjm/content/ 

924 survey-article-greek-geometry-rational-trigonometry-and-snellius-–-pothenot-surveying>}, 

925 U{Devine Proportions, page 252<http://www.MS.LT/derlius/WildbergerDivineProportions.pdf>} 

926 and function L{pygeodesy.snellius3}. 

927 ''' 

928 def _s(x): 

929 return sin(x)**2 

930 

931 def _vpa(r1, r3, q2, q3, s3): 

932 r = r1 * r3 * _4_0 

933 n = (r - Fsum(r1, r3, -q2).fpow(2)).fover(s3) 

934 if n < 0 or isnear0(r): 

935 raise ValueError(_coincident_) 

936 return sqrt((n / r) * q3) if n else _0_0 

937 

938 try: 

939 a, b, c, da, db = t = map1(float, a, b, c, alpha, beta) 

940 if min(t) < 0: 

941 raise ValueError(_negative_) 

942 

943 ra, rb = radians(da), radians(db) 

944 s1, s2, s3 = s = map1(_s, rb, ra, ra + rb) # rb, ra! 

945 if min(s) < EPS02: 

946 raise ValueError(_or(_coincident_, _colinear_)) 

947 

948 q1, q2, q3 = q = a**2, b**2, c**2 

949 if min(q) < EPS02: 

950 raise ValueError(_coincident_) 

951 

952 r1 = s2 * q3 / s3 # s2! 

953 r2 = s1 * q3 / s3 # s1! 

954 Qs = Fsum(*q) # == hypot2_(a, b, c) 

955 Ss = Fsum(*s) # == fsum1(s, floats=True) 

956 s += Qs * _0_5, # tuple! 

957 C0 = Fdot(s, q1, q2, q3, -Ss) 

958 r3 = C0.fover(-s3) 

959 d0 = Qs.fpow(2).fsub_(hypot2_(*q) * _2_0).fmul(s1 * s2).fover(s3) 

960 if d0 > EPS02: # > c0 

961 d0 = sqrt(d0) 

962 if not callable(R3): 

963 raise ValueError(_R3__ + _not_(callable.__name__)) 

964 r3 = R3(float(C0 + d0), float(C0 - d0)) # XXX min or max 

965 elif d0 < 0: 

966 raise ValueError(_negative_) 

967 

968 pa = _vpa(r1, r3, q2, q3, s3) 

969 pb = _vpa(r2, r3, q1, q3, s3) 

970 pc = favg(_triSide2(b, pa, ra).a, 

971 _triSide2(a, pb, rb).a) 

972 return Survey3Tuple(pa, pb, pc, name=wildberger3.__name__) 

973 

974 except (TypeError, ValueError) as x: 

975 raise TriangleError(a=a, b=b, c=c, alpha=alpha, beta=beta, R3=R3, cause=x) 

976 

977 

978def _zidw(x, y, useZ, *ABC): 

979 if useZ: # interpolate z or coplanar with A, B and C? 

980 t = tuple(_.z for _ in ABC) 

981 v = Vector3d(x, y, fmean(t)) 

982 z = fidw(t, (v.minus(T).length for T in ABC)) 

983 else: 

984 z = INT0 

985 return z 

986 

987# **) MIT License 

988# 

989# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved. 

990# 

991# Permission is hereby granted, free of charge, to any person obtaining a 

992# copy of this software and associated documentation files (the "Software"), 

993# to deal in the Software without restriction, including without limitation 

994# the rights to use, copy, modify, merge, publish, distribute, sublicense, 

995# and/or sell copies of the Software, and to permit persons to whom the 

996# Software is furnished to do so, subject to the following conditions: 

997# 

998# The above copyright notice and this permission notice shall be included 

999# in all copies or substantial portions of the Software. 

1000# 

1001# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS 

1002# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 

1003# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 

1004# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR 

1005# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, 

1006# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR 

1007# OTHER DEALINGS IN THE SOFTWARE.