1

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

3

4 u'''3-Point resection functions L{cassini}, L{collins5}, L{pierlot}, L{pierlotx} and

5 L{tienstra7}, survey functions L{snellius3} and L{wildberger3} and triangle functions

6 L{triAngle}, L{triAngle5}, 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

13 from __future__ import division as _; del _ # PYCHOK semicolon

14

15 from pygeodesy.basics import map1, map2, _zip, _ALL_LAZY

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

17 _0_0, _0_5, _1_0, _N_1_0, _2_0, _N_2_0, _4_0, _16_0, \

18 _180_0, _360_0, isnear0, _over, _umod_360

19 from pygeodesy.errors import _and, _or, TriangleError, _ValueError, _xcallable, \

20 _xkwds, _xkwds_pop2

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

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

23 from pygeodesy.interns import _a_, _A_, _area_, _b_, _B_, _c_, _C_, _coincident_, \

24 _colinear_, _d_, _eps_, _invalid_, _negative_, _not_, \

25 _rIn_, _SPACE_

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

27 from pygeodesy.named import _NamedTuple, _Pass, Fmt

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

29 from pygeodesy.units import Degrees, Distance, Radians

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

31 from pygeodesy.vector3d import _otherV3d, Vector3d

32

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

34

35__all__ = _ALL_LAZY.resections

36__version__ = '24.03.24'

37

38_concyclic_ = 'concyclic'

39_PA_ = 'PA'

40_PB_ = 'PB'

41_PC_ = 'PC'

42_pointH_ = 'pointH'

43_pointP_ = 'pointP'

44_positive_ = 'positive'

45_radA_ = 'radA'

46_radB_ = 'radB'

47_radC_ = 'radC'

48

49

50 class Collins5Tuple(_NamedTuple):

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

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

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

54 '''

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

56 _Units_ = (_Pass, _Pass, Distance, Distance, Distance)

57

58

59 def _F1(*xs): # class

60 '''(INTERNAL) An L{Fsum}, 1-primed.

61 '''

62 F = Fsum(_1_0, *xs)

63 F += _N_1_0

64 return F

65

66

67 class ResectionError(_ValueError):

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

69 '''

70 pass

71

72

73 class Survey3Tuple(_NamedTuple):

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

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

76 '''

77 _Names_ = (_PA_, _PB_, _PC_)

78 _Units_ = ( Distance, Distance, Distance)

79

80

81 class Tienstra7Tuple(_NamedTuple):

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

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

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

85 '''

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

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

88

89

90 class TriAngle5Tuple(_NamedTuple):

91 '''5-Tuple C{(radA, radB, radC, rIn, area)} with the interior angles at

92 triangle corners C{A}, C{B} and C{C} in C{radians}, the C{InCircle}

93 radius C{rIn} aka C{inradius} in C{meter} and the triangle C{area}

94 in C{meter} I{squared}, conventionally.

95 '''

96 _Names_ = (_radA_, _radB_, _radC_, _rIn_, _area_)

97 _Units_ = ( Radians, Radians, Radians, Distance, _Pass)

98

99

100 class TriSide2Tuple(_NamedTuple):

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

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

103 '''

104 _Names_ = (_a_, _radA_)

105 _Units_ = ( Distance, Radians)

106

107

108 class TriSide4Tuple(_NamedTuple):

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

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

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

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

113 '''

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

115 _Units_ = ( Distance, Distance, Radians, Distance)

116

117

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

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

120 '''

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

122 _otherV3d(useZ=useZ, pointB=pointB),

123 _otherV3d(useZ=useZ, pointC=pointC))

124

125

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

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

128 '''

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

130 _otherV3d(useZ=useZ, point2=point2),

131 _otherV3d(useZ=useZ, point3=point3))

132

133

134 def cassini(pointA, pointB, pointC, alpha, beta, useZ=False, **Clas_and_kwds):

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

136

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

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

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

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

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

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

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

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

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

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

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

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

149 @kwarg Clas_and_kwds: Optional class C{B{Clas}=B{pointA}.classof} to

150 return the survey point with optionally other B{C{Clas}}

151 keyword arguments to instantiate the survey point.

152

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

154

155 @return: The survey point, an instance of B{C{Clas}} or B{C{pointA}}'s

156 (sub-)class.

157

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

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

160

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

162

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

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

165 and functions L{collins5}, L{pierlot}, L{pierlotx} and L{tienstra7}.

166 '''

167

168 def _H(A, C, sa):

169 s, c = sincos2d(sa)

170 if isnear0(s):

171 raise ValueError(_or(_coincident_, _colinear_))

172 t = s, c, c

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

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

175 return x, y

176

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

178 try:

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

180 if min(sa, sb) < 0:

181 raise ValueError(_negative_)

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

183 raise ValueError()

184

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

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

187

188 x = x1 - x2

189 y = y1 - y2

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

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

192

193 m = y / x

194 n = x / y

195 N = n + m

196 if isnear0(N):

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

198

199 t = n, m, _1_0, _N_1_0

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

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

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

203 return _Clas(cassini, pointA, Clas_and_kwds, x, y, z)

204

205 except (TypeError, ValueError) as x:

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

207 alpha=alpha, beta=beta, cause=x)

208

209

210 def _Clas(where, point, Clas_and_kwds, *args):

211 '''(INTERNAL) Return a C{B{Clas}=point.classof} survey point.

212 '''

213 Clas, kwds = _xkwds_pop2(Clas_and_kwds, Clas=point.classof)

214 return Clas(*args, **_xkwds(kwds, name=where.__name__))

215

216

217 def collins5(pointA, pointB, pointC, alpha, beta, useZ=False, **Clas_and_kwds):

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

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

220

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

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

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

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

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

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

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

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

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

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

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

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

233 @kwarg Clas_and_kwds: Optional class C{B{Clas}=B{pointA}.classof} to

234 return the survey point with optionally other B{C{Clas}}

235 keyword arguments to instantiate the survey point.

236

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

238

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

240 auxiliary C{pointH}, each an instance of B{C{Clas}} or B{C{pointA}}'s

241 (sub-)class and triangle sides C{a}, C{b} and C{c} in C{meter},

242 conventionally.

243

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

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

246

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

248

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

250 and functions L{cassini}, L{pierlot}, L{pierlotx} and L{tienstra7}.

251 '''

252

253 def _azi_len2(A, B, pi2=PI2):

254 v = B.minus(A)

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

256 if r < 0 and pi2:

257 r += pi2

258 return r, v.length

259

260 def _xyz(d, r, A, B, C, useZ):

261 s, c = sincos2(r)

262 x = A.x + d * s

263 y = A.y + d * c

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

265 return x, y, z

266

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

268 try:

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

270 if min(ra, rb) < 0:

271 raise ValueError(_negative_)

272

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

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

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

276

277# za, a = _azi_len2(C, B)

278 zb, b = _azi_len2(C, A)

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

280

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

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

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

284 H = _xyz(d, r, A, B, C, useZ)

285

286 zh, _ = _azi_len2(C, Vector3d(*H))

287

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

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

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

291 P = _xyz(d, r, A, B, C, useZ)

292 P = _Clas(collins5, pointA, Clas_and_kwds, *P)

293

294 H = _Clas(collins5, pointA, Clas_and_kwds, *H)

295 a = B.minus(C).length

296

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

298

299 except (TypeError, ValueError) as x:

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

301 alpha=alpha, beta=beta, cause=x)

302

303

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

305 **Clas_and_kwds):

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

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

308 the (pseudo-)singularities.

309

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

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

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

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

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

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

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

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

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

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

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

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

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

323 @kwarg Clas_and_kwds: Optional class C{B{Clas}=B{point1}.classof} to

324 return the survey point with optionally other B{C{Clas}}

325 keyword arguments to instantiate the survey point.

326

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

328 by angle, modulo 360, counter-clockwise.

329

330 @return: The survey (or robot) point, an instance of B{C{Clas}} or B{C{point1}}'s

331 (sub-)class.

332

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

334 or invalid B{C{alpha12}} or B{C{alpha23}} or

335 non-positive B{C{eps}}.

336

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

338

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

340 triangulation/doc/total_8c_source.html>}, U{V. Pierlot, M. Van Droogenbroeck,

341 "A New Three Object Triangulation Algorithm for Mobile Robot Positioning"

342 <https://ORBi.ULiege.BE/bitstream/2268/157469/1/Pierlot2014ANewThree.pdf>},

343 U{Vincent Pierlot, Marc Van Droogenbroeck, "18 Triangulation Algorithms for 2D

344 Positioning (also known as the Resection Problem)"<http://www.Telecom.ULg.ac.BE/

345 triangulation>} and functions L{pierlotx}, L{cassini}, L{collins5} and L{tienstra7}.

346 '''

347

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

349 if eps > 0:

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

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

352

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

354 try:

355 xyz = _pierlot3(B1, B2, B3, alpha12, alpha23, useZ, _cot)

356 return _Clas(pierlot, point1, Clas_and_kwds, *xyz)

357

358 except (TypeError, ValueError) as x:

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

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

361

362

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

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

365 '''

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

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

368

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

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

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

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

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

374 # = c31 / s31

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

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

377

378 K = _F1(x3_ * x1_, cot31 * (y3_ * x1_),

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

380 if K:

381 cot12 = cot(s12, c12)

382 cot23 = cot(s23, c23)

383

384 # x12 = x1_ + cot12 * y1_

385 # y12 = y1_ - cot12 * x1_

386

387 # x23 = x3_ - cot23 * y3_

388 # y23 = y3_ + cot23 * x3_

389

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

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

392

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

394 X12_23 = _F1(x1_, cot12 * y1_, -x3_, cot23 * y3_)

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

396 Y12_23 = _F1(y1_, -cot12 * x1_, -y3_, -cot23 * x3_)

397

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

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

400 X31_23 = _F1(x1_, -cot31 * y1_, cot31 * y3_, cot23 * y3_)

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

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

403 Y31_23 = _F1(y1_, cot31 * x1_, -cot31 * x3_, -cot23 * x3_)

404

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

406 # = (x31 - x23) * (y12 - y23) - (x12 - x23) * (y31 - y23)

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

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

409 x, y = _pierlotxy2(B2, -K, Y12_23, X12_23, (X31_23 * Y12_23 -

410 X12_23 * Y31_23))

411 else:

412 x, y, _ = B2.xyz

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

414

415

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

417 **Clas_and_kwds):

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

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

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

421

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

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

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

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

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

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

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

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

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

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

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

433 @kwarg Clas_and_kwds: Optional class C{B{Clas}=B{point1}.classof} to

434 return the survey point with optionally other B{C{Clas}}

435 keyword arguments to instantiate the survey point.

436

437 @return: The survey (or robot) point, an instance of B{C{Clas}} or B{C{point1}}'s

438 (sub-)class.

439

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

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

442

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

444

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

446 triangulation/doc/total_8c_source.html>} and function L{pierlot}, L{cassini},

447 L{collins5} and L{tienstra7}.

448 '''

449

450 def _a_z_Bs(Bs, *alphas):

451 ds = map2(_umod_360, alphas) # 0 <= alphas < 360

452 ds, Bs = zip(*sorted(_zip(ds, Bs))) # unzip

453 for p, d, B in _zip(ds, _rotate(ds), Bs):

454 d -= p # a12 = a2 - a1, ...

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

456 yield d, z, B

457

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

459 try:

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

461 except ZeroDivisionError:

462 raise ValueError(_or(_coincident_, _colinear_))

463

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

465 try:

466 Cs = [0] # pseudo-global, passing the exception Case

467 xyz = _pierlotx3(_a_z_Bs(Bs, alpha1, alpha2, alpha3),

468 useZ, _cot, Cs.append)

469 return _Clas(pierlotx, point1, Clas_and_kwds, *xyz)

470

471 except (TypeError, ValueError) as x:

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

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

474

475

476 def _pierlotx3(a_z_Bs, useZ, cot, Cs):

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

478 '''

479 (a12, z12, B1), \

480 (a23, z23, B2), \

481 (a31, z31, B3) = a_z_Bs

482 if z12 and not z23:

483 Cs(1)

484 elif z23 and not z31:

485 Cs(2)

486 a23, B1, B2, B3 = a31, B2, B3, B1

487 elif z31 and not z12:

488 Cs(3)

489 a23, B2, B3 = a12, B3, B2

490 else:

491 Cs(4)

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

493

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

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

496

497 K = _F1(y1_ * x2_, -x1_ * y2_)

498 if K:

499 cot23 = cot(*sincos2d(a23))

500

501 # x23 = x2_ + cot23 * y2_

502 # y23 = y2_ - cot23 * x2_

503

504 # x31 = x1_ + cot23 * y1_

505 # y31 = y1_ - cot23 * x1_

506

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

508 X31_23 = _F1(x1_, cot23 * y1_, -x2_, -cot23 * y2_)

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

510 Y31_23 = _F1(y1_, -cot23 * x1_, -y2_, cot23 * x2_)

511

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

513 # x = (D * B3.x - K * Y31_23).fover(d)

514 # y = (D * B3.y + K * X31_23).fover(d)

515 x, y = _pierlotxy2(B3, K, Y31_23, X31_23, (X31_23 * _F1(x2_, -x1_) +

516 Y31_23 * _F1(y2_, -y1_)))

517 else:

518 x, y, _ = B3.xyz

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

520

521

522 def _pierlotxy2(B, K, X, Y, D):

523 '''(INTERNAL) Helper for C{_pierlot3} and C{_pierlotx3}.

524 '''

525 d = float(D)

526 if isnear0(d):

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

528 x = (D * B.x - K * X).fover(d)

529 y = (D * B.y + K * Y).fover(d)

530 return x, y

531

532

533 def _rotate(xs, n=1):

534 '''Rotate list or tuple C{xs} by C{n} items, right if C{n > 0} else left.

535 '''

536 return xs[n:] + xs[:n]

537

538

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

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

541

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

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

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

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

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

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

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

549

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

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

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

553

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

555 or B{C{beta}}.

556

557 @see: Function L{wildberger3}.

558 '''

559 try:

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

561 if min(t) < 0:

562 raise ValueError(_negative_)

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

564

565 r = fsum1f_(ra, rb, rC) * _0_5

566 k = PI - r

567 if min(k, r) < 0:

568 raise ValueError(_or(_coincident_, _colinear_))

569

570 sa, sb = map1(sin, ra, rb)

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

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

573 p = atan2(sp * sr, cp * cr)

574 pa = k + p

575 pb = k - p

576

577 if fabs(sb) > fabs(sa):

578 pc = fabs(a * sin(pb) / sb)

579 elif sa:

580 pc = fabs(b * sin(pa) / sa)

581 else:

582 raise ValueError(_or(_colinear_, _coincident_))

583

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

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

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

587

588 except (TypeError, ValueError) as x:

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

590

591

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

593 useZ=False, **Clas_and_kwds):

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

595

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

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

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

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

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

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

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

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

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

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

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

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

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

609 (C{bool}).

610 @kwarg Clas_and_kwds: Optional class C{B{Clas}=B{pointA}.classof} to return the survey

611 point with optionally other B{C{Clas}} keyword arguments to instantiate

612 the survey point.

613

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

615

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

617 instance of B{C{Clas}} or B{C{pointA}}'s (sub-)class, with triangle angles C{A}

618 at B{C{pointA}}, C{B} at B{C{pointB}} and C{C} at B{C{pointC}} in C{degrees}

619 and with triangle sides C{a}, C{b} and C{c} in C{meter}, conventionally.

620

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

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

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

624

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

626

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

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

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

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

631 functions L{cassini}, L{collins5}, L{pierlot} and L{pierlotx}.

632 '''

633

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

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

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

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

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

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

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

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

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

643 if isnear0(c):

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

645 ks.append(sr * ss / c)

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

647

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

649 try:

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

651 if beta is None:

652 if gamma is None:

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

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

655 elif gamma is None:

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

657 else: # subtended angles must add to 360 degrees

658 r = fsum1f_(sa, sb, sc)

659 if fabs(r - PI2) > EPS:

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

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

662 raise ValueError(_negative_)

663

664 # triangle sides

665 a = B.minus(C).length

666 b = A.minus(C).length

667 c = A.minus(B).length

668

669 ks = [] # 3 Ks and triangle angles

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

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

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

673

674 k = fsum1(ks, floats=True)

675 if isnear0(k):

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

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

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

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

680

681 P = _Clas(tienstra7, pointA, Clas_and_kwds, x, y, z)

682 return Tienstra7Tuple(P, dA, dB, dC, a, b, c, name=tienstra7.__name__)

683

684 except (TypeError, ValueError) as x:

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

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

687

688

689 def triAngle(a, b, c):

690 '''Compute one angle of a triangle.

691

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

693 C{meter}, conventionally).

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

695 C{meter}, conventionally).

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

697 C{meter}, conventionally).

698

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

700 triangle side B{C{c}}.

701

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

703

704 @see: Functions L{triAngle5} and L{triSide}.

705 '''

706 try:

707 return _triAngle(a, b, c)

708 except (TypeError, ValueError) as x:

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

710

711

712 def _triAngle(a, b, c):

713 # (INTERNAL) To allow callers to embellish errors

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

715 if a < b:

716 a, b = b, a

717 if b < 0 or c < 0:

718 raise ValueError(_negative_)

719 if a < EPS0:

720 raise ValueError(_coincident_)

721 b_a = b / a

722 if b_a < EPS0:

723 raise ValueError(_coincident_)

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

725 return acos1(t)

726

727

728 def triAngle5(a, b, c):

729 '''Compute the angles of a triangle.

730

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

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

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

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

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

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

737

738 @return: L{TriAngle5Tuple}C{(radA, radB, radC, rIn, area)} with angles

739 C{radA}, C{radB} and C{radC} at triangle corners C{A}, C{B}

740 and C{C}, all in C{radians}, the C{InCircle} radius C{rIn}

741 aka C{inradius}, same units as triangle sides B{C{a}},

742 B{C{b}} and B{C{c}} and the triangle C{area} in those same

743 units I{squared}.

744

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

746

747 @see: Functions L{triAngle} and L{triArea}.

748 '''

749 try:

750 x, y, z = map1(float, a, b, c)

751 ab = x < y

752 if ab:

753 x, y = y, x

754 bc = y < z

755 if bc:

756 y, z = z, y

757

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

759 s = fsum1f_(z, y, x) * _0_5

760 sa, sb, r = (s - x), (s - y), (s - z)

761 r *= _over(sa * sb, s)

762 if r < EPS02:

763 raise ValueError(_coincident_)

764 r = sqrt(r)

765 rA = atan2(r, sa) * _2_0

766 rB = atan2(r, sb) * _2_0

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

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

769 raise ValueError(_colinear_)

770 s *= r # Heron's area

771 elif z < 0:

772 raise ValueError(_negative_)

773 else: # 0 <= c <= EPS0

774 rA = rB = PI_2

775 rC = r = s = _0_0

776

777 if bc:

778 rB, rC = rC, rB

779 if ab:

780 rA, rB = rB, rA

781 return TriAngle5Tuple(rA, rB, rC, r, s, name=triAngle5.__name__)

782

783 except (TypeError, ValueError) as x:

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

785

786

787 def triArea(a, b, c):

788 '''Compute the area of a triangle using U{Heron's<https://

789 WikiPedia.org/wiki/Heron%27s_formula>} C{stable} formula.

790

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

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

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

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

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

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

797

798 @return: The triangle area (C{float}, conventionally C{meter} or

799 same units as B{C{a}}, B{C{b}} and B{C{c}} I{squared}).

800

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

802 '''

803 try:

804 r, y, x = sorted(map1(float, a, b, c))

805 if r > 0: # r = min(a, b, c)

806 ab = x - y

807 bc = y - r

808 y += r

809 r = (x + y) * (r - ab) * (r + ab) * (x + bc)

810 if r:

811 r = sqrt(r / _16_0)

812 elif r < 0:

813 raise ValueError(_negative_)

814 return r

815

816 except (TypeError, ValueError) as x:

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

818

819

820 def triSide(a, b, radC):

821 '''Compute one side of a triangle.

822

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

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

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

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

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

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

829

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

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

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

833

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

835

836 @see: Functions L{sqrt_a}, L{triAngle}, L{triSide2} and L{triSide4}.

837 '''

838 try:

839 return _triSide(a, b, radC)

840 except (TypeError, ValueError) as x:

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

842

843

844 def _triSide(a, b, radC):

845 # (INTERNAL) To allow callers to embellish errors

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

847 if min(t) < 0:

848 raise ValueError(_negative_)

849

850 if a < b:

851 a, b = b, a

852 if a > EPS0:

853 ba = b / a

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

855 if c2 > EPS02:

856 return a * sqrt(c2)

857 elif c2 < 0:

858 raise ValueError(_invalid_)

859 return hypot(a, b)

860

861

862 def triSide2(b, c, radB):

863 '''Compute a side and its opposite angle of a triangle.

864

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

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

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

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

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

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

871

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

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

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

875

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

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

878

879 @see: Functions L{sqrt_a}, L{triSide} and L{triSide4}.

880 '''

881 try:

882 return _triSide2(b, c, radB)

883 except (TypeError, ValueError) as x:

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

885

886

887 def _triSide2(b, c, radB):

888 # (INTERNAL) To allow callers to embellish errors

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

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

891 raise ValueError(_negative_)

892 sB, cB = sincos2(rB)

893 if isnear0(sB):

894 if not isnear0(b):

895 raise ValueError(_invalid_)

896 a, rA = ((b + c), PI) if cB < 0 else (fabs(b - c), _0_0)

897 elif isnear0(b):

898 raise ValueError(_invalid_)

899 else:

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

901 a = sin(rA) * b / sB

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

903

904

905 def triSide4(radA, radB, c):

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

907

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

909 (non-negative C{radians}).

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

911 (non-negative C{radians}).

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

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

914

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

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

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

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

919 triangle side B{C{c}}.

920

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

922

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

924 and functions L{sqrt_a}, L{triSide} and L{triSide2}.

925 '''

926 try:

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

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

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

930 raise ValueError(_negative_)

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

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

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

934 raise ValueError(_invalid_)

935 sc = c / sc

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

937 name=triSide4.__name__)

938

939 except (TypeError, ValueError) as x:

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

941

942

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

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

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

946

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

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

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

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

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

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

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

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

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

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

957

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

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

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

961

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

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

964

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

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

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

968 and function L{snellius3}.

969 '''

970 def _s(x):

971 return sin(x)**2

972

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

974 r = r1 * r3 * _4_0

975 n = (r - _F1(r1, r3, -q2).fpow(2)).fover(s3)

976 if n < 0 or isnear0(r):

977 raise ValueError(_coincident_)

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

979

980 try:

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

982 if min(t) < 0:

983 raise ValueError(_negative_)

984

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

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

987 if min(s) < EPS02:

988 raise ValueError(_or(_coincident_, _colinear_))

989

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

991 if min(q) < EPS02:

992 raise ValueError(_coincident_)

993

994 r1 = s2 * q3 / s3 # s2!

995 r2 = s1 * q3 / s3 # s1!

996 Qs = _F1(*q) # == hypot2_(a, b, c)

997 Ss = _F1(*s) # == fsum1(s)

998 s += (Qs * _0_5), # tuple!

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

1000 r3 = C0.fover(-s3)

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

1002 if d0 > EPS02: # > c0

1003 _xcallable(R3=R3)

1004 d0 = sqrt(d0)

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

1006 elif d0 < 0:

1007 raise ValueError(_negative_)

1008

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

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

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

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

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

1014

1015 except (TypeError, ValueError) as x:

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

1017

1018

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

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

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

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

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

1024 else:

1025 z = INT0

1026 return z

1027

1028# **) MIT License

1029#

1030# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.

1031#

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

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

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

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

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

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

1038#

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

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

1041#

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

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

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

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

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

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

1048# OTHER DEALINGS IN THE SOFTWARE.