C library for Geodesics: geodesic.c Source File (original) (raw)

27 #include <math.h>

28 #include <float.h>

30 #if !defined(__cplusplus)

31 #define nullptr 0

32 #endif

34 #define GEOGRAPHICLIB_GEODESIC_ORDER 6

35 #define nA1 GEOGRAPHICLIB_GEODESIC_ORDER

36 #define nC1 GEOGRAPHICLIB_GEODESIC_ORDER

37 #define nC1p GEOGRAPHICLIB_GEODESIC_ORDER

38 #define nA2 GEOGRAPHICLIB_GEODESIC_ORDER

39 #define nC2 GEOGRAPHICLIB_GEODESIC_ORDER

40 #define nA3 GEOGRAPHICLIB_GEODESIC_ORDER

41 #define nA3x nA3

42 #define nC3 GEOGRAPHICLIB_GEODESIC_ORDER

43 #define nC3x ((nC3 * (nC3 - 1)) / 2)

44 #define nC4 GEOGRAPHICLIB_GEODESIC_ORDER

45 #define nC4x ((nC4 * (nC4 + 1)) / 2)

46 #define nC (GEOGRAPHICLIB_GEODESIC_ORDER + 1)

48 typedef int boolx;

49 enum booly { FALSE = 0, TRUE = 1 };

53 enum dms { qd = 90, hd = 2 * qd, td = 2 * hd };

55 static unsigned init = 0;

56 static unsigned digits, maxit1, maxit2;

57 static double epsilon, realmin, pi, degree, NaN,

58 tiny, tol0, tol1, tol2, tolb, xthresh;

60 static void Init(void) {

61 if (!init) {

62 digits = DBL_MANT_DIG;

63 epsilon = DBL_EPSILON;

64 realmin = DBL_MIN;

65 #if defined(M_PI)

66 pi = M_PI;

67 #else

68 pi = atan2(0.0, -1.0);

69 #endif

70 maxit1 = 20;

71 maxit2 = maxit1 + digits + 10;

72 tiny = sqrt(realmin);

73 tol0 = epsilon;

77 tol1 = 200 * tol0;

78 tol2 = sqrt(tol0);

80 tolb = tol0;

81 xthresh = 1000 * tol2;

82 degree = pi/hd;

83 NaN = nan("0");

84 init = 1;

85 }

86 }

88 enum captype {

89 CAP_NONE = 0U,

90 CAP_C1 = 1U<<0,

91 CAP_C1p = 1U<<1,

92 CAP_C2 = 1U<<2,

93 CAP_C3 = 1U<<3,

94 CAP_C4 = 1U<<4,

95 CAP_ALL = 0x1FU,

96 OUT_ALL = 0x7F80U

97 };

99 static double sq(double x) { return x * x; }

100

101 static double sumx(double u, double v, double* t) {

102 volatile double s = u + v;

103 volatile double up = s - v;

104 volatile double vpp = s - up;

105 up -= u;

106 vpp -= v;

107 if (t) *t = s != 0 ? 0 - (up + vpp) : s;

108

109

110

111 return s;

112 }

113

114 static double polyvalx(int N, const double p[], double x) {

115 double y = N < 0 ? 0 : *p++;

116 while (--N >= 0) y = y * x + *p++;

117 return y;

118 }

119

120 static void swapx(double* x, double* y)

121 { double t = *x; *x = *y; *y = t; }

122

123 static void norm2(double* sinx, double* cosx) {

124 #if defined(_MSC_VER) && defined(_M_IX86)

125

126

127

128

129

130

131

132

133 double r = sqrt(*sinx * *sinx + *cosx * *cosx);

134 #else

135 double r = hypot(*sinx, *cosx);

136 #endif

137 *sinx /= r;

138 *cosx /= r;

139 }

140

141 static double AngNormalize(double x) {

142 double y = remainder(x, (double)td);

143 return fabs(y) == hd ? copysign((double)hd, x) : y;

144 }

145

146 static double LatFix(double x)

147 { return fabs(x) > qd ? NaN : x; }

148

149 static double AngDiff(double x, double y, double* e) {

150

151

152 double t, d = sumx(remainder(-x, (double)td), remainder( y, (double)td), &t);

153

154

155 d = sumx(remainder(d, (double)td), t, &t);

156

157 if (d == 0 || fabs(d) == hd)

158

159

160 d = copysign(d, t == 0 ? y - x : -t);

161 if (e) *e = t;

162 return d;

163 }

164

165 static double AngRound(double x) {

166

167 const double z = 1.0/16.0;

168 volatile double y = fabs(x);

169 volatile double w = z - y;

170

171 y = w > 0 ? z - w : y;

172 return copysign(y, x);

173 }

174

175 static void sincosdx(double x, double* sinx, double* cosx) {

176

177

178 double r, s, c; int q = 0;

179 r = remquo(x, (double)qd, &q);

180

181 r *= degree;

182

183 s = sin(r); c = cos(r);

184 switch ((unsigned)q & 3U) {

185 case 0U: *sinx = s; *cosx = c; break;

186 case 1U: *sinx = c; *cosx = -s; break;

187 case 2U: *sinx = -s; *cosx = -c; break;

188 default: *sinx = -c; *cosx = s; break;

189 }

190

191 *cosx += 0;

192

193 if (*sinx == 0) *sinx = copysign(*sinx, x);

194 }

195

196 static void sincosde(double x, double t, double* sinx, double* cosx) {

197

198

199 double r, s, c; int q = 0;

200 r = AngRound(remquo(x, (double)qd, &q) + t);

201

202 r *= degree;

203

204 s = sin(r); c = cos(r);

205 switch ((unsigned)q & 3U) {

206 case 0U: *sinx = s; *cosx = c; break;

207 case 1U: *sinx = c; *cosx = -s; break;

208 case 2U: *sinx = -s; *cosx = -c; break;

209 default: *sinx = -c; *cosx = s; break;

210 }

211

212 *cosx += 0;

213

214 if (*sinx == 0) *sinx = copysign(*sinx, x);

215 }

216

217 static double atan2dx(double y, double x) {

218

219

220

221

222 int q = 0; double ang;

223 if (fabs(y) > fabs(x)) { swapx(&x, &y); q = 2; }

224 if (signbit(x)) { x = -x; ++q; }

225

226 ang = atan2(y, x) / degree;

227 switch (q) {

228 case 1: ang = copysign((double)hd, y) - ang; break;

229 case 2: ang = qd - ang; break;

230 case 3: ang = -qd + ang; break;

231 default: break;

232 }

233 return ang;

234 }

235

239 static double SinCosSeries(boolx sinp,

240 double sinx, double cosx,

241 const double c[], int n);

242 static void Lengths(const struct geod_geodesic* g,

243 double eps, double sig12,

244 double ssig1, double csig1, double dn1,

245 double ssig2, double csig2, double dn2,

246 double cbet1, double cbet2,

247 double* ps12b, double* pm12b, double* pm0,

248 double* pM12, double* pM21,

249

250 double Ca[]);

251 static double Astroid(double x, double y);

252 static double InverseStart(const struct geod_geodesic* g,

253 double sbet1, double cbet1, double dn1,

254 double sbet2, double cbet2, double dn2,

255 double lam12, double slam12, double clam12,

256 double* psalp1, double* pcalp1,

257

258 double* psalp2, double* pcalp2,

259

260 double* pdnm,

261

262 double Ca[]);

263 static double Lambda12(const struct geod_geodesic* g,

264 double sbet1, double cbet1, double dn1,

265 double sbet2, double cbet2, double dn2,

266 double salp1, double calp1,

267 double slam120, double clam120,

268 double* psalp2, double* pcalp2,

269 double* psig12,

270 double* pssig1, double* pcsig1,

271 double* pssig2, double* pcsig2,

272 double* peps,

273 double* pdomg12,

274 boolx diffp, double* pdlam12,

275

276 double Ca[]);

277 static double A3f(const struct geod_geodesic* g, double eps);

278 static void C3f(const struct geod_geodesic* g, double eps, double c[]);

279 static void C4f(const struct geod_geodesic* g, double eps, double c[]);

280 static double A1m1f(double eps);

281 static void C1f(double eps, double c[]);

282 static void C1pf(double eps, double c[]);

283 static double A2m1f(double eps);

284 static void C2f(double eps, double c[]);

285 static int transit(double lon1, double lon2);

286 static int transitdirect(double lon1, double lon2);

287 static void accini(double s[]);

288 static void acccopy(const double s[], double t[]);

289 static void accadd(double s[], double y);

290 static double accsum(const double s[], double y);

291 static void accneg(double s[]);

292 static void accrem(double s[], double y);

293 static double areareduceA(double area[], double area0,

294 int crossings, boolx reverse, boolx sign);

295 static double areareduceB(double area, double area0,

296 int crossings, boolx reverse, boolx sign);

297

299 if (!init) Init();

300 g->a = a;

301 g->f = f;

302 g->f1 = 1 - g->f;

303 g->e2 = g->f * (2 - g->f);

304 g->ep2 = g->e2 / sq(g->f1);

305 g->n = g->f / ( 2 - g->f);

306 g->b = g->a * g->f1;

307 g->c2 = (sq(g->a) + sq(g->b) *

308 (g->e2 == 0 ? 1 :

309 (g->e2 > 0 ? atanh(sqrt(g->e2)) : atan(sqrt(-g->e2))) /

310 sqrt(fabs(g->e2))))/2;

311

312

313

314

315

316

317

318

319

320 g->etol2 = 0.1 * tol2 /

321 sqrt( fmax(0.001, fabs(g->f)) * fmin(1.0, 1 - g->f/2) / 2 );

322

323 A3coeff(g);

324 C3coeff(g);

325 C4coeff(g);

326 }

327

330 double lat1, double lon1,

331 double azi1, double salp1, double calp1,

332 unsigned caps) {

333 double cbet1, sbet1, eps;

334 l->a = g->a;

335 l->f = g->f;

336 l->b = g->b;

337 l->c2 = g->c2;

338 l->f1 = g->f1;

339

341

343

344 l->lat1 = LatFix(lat1);

345 l->lon1 = lon1;

346 l->azi1 = azi1;

347 l->salp1 = salp1;

348 l->calp1 = calp1;

349

350 sincosdx(AngRound(l->lat1), &sbet1, &cbet1); sbet1 *= l->f1;

351

352 norm2(&sbet1, &cbet1); cbet1 = fmax(tiny, cbet1);

353 l->dn1 = sqrt(1 + g->ep2 * sq(sbet1));

354

355

356 l->salp0 = l->salp1 * cbet1;

357

358

359 l->calp0 = hypot(l->calp1, l->salp1 * sbet1);

360

361

362

363

364

365

366

367

368

369 l->ssig1 = sbet1; l->somg1 = l->salp0 * sbet1;

370 l->csig1 = l->comg1 = sbet1 != 0 || l->calp1 != 0 ? cbet1 * l->calp1 : 1;

371 norm2(&l->ssig1, &l->csig1);

372

373

374 l->k2 = sq(l->calp0) * g->ep2;

375 eps = l->k2 / (2 * (1 + sqrt(1 + l->k2)) + l->k2);

376

377 if (l->caps & CAP_C1) {

378 double s, c;

379 l->A1m1 = A1m1f(eps);

380 C1f(eps, l->C1a);

381 l->B11 = SinCosSeries(TRUE, l->ssig1, l->csig1, l->C1a, nC1);

382 s = sin(l->B11); c = cos(l->B11);

383

384 l->stau1 = l->ssig1 * c + l->csig1 * s;

385 l->ctau1 = l->csig1 * c - l->ssig1 * s;

386

387

388 }

389

390 if (l->caps & CAP_C1p)

391 C1pf(eps, l->C1pa);

392

393 if (l->caps & CAP_C2) {

394 l->A2m1 = A2m1f(eps);

395 C2f(eps, l->C2a);

396 l->B21 = SinCosSeries(TRUE, l->ssig1, l->csig1, l->C2a, nC2);

397 }

398

399 if (l->caps & CAP_C3) {

400 C3f(g, eps, l->C3a);

401 l->A3c = -l->f * l->salp0 * A3f(g, eps);

402 l->B31 = SinCosSeries(TRUE, l->ssig1, l->csig1, l->C3a, nC3-1);

403 }

404

405 if (l->caps & CAP_C4) {

406 C4f(g, eps, l->C4a);

407

408 l->A4 = sq(l->a) * l->calp0 * l->salp0 * g->e2;

409 l->B41 = SinCosSeries(FALSE, l->ssig1, l->csig1, l->C4a, nC4);

410 }

411

412 l->a13 = l->s13 = NaN;

413 }

414

417 double lat1, double lon1, double azi1, unsigned caps) {

418 double salp1, calp1;

419 azi1 = AngNormalize(azi1);

420

421 sincosdx(AngRound(azi1), &salp1, &calp1);

422 geod_lineinit_int(l, g, lat1, lon1, azi1, salp1, calp1, caps);

423 }

424

427 double lat1, double lon1, double azi1,

428 unsigned flags, double s12_a12,

429 unsigned caps) {

432 }

433

436 double lat1, double lon1, double azi1,

437 double s12, unsigned caps) {

439 }

440

442 unsigned flags, double s12_a12,

443 double* plat2, double* plon2, double* pazi2,

444 double* ps12, double* pm12,

445 double* pM12, double* pM21,

446 double* pS12) {

447 double lat2 = 0, lon2 = 0, azi2 = 0, s12 = 0,

448 m12 = 0, M12 = 0, M21 = 0, S12 = 0;

449

450 double sig12, ssig12, csig12, B12 = 0, AB1 = 0;

451 double omg12, lam12, lon12;

452 double ssig2, csig2, sbet2, cbet2, somg2, comg2, salp2, calp2, dn2;

453 unsigned outmask =

461

462 outmask &= l->caps & OUT_ALL;

464

465 return NaN;

466

468

469 sig12 = s12_a12 * degree;

470 sincosdx(s12_a12, &ssig12, &csig12);

471 } else {

472

473 double

474 tau12 = s12_a12 / (l->b * (1 + l->A1m1)),

475 s = sin(tau12),

476 c = cos(tau12);

477

478 B12 = - SinCosSeries(TRUE,

479 l->stau1 * c + l->ctau1 * s,

480 l->ctau1 * c - l->stau1 * s,

481 l->C1pa, nC1p);

482 sig12 = tau12 - (B12 - l->B11);

483 ssig12 = sin(sig12); csig12 = cos(sig12);

484 if (fabs(l->f) > 0.01) {

485

486

487

488

489

490

491

492

493

494

495

496

497

498

499

500

501

502

503

504

505

506 double serr;

507 ssig2 = l->ssig1 * csig12 + l->csig1 * ssig12;

508 csig2 = l->csig1 * csig12 - l->ssig1 * ssig12;

509 B12 = SinCosSeries(TRUE, ssig2, csig2, l->C1a, nC1);

510 serr = (1 + l->A1m1) * (sig12 + (B12 - l->B11)) - s12_a12 / l->b;

511 sig12 = sig12 - serr / sqrt(1 + l->k2 * sq(ssig2));

512 ssig12 = sin(sig12); csig12 = cos(sig12);

513

514 }

515 }

516

517

518 ssig2 = l->ssig1 * csig12 + l->csig1 * ssig12;

519 csig2 = l->csig1 * csig12 - l->ssig1 * ssig12;

520 dn2 = sqrt(1 + l->k2 * sq(ssig2));

523 B12 = SinCosSeries(TRUE, ssig2, csig2, l->C1a, nC1);

524 AB1 = (1 + l->A1m1) * (B12 - l->B11);

525 }

526

527 sbet2 = l->calp0 * ssig2;

528

529 cbet2 = hypot(l->salp0, l->calp0 * csig2);

530 if (cbet2 == 0)

531

532 cbet2 = csig2 = tiny;

533

534 salp2 = l->salp0; calp2 = l->calp0 * csig2;

535

538 l->b * ((1 + l->A1m1) * sig12 + AB1) :

539 s12_a12;

540

542 double E = copysign(1, l->salp0);

543

544 somg2 = l->salp0 * ssig2; comg2 = csig2;

545

547 ? E * (sig12

548 - (atan2( ssig2, csig2) - atan2( l->ssig1, l->csig1))

549 + (atan2(E * somg2, comg2) - atan2(E * l->somg1, l->comg1)))

550 : atan2(somg2 * l->comg1 - comg2 * l->somg1,

551 comg2 * l->comg1 + somg2 * l->somg1);

552 lam12 = omg12 + l->A3c *

553 ( sig12 + (SinCosSeries(TRUE, ssig2, csig2, l->C3a, nC3-1)

554 - l->B31));

555 lon12 = lam12 / degree;

557 AngNormalize(AngNormalize(l->lon1) + AngNormalize(lon12));

558 }

559

561 lat2 = atan2dx(sbet2, l->f1 * cbet2);

562

564 azi2 = atan2dx(salp2, calp2);

565

567 double

568 B22 = SinCosSeries(TRUE, ssig2, csig2, l->C2a, nC2),

569 AB2 = (1 + l->A2m1) * (B22 - l->B21),

570 J12 = (l->A1m1 - l->A2m1) * sig12 + (AB1 - AB2);

572

573

574 m12 = l->b * ((dn2 * (l->csig1 * ssig2) - l->dn1 * (l->ssig1 * csig2))

575 - l->csig1 * csig2 * J12);

577 double t = l->k2 * (ssig2 - l->ssig1) * (ssig2 + l->ssig1) /

578 (l->dn1 + dn2);

579 M12 = csig12 + (t * ssig2 - csig2 * J12) * l->ssig1 / l->dn1;

580 M21 = csig12 - (t * l->ssig1 - l->csig1 * J12) * ssig2 / dn2;

581 }

582 }

583

585 double

586 B42 = SinCosSeries(FALSE, ssig2, csig2, l->C4a, nC4);

587 double salp12, calp12;

588 if (l->calp0 == 0 || l->salp0 == 0) {

589

590 salp12 = salp2 * l->calp1 - calp2 * l->salp1;

591 calp12 = calp2 * l->calp1 + salp2 * l->salp1;

592 } else {

593

594

595

596

597

598

599

600

601 salp12 = l->calp0 * l->salp0 *

602 (csig12 <= 0 ? l->csig1 * (1 - csig12) + ssig12 * l->ssig1 :

603 ssig12 * (l->csig1 * ssig12 / (1 + csig12) + l->ssig1));

604 calp12 = sq(l->salp0) + sq(l->calp0) * l->csig1 * csig2;

605 }

606 S12 = l->c2 * atan2(salp12, calp12) + l->A4 * (B42 - l->B41);

607 }

608

609

610

611

612

613

614

615

616

618 *plat2 = lat2;

620 *plon2 = lon2;

622 *pazi2 = azi2;

624 *ps12 = s12;

626 *pm12 = m12;

628 if (pM12) *pM12 = M12;

629 if (pM21) *pM21 = M21;

630 }

631 if ((outmask & GEOD_AREA) && pS12)

632 *pS12 = S12;

633

634 return (flags & GEOD_ARCMODE) ? s12_a12 : sig12 / degree;

635 }

636

638 l->s13 = s13;

640 nullptr, nullptr, nullptr, nullptr, nullptr);

641 }

642

643 static void geod_setarc(struct geod_geodesicline* l, double a13) {

644 l->a13 = a13; l->s13 = NaN;

646 nullptr, nullptr, nullptr, nullptr);

647 }

648

650 unsigned flags, double s13_a13) {

652 geod_setarc(l, s13_a13) :

654 }

655

657 double* plat2, double* plon2, double* pazi2) {

659 nullptr, nullptr, nullptr, nullptr, nullptr);

660 }

661

663 double lat1, double lon1, double azi1,

664 unsigned flags, double s12_a12,

665 double* plat2, double* plon2, double* pazi2,

666 double* ps12, double* pm12, double* pM12, double* pM21,

667 double* pS12) {

669 unsigned outmask =

677

679

680 outmask |

683 plat2, plon2, pazi2, ps12, pm12, pM12, pM21, pS12);

684 }

685

688 double s12,

689 double* plat2, double* plon2, double* pazi2) {

691 nullptr, nullptr, nullptr, nullptr, nullptr);

692 }

693

694 static double geod_geninverse_int(const struct geod_geodesic* g,

696 double lat2, double lon2,

697 double* ps12,

698 double* psalp1, double* pcalp1,

699 double* psalp2, double* pcalp2,

700 double* pm12, double* pM12, double* pM21,

701 double* pS12) {

702 double s12 = 0, m12 = 0, M12 = 0, M21 = 0, S12 = 0;

703 double lon12, lon12s;

704 int latsign, lonsign, swapp;

705 double sbet1, cbet1, sbet2, cbet2, s12x = 0, m12x = 0;

706 double dn1, dn2, lam12, slam12, clam12;

707 double a12 = 0, sig12, calp1 = 0, salp1 = 0, calp2 = 0, salp2 = 0;

708 double Ca[nC];

709 boolx meridian;

710

711 double omg12 = 0, somg12 = 2, comg12 = 0;

712

713 unsigned outmask =

718

719 outmask &= OUT_ALL;

720

721

722

723 lon12 = AngDiff(lon1, lon2, &lon12s);

724

725 lonsign = signbit(lon12) ? -1 : 1;

726 lon12 *= lonsign; lon12s *= lonsign;

727 lam12 = lon12 * degree;

728

729 sincosde(lon12, lon12s, &slam12, &clam12);

730 lon12s = (hd - lon12) - lon12s;

731

732

733 lat1 = AngRound(LatFix(lat1));

734 lat2 = AngRound(LatFix(lat2));

735

736

737 swapp = fabs(lat1) < fabs(lat2) || lat2 != lat2 ? -1 : 1;

738 if (swapp < 0) {

739 lonsign *= -1;

740 swapx(&lat1, &lat2);

741 }

742

743 latsign = signbit(lat1) ? 1 : -1;

744 lat1 *= latsign;

745 lat2 *= latsign;

746

747

748

749

750

751

752

753

754

755

756

757

758 sincosdx(lat1, &sbet1, &cbet1); sbet1 *= g->f1;

759

760 norm2(&sbet1, &cbet1); cbet1 = fmax(tiny, cbet1);

761

762 sincosdx(lat2, &sbet2, &cbet2); sbet2 *= g->f1;

763

764 norm2(&sbet2, &cbet2); cbet2 = fmax(tiny, cbet2);

765

766

767

768

769

770

771

772

773

774 if (cbet1 < -sbet1) {

775 if (cbet2 == cbet1)

776 sbet2 = copysign(sbet1, sbet2);

777 } else {

778 if (fabs(sbet2) == -sbet1)

779 cbet2 = cbet1;

780 }

781

782 dn1 = sqrt(1 + g->ep2 * sq(sbet1));

783 dn2 = sqrt(1 + g->ep2 * sq(sbet2));

784

785 meridian = lat1 == -qd || slam12 == 0;

786

787 if (meridian) {

788

789

790

791

792 double ssig1, csig1, ssig2, csig2;

793 calp1 = clam12; salp1 = slam12;

794 calp2 = 1; salp2 = 0;

795

796

797 ssig1 = sbet1; csig1 = calp1 * cbet1;

798 ssig2 = sbet2; csig2 = calp2 * cbet2;

799

800

801 sig12 = atan2(fmax(0.0, csig1 * ssig2 - ssig1 * csig2) + 0,

802 csig1 * csig2 + ssig1 * ssig2);

803 Lengths(g, g->n, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,

804 cbet1, cbet2, &s12x, &m12x, nullptr,

807 Ca);

808

809

810

811

812

813

814

815 if (sig12 < 1 || m12x >= 0) {

816

817 if (sig12 < 3 * tiny ||

818

819 (sig12 < tol0 && (s12x < 0 || m12x < 0)))

820 sig12 = m12x = s12x = 0;

821 m12x *= g->b;

822 s12x *= g->b;

823 a12 = sig12 / degree;

824 } else

825

826 meridian = FALSE;

827 }

828

829 if (!meridian &&

830 sbet1 == 0 &&

831

832 (g->f <= 0 || lon12s >= g->f * hd)) {

833

834

835 calp1 = calp2 = 0; salp1 = salp2 = 1;

836 s12x = g->a * lam12;

837 sig12 = omg12 = lam12 / g->f1;

838 m12x = g->b * sin(sig12);

840 M12 = M21 = cos(sig12);

841 a12 = lon12 / g->f1;

842

843 } else if (!meridian) {

844

845

846

847

848

849 double dnm = 0;

850 sig12 = InverseStart(g, sbet1, cbet1, dn1, sbet2, cbet2, dn2,

851 lam12, slam12, clam12,

852 &salp1, &calp1, &salp2, &calp2, &dnm,

853 Ca);

854

855 if (sig12 >= 0) {

856

857 s12x = sig12 * g->b * dnm;

858 m12x = sq(dnm) * g->b * sin(sig12 / dnm);

860 M12 = M21 = cos(sig12 / dnm);

861 a12 = sig12 / degree;

862 omg12 = lam12 / (g->f1 * dnm);

863 } else {

864

865

866

867

868

869

870

871

872

873

874

875

876 double ssig1 = 0, csig1 = 0, ssig2 = 0, csig2 = 0, eps = 0, domg12 = 0;

877 unsigned numit = 0;

878

879 double salp1a = tiny, calp1a = 1, salp1b = tiny, calp1b = -1;

880 boolx tripn = FALSE;

881 boolx tripb = FALSE;

882 for (;; ++numit) {

883

884

885 double dv = 0,

886 v = Lambda12(g, sbet1, cbet1, dn1, sbet2, cbet2, dn2, salp1, calp1,

887 slam12, clam12,

888 &salp2, &calp2, &sig12, &ssig1, &csig1, &ssig2, &csig2,

889 &eps, &domg12, numit < maxit1, &dv, Ca);

890 if (tripb ||

891

892 !(fabs(v) >= (tripn ? 8 : 1) * tol0) ||

893

894 numit == maxit2)

895 break;

896

897 if (v > 0 && (numit > maxit1 || calp1/salp1 > calp1b/salp1b))

899 else if (v < 0 && (numit > maxit1 || calp1/salp1 < calp1a/salp1a))

901 if (numit < maxit1 && dv > 0) {

902 double

903 dalp1 = -v/dv;

904 if (fabs(dalp1) < pi) {

905 double

906 sdalp1 = sin(dalp1), cdalp1 = cos(dalp1),

907 nsalp1 = salp1 * cdalp1 + calp1 * sdalp1;

908 if (nsalp1 > 0) {

912

913

914

915 tripn = fabs(v) <= 16 * tol0;

916 continue;

917 }

918 }

919 }

920

921

922

923

924

925

926

927

928 salp1 = (salp1a + salp1b)/2;

929 calp1 = (calp1a + calp1b)/2;

931 tripn = FALSE;

932 tripb = (fabs(salp1a - salp1) + (calp1a - calp1) < tolb ||

933 fabs(salp1 - salp1b) + (calp1 - calp1b) < tolb);

934 }

935 Lengths(g, eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,

936 cbet1, cbet2, &s12x, &m12x, nullptr,

939 m12x *= g->b;

940 s12x *= g->b;

941 a12 = sig12 / degree;

943

944 double sdomg12 = sin(domg12), cdomg12 = cos(domg12);

945 somg12 = slam12 * cdomg12 - clam12 * sdomg12;

946 comg12 = clam12 * cdomg12 + slam12 * sdomg12;

947 }

948 }

949 }

950

952 s12 = 0 + s12x;

953

955 m12 = 0 + m12x;

956

958 double

959

960 salp0 = salp1 * cbet1,

961 calp0 = hypot(calp1, salp1 * sbet1);

962 double alp12;

963 if (calp0 != 0 && salp0 != 0) {

964 double

965

966 ssig1 = sbet1, csig1 = calp1 * cbet1,

967 ssig2 = sbet2, csig2 = calp2 * cbet2,

968 k2 = sq(calp0) * g->ep2,

969 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2),

970

971 A4 = sq(g->a) * calp0 * salp0 * g->e2;

972 double B41, B42;

973 norm2(&ssig1, &csig1);

974 norm2(&ssig2, &csig2);

975 C4f(g, eps, Ca);

976 B41 = SinCosSeries(FALSE, ssig1, csig1, Ca, nC4);

977 B42 = SinCosSeries(FALSE, ssig2, csig2, Ca, nC4);

978 S12 = A4 * (B42 - B41);

979 } else

980

981 S12 = 0;

982

983 if (!meridian && somg12 == 2) {

984 somg12 = sin(omg12); comg12 = cos(omg12);

985 }

986

987 if (!meridian &&

988

989 comg12 > -0.7071 &&

990 sbet2 - sbet1 < 1.75) {

991

992

993

994 double

995 domg12 = 1 + comg12, dbet1 = 1 + cbet1, dbet2 = 1 + cbet2;

996 alp12 = 2 * atan2( somg12 * ( sbet1 * dbet2 + sbet2 * dbet1 ),

997 domg12 * ( sbet1 * sbet2 + dbet1 * dbet2 ) );

998 } else {

999

1000 double

1001 salp12 = salp2 * calp1 - calp2 * salp1,

1002 calp12 = calp2 * calp1 + salp2 * salp1;

1003

1004

1005

1006

1007 if (salp12 == 0 && calp12 < 0) {

1008 salp12 = tiny * calp1;

1009 calp12 = -1;

1010 }

1011 alp12 = atan2(salp12, calp12);

1012 }

1013 S12 += g->c2 * alp12;

1014 S12 *= swapp * lonsign * latsign;

1015

1016 S12 += 0;

1017 }

1018

1019

1020 if (swapp < 0) {

1021 swapx(&salp1, &salp2);

1022 swapx(&calp1, &calp2);

1024 swapx(&M12, &M21);

1025 }

1026

1027 salp1 *= swapp * lonsign; calp1 *= swapp * latsign;

1028 salp2 *= swapp * lonsign; calp2 *= swapp * latsign;

1029

1030 if (psalp1) *psalp1 = salp1;

1031 if (pcalp1) *pcalp1 = calp1;

1032 if (psalp2) *psalp2 = salp2;

1033 if (pcalp2) *pcalp2 = calp2;

1034

1036 *ps12 = s12;

1038 *pm12 = m12;

1040 if (pM12) *pM12 = M12;

1041 if (pM21) *pM21 = M21;

1042 }

1044 *pS12 = S12;

1045

1046

1047 return a12;

1048 }

1049

1051 double lat1, double lon1, double lat2, double lon2,

1052 double* ps12, double* pazi1, double* pazi2,

1053 double* pm12, double* pM12, double* pM21,

1054 double* pS12) {

1056 a12 = geod_geninverse_int(g, lat1, lon1, lat2, lon2, ps12,

1058 pm12, pM12, pM21, pS12);

1059 if (pazi1) *pazi1 = atan2dx(salp1, calp1);

1060 if (pazi2) *pazi2 = atan2dx(salp2, calp2);

1061 return a12;

1062 }

1063

1066 double lat1, double lon1, double lat2, double lon2,

1067 unsigned caps) {

1069 a12 = geod_geninverse_int(g, lat1, lon1, lat2, lon2, nullptr,

1071 nullptr, nullptr, nullptr, nullptr),

1074

1077 geod_setarc(l, a12);

1078 }

1079

1081 double lat1, double lon1, double lat2, double lon2,

1082 double* ps12, double* pazi1, double* pazi2) {

1084 nullptr, nullptr, nullptr, nullptr);

1085 }

1086

1087 double SinCosSeries(boolx sinp, double sinx, double cosx,

1088 const double c[], int n) {

1089

1090

1091

1092

1093

1094 double ar, y0, y1;

1095 c += (n + sinp);

1096 ar = 2 * (cosx - sinx) * (cosx + sinx);

1097 y0 = (n & 1) ? *--c : 0; y1 = 0;

1098

1099 n /= 2;

1100 while (n--) {

1101

1102 y1 = ar * y0 - y1 + *--c;

1103 y0 = ar * y1 - y0 + *--c;

1104 }

1105 return sinp

1106 ? 2 * sinx * cosx * y0

1107 : cosx * (y0 - y1);

1108 }

1109

1111 double eps, double sig12,

1112 double ssig1, double csig1, double dn1,

1113 double ssig2, double csig2, double dn2,

1114 double cbet1, double cbet2,

1115 double* ps12b, double* pm12b, double* pm0,

1116 double* pM12, double* pM21,

1117

1118 double Ca[]) {

1119 double m0 = 0, J12 = 0, A1 = 0, A2 = 0;

1120 double Cb[nC];

1121

1122

1123

1124 boolx redlp = pm12b || pm0 || pM12 || pM21;

1125 if (ps12b || redlp) {

1126 A1 = A1m1f(eps);

1127 C1f(eps, Ca);

1128 if (redlp) {

1129 A2 = A2m1f(eps);

1130 C2f(eps, Cb);

1131 m0 = A1 - A2;

1132 A2 = 1 + A2;

1133 }

1134 A1 = 1 + A1;

1135 }

1136 if (ps12b) {

1137 double B1 = SinCosSeries(TRUE, ssig2, csig2, Ca, nC1) -

1138 SinCosSeries(TRUE, ssig1, csig1, Ca, nC1);

1139

1140 *ps12b = A1 * (sig12 + B1);

1141 if (redlp) {

1142 double B2 = SinCosSeries(TRUE, ssig2, csig2, Cb, nC2) -

1143 SinCosSeries(TRUE, ssig1, csig1, Cb, nC2);

1144 J12 = m0 * sig12 + (A1 * B1 - A2 * B2);

1145 }

1146 } else if (redlp) {

1147

1148 int l;

1149 for (l = 1; l <= nC2; ++l)

1150 Cb[l] = A1 * Ca[l] - A2 * Cb[l];

1151 J12 = m0 * sig12 + (SinCosSeries(TRUE, ssig2, csig2, Cb, nC2) -

1152 SinCosSeries(TRUE, ssig1, csig1, Cb, nC2));

1153 }

1154 if (pm0) *pm0 = m0;

1155 if (pm12b)

1156

1157

1158

1159 *pm12b = dn2 * (csig1 * ssig2) - dn1 * (ssig1 * csig2) -

1160 csig1 * csig2 * J12;

1161 if (pM12 || pM21) {

1162 double csig12 = csig1 * csig2 + ssig1 * ssig2;

1163 double t = g->ep2 * (cbet1 - cbet2) * (cbet1 + cbet2) / (dn1 + dn2);

1164 if (pM12)

1165 *pM12 = csig12 + (t * ssig2 - csig2 * J12) * ssig1 / dn1;

1166 if (pM21)

1167 *pM21 = csig12 - (t * ssig1 - csig1 * J12) * ssig2 / dn2;

1168 }

1169 }

1170

1171 double Astroid(double x, double y) {

1172

1173

1174 double k;

1175 double

1176 p = sq(x),

1177 q = sq(y),

1178 r = (p + q - 1) / 6;

1179 if ( !(q == 0 && r <= 0) ) {

1180 double

1181

1182

1183 S = p * q / 4,

1184 r2 = sq(r),

1185 r3 = r * r2,

1186

1187

1188 disc = S * (S + 2 * r3);

1189 double u = r;

1190 double v, uv, w;

1191 if (disc >= 0) {

1192 double T3 = S + r3, T;

1193

1194

1195

1196 T3 += T3 < 0 ? -sqrt(disc) : sqrt(disc);

1197

1198 T = cbrt(T3);

1199

1200 u += T + (T != 0 ? r2 / T : 0);

1201 } else {

1202

1203 double ang = atan2(sqrt(-disc), -(S + r3));

1204

1205

1206 u += 2 * r * cos(ang / 3);

1207 }

1208 v = sqrt(sq(u) + q);

1209

1210 uv = u < 0 ? q / (v - u) : u + v;

1211 w = (uv - q) / (2 * v);

1212

1213

1214 k = uv / (sqrt(uv + sq(w)) + w);

1215 } else {

1216

1217

1218 k = 0;

1219 }

1220 return k;

1221 }

1222

1223 double InverseStart(const struct geod_geodesic* g,

1224 double sbet1, double cbet1, double dn1,

1225 double sbet2, double cbet2, double dn2,

1226 double lam12, double slam12, double clam12,

1227 double* psalp1, double* pcalp1,

1228

1229 double* psalp2, double* pcalp2,

1230

1231 double* pdnm,

1232

1233 double Ca[]) {

1234 double salp1 = 0, calp1 = 0, salp2 = 0, calp2 = 0, dnm = 0;

1235

1236

1237

1238

1239 double

1240 sig12 = -1,

1241

1242 sbet12 = sbet2 * cbet1 - cbet2 * sbet1,

1243 cbet12 = cbet2 * cbet1 + sbet2 * sbet1;

1244 double sbet12a;

1245 boolx shortline = cbet12 >= 0 && sbet12 < 0.5 && cbet2 * lam12 < 0.5;

1246 double somg12, comg12, ssig12, csig12;

1247 sbet12a = sbet2 * cbet1 + cbet2 * sbet1;

1248 if (shortline) {

1249 double sbetm2 = sq(sbet1 + sbet2), omg12;

1250

1251

1252 sbetm2 /= sbetm2 + sq(cbet1 + cbet2);

1253 dnm = sqrt(1 + g->ep2 * sbetm2);

1254 omg12 = lam12 / (g->f1 * dnm);

1255 somg12 = sin(omg12); comg12 = cos(omg12);

1256 } else {

1257 somg12 = slam12; comg12 = clam12;

1258 }

1259

1260 salp1 = cbet2 * somg12;

1261 calp1 = comg12 >= 0 ?

1262 sbet12 + cbet2 * sbet1 * sq(somg12) / (1 + comg12) :

1263 sbet12a - cbet2 * sbet1 * sq(somg12) / (1 - comg12);

1264

1266 csig12 = sbet1 * sbet2 + cbet1 * cbet2 * comg12;

1267

1268 if (shortline && ssig12 < g->etol2) {

1269

1270 salp2 = cbet1 * somg12;

1271 calp2 = sbet12 - cbet1 * sbet2 *

1272 (comg12 >= 0 ? sq(somg12) / (1 + comg12) : 1 - comg12);

1273 norm2(&salp2, &calp2);

1274

1275 sig12 = atan2(ssig12, csig12);

1276 } else if (fabs(g->n) > 0.1 ||

1277 csig12 >= 0 ||

1278 ssig12 >= 6 * fabs(g->n) * pi * sq(cbet1)) {

1279

1280 } else {

1281

1282

1283 double x, y, lamscale, betscale;

1284 double lam12x = atan2(-slam12, -clam12);

1285 if (g->f >= 0) {

1286

1287 {

1288 double

1289 k2 = sq(sbet1) * g->ep2,

1290 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2);

1291 lamscale = g->f * cbet1 * A3f(g, eps) * pi;

1292 }

1293 betscale = lamscale * cbet1;

1294

1295 x = lam12x / lamscale;

1296 y = sbet12a / betscale;

1297 } else {

1298

1299 double

1300 cbet12a = cbet2 * cbet1 - sbet2 * sbet1,

1301 bet12a = atan2(sbet12a, cbet12a);

1302 double m12b, m0;

1303

1304

1305 Lengths(g, g->n, pi + bet12a,

1306 sbet1, -cbet1, dn1, sbet2, cbet2, dn2,

1307 cbet1, cbet2, nullptr, &m12b, &m0, nullptr, nullptr, Ca);

1308 x = -1 + m12b / (cbet1 * cbet2 * m0 * pi);

1309 betscale = x < -0.01 ? sbet12a / x :

1310 -g->f * sq(cbet1) * pi;

1311 lamscale = betscale / cbet1;

1312 y = lam12x / lamscale;

1313 }

1314

1315 if (y > -tol1 && x > -1 - xthresh) {

1316

1317 if (g->f >= 0) {

1319 } else {

1320 calp1 = fmax(x > -tol1 ? 0.0 : -1.0, x);

1322 }

1323 } else {

1324

1325

1326

1327

1328

1329

1330

1331

1332

1333

1334

1335

1336

1337

1338

1339

1340

1341

1342

1343

1344

1345

1346

1347

1348

1349

1350

1351

1352

1353

1354

1355

1356

1357

1358 double k = Astroid(x, y);

1359 double

1360 omg12a = lamscale * ( g->f >= 0 ? -x * k/(1 + k) : -y * (1 + k)/k );

1361 somg12 = sin(omg12a); comg12 = -cos(omg12a);

1362

1363 salp1 = cbet2 * somg12;

1364 calp1 = sbet12a - cbet2 * sbet1 * sq(somg12) / (1 - comg12);

1365 }

1366 }

1367

1368 if (!(salp1 <= 0))

1370 else {

1372 }

1373

1374 *psalp1 = salp1;

1375 *pcalp1 = calp1;

1376 if (shortline)

1377 *pdnm = dnm;

1378 if (sig12 >= 0) {

1379 *psalp2 = salp2;

1380 *pcalp2 = calp2;

1381 }

1382 return sig12;

1383 }

1384

1386 double sbet1, double cbet1, double dn1,

1387 double sbet2, double cbet2, double dn2,

1389 double slam120, double clam120,

1390 double* psalp2, double* pcalp2,

1391 double* psig12,

1392 double* pssig1, double* pcsig1,

1393 double* pssig2, double* pcsig2,

1394 double* peps,

1395 double* pdomg12,

1396 boolx diffp, double* pdlam12,

1397

1398 double Ca[]) {

1399 double salp2 = 0, calp2 = 0, sig12 = 0,

1400 ssig1 = 0, csig1 = 0, ssig2 = 0, csig2 = 0, eps = 0,

1401 domg12 = 0, dlam12 = 0;

1402 double salp0, calp0;

1403 double somg1, comg1, somg2, comg2, somg12, comg12, lam12;

1404 double B312, eta, k2;

1405

1406 if (sbet1 == 0 && calp1 == 0)

1407

1408

1410

1411

1412 salp0 = salp1 * cbet1;

1413 calp0 = hypot(calp1, salp1 * sbet1);

1414

1415

1416

1417 ssig1 = sbet1; somg1 = salp0 * sbet1;

1418 csig1 = comg1 = calp1 * cbet1;

1419 norm2(&ssig1, &csig1);

1420

1421

1422

1423

1424

1425

1426 salp2 = cbet2 != cbet1 ? salp0 / cbet2 : salp1;

1427

1428

1429

1430

1431 calp2 = cbet2 != cbet1 || fabs(sbet2) != -sbet1 ?

1432 sqrt(sq(calp1 * cbet1) +

1433 (cbet1 < -sbet1 ?

1434 (cbet2 - cbet1) * (cbet1 + cbet2) :

1435 (sbet1 - sbet2) * (sbet1 + sbet2))) / cbet2 :

1437

1438

1439 ssig2 = sbet2; somg2 = salp0 * sbet2;

1440 csig2 = comg2 = calp2 * cbet2;

1441 norm2(&ssig2, &csig2);

1442

1443

1444

1445 sig12 = atan2(fmax(0.0, csig1 * ssig2 - ssig1 * csig2) + 0,

1446 csig1 * csig2 + ssig1 * ssig2);

1447

1448

1449 somg12 = fmax(0.0, comg1 * somg2 - somg1 * comg2) + 0;

1450 comg12 = comg1 * comg2 + somg1 * somg2;

1451

1452 eta = atan2(somg12 * clam120 - comg12 * slam120,

1453 comg12 * clam120 + somg12 * slam120);

1454 k2 = sq(calp0) * g->ep2;

1455 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2);

1456 C3f(g, eps, Ca);

1457 B312 = (SinCosSeries(TRUE, ssig2, csig2, Ca, nC3-1) -

1458 SinCosSeries(TRUE, ssig1, csig1, Ca, nC3-1));

1459 domg12 = -g->f * A3f(g, eps) * salp0 * (sig12 + B312);

1460 lam12 = eta + domg12;

1461

1462 if (diffp) {

1463 if (calp2 == 0)

1464 dlam12 = - 2 * g->f1 * dn1 / sbet1;

1465 else {

1466 Lengths(g, eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,

1467 cbet1, cbet2, nullptr, &dlam12, nullptr, nullptr, nullptr, Ca);

1468 dlam12 *= g->f1 / (calp2 * cbet2);

1469 }

1470 }

1471

1472 *psalp2 = salp2;

1473 *pcalp2 = calp2;

1474 *psig12 = sig12;

1475 *pssig1 = ssig1;

1476 *pcsig1 = csig1;

1477 *pssig2 = ssig2;

1478 *pcsig2 = csig2;

1479 *peps = eps;

1480 *pdomg12 = domg12;

1481 if (diffp)

1482 *pdlam12 = dlam12;

1483

1484 return lam12;

1485 }

1486

1487 double A3f(const struct geod_geodesic* g, double eps) {

1488

1489 return polyvalx(nA3 - 1, g->A3x, eps);

1490 }

1491

1492 void C3f(const struct geod_geodesic* g, double eps, double c[]) {

1493

1494

1495 double mult = 1;

1496 int o = 0, l;

1497 for (l = 1; l < nC3; ++l) {

1498 int m = nC3 - l - 1;

1499 mult *= eps;

1500 c[l] = mult * polyvalx(m, g->C3x + o, eps);

1501 o += m + 1;

1502 }

1503 }

1504

1505 void C4f(const struct geod_geodesic* g, double eps, double c[]) {

1506

1507

1508 double mult = 1;

1509 int o = 0, l;

1510 for (l = 0; l < nC4; ++l) {

1511 int m = nC4 - l - 1;

1512 c[l] = mult * polyvalx(m, g->C4x + o, eps);

1513 o += m + 1;

1514 mult *= eps;

1515 }

1516 }

1517

1518

1519 double A1m1f(double eps) {

1520 static const double coeff[] = {

1521

1522 1, 4, 64, 0, 256,

1523 };

1524 int m = nA1/2;

1525 double t = polyvalx(m, coeff, sq(eps)) / coeff[m + 1];

1526 return (t + eps) / (1 - eps);

1527 }

1528

1529

1530 void C1f(double eps, double c[]) {

1531 static const double coeff[] = {

1532

1533 -1, 6, -16, 32,

1534

1535 -9, 64, -128, 2048,

1536

1537 9, -16, 768,

1538

1539 3, -5, 512,

1540

1541 -7, 1280,

1542

1543 -7, 2048,

1544 };

1545 double

1546 eps2 = sq(eps),

1547 d = eps;

1548 int o = 0, l;

1549 for (l = 1; l <= nC1; ++l) {

1550 int m = (nC1 - l) / 2;

1551 c[l] = d * polyvalx(m, coeff + o, eps2) / coeff[o + m + 1];

1552 o += m + 2;

1553 d *= eps;

1554 }

1555 }

1556

1557

1558 void C1pf(double eps, double c[]) {

1559 static const double coeff[] = {

1560

1561 205, -432, 768, 1536,

1562

1563 4005, -4736, 3840, 12288,

1564

1565 -225, 116, 384,

1566

1567 -7173, 2695, 7680,

1568

1569 3467, 7680,

1570

1571 38081, 61440,

1572 };

1573 double

1574 eps2 = sq(eps),

1575 d = eps;

1576 int o = 0, l;

1577 for (l = 1; l <= nC1p; ++l) {

1578 int m = (nC1p - l) / 2;

1579 c[l] = d * polyvalx(m, coeff + o, eps2) / coeff[o + m + 1];

1580 o += m + 2;

1581 d *= eps;

1582 }

1583 }

1584

1585

1586 double A2m1f(double eps) {

1587 static const double coeff[] = {

1588

1589 -11, -28, -192, 0, 256,

1590 };

1591 int m = nA2/2;

1592 double t = polyvalx(m, coeff, sq(eps)) / coeff[m + 1];

1593 return (t - eps) / (1 + eps);

1594 }

1595

1596

1597 void C2f(double eps, double c[]) {

1598 static const double coeff[] = {

1599

1600 1, 2, 16, 32,

1601

1602 35, 64, 384, 2048,

1603

1604 15, 80, 768,

1605

1606 7, 35, 512,

1607

1608 63, 1280,

1609

1610 77, 2048,

1611 };

1612 double

1613 eps2 = sq(eps),

1614 d = eps;

1615 int o = 0, l;

1616 for (l = 1; l <= nC2; ++l) {

1617 int m = (nC2 - l) / 2;

1618 c[l] = d * polyvalx(m, coeff + o, eps2) / coeff[o + m + 1];

1619 o += m + 2;

1620 d *= eps;

1621 }

1622 }

1623

1624

1626 static const double coeff[] = {

1627

1628 -3, 128,

1629

1630 -2, -3, 64,

1631

1632 -1, -3, -1, 16,

1633

1634 3, -1, -2, 8,

1635

1636 1, -1, 2,

1637

1638 1, 1,

1639 };

1640 int o = 0, k = 0, j;

1641 for (j = nA3 - 1; j >= 0; --j) {

1642 int m = nA3 - j - 1 < j ? nA3 - j - 1 : j;

1643 g->A3x[k++] = polyvalx(m, coeff + o, g->n) / coeff[o + m + 1];

1644 o += m + 2;

1645 }

1646 }

1647

1648

1650 static const double coeff[] = {

1651

1652 3, 128,

1653

1654 2, 5, 128,

1655

1656 -1, 3, 3, 64,

1657

1658 -1, 0, 1, 8,

1659

1660 -1, 1, 4,

1661

1662 5, 256,

1663

1664 1, 3, 128,

1665

1666 -3, -2, 3, 64,

1667

1668 1, -3, 2, 32,

1669

1670 7, 512,

1671

1672 -10, 9, 384,

1673

1674 5, -9, 5, 192,

1675

1676 7, 512,

1677

1678 -14, 7, 512,

1679

1680 21, 2560,

1681 };

1682 int o = 0, k = 0, l, j;

1683 for (l = 1; l < nC3; ++l) {

1684 for (j = nC3 - 1; j >= l; --j) {

1685 int m = nC3 - j - 1 < j ? nC3 - j - 1 : j;

1686 g->C3x[k++] = polyvalx(m, coeff + o, g->n) / coeff[o + m + 1];

1687 o += m + 2;

1688 }

1689 }

1690 }

1691

1692

1694 static const double coeff[] = {

1695

1696 97, 15015,

1697

1698 1088, 156, 45045,

1699

1700 -224, -4784, 1573, 45045,

1701

1702 -10656, 14144, -4576, -858, 45045,

1703

1704 64, 624, -4576, 6864, -3003, 15015,

1705

1706 100, 208, 572, 3432, -12012, 30030, 45045,

1707

1708 1, 9009,

1709

1710 -2944, 468, 135135,

1711

1712 5792, 1040, -1287, 135135,

1713

1714 5952, -11648, 9152, -2574, 135135,

1715

1716 -64, -624, 4576, -6864, 3003, 135135,

1717

1718 8, 10725,

1719

1720 1856, -936, 225225,

1721

1722 -8448, 4992, -1144, 225225,

1723

1724 -1440, 4160, -4576, 1716, 225225,

1725

1726 -136, 63063,

1727

1728 1024, -208, 105105,

1729

1730 3584, -3328, 1144, 315315,

1731

1732 -128, 135135,

1733

1734 -2560, 832, 405405,

1735

1736 128, 99099,

1737 };

1738 int o = 0, k = 0, l, j;

1739 for (l = 0; l < nC4; ++l) {

1740 for (j = nC4 - 1; j >= l; --j) {

1741 int m = nC4 - j - 1;

1742 g->C4x[k++] = polyvalx(m, coeff + o, g->n) / coeff[o + m + 1];

1743 o += m + 2;

1744 }

1745 }

1746 }

1747

1748 int transit(double lon1, double lon2) {

1749 double lon12;

1750

1751

1752

1753 lon12 = AngDiff(lon1, lon2, nullptr);

1755 lon2 = AngNormalize(lon2);

1756 return

1757 lon12 > 0 && ((lon1 < 0 && lon2 >= 0) ||

1758 (lon1 > 0 && lon2 == 0)) ? 1 :

1759 (lon12 < 0 && lon1 >= 0 && lon2 < 0 ? -1 : 0);

1760 }

1761

1762 int transitdirect(double lon1, double lon2) {

1763

1764

1765 lon1 = remainder(lon1, 2.0 * td); lon2 = remainder(lon2, 2.0 * td);

1766 return ( (lon2 >= 0 && lon2 < td ? 0 : 1) -

1767 (lon1 >= 0 && lon1 < td ? 0 : 1) );

1768 }

1769

1770 void accini(double s[]) {

1771

1772 s[0] = s[1] = 0;

1773 }

1774

1775 void acccopy(const double s[], double t[]) {

1776

1777 t[0] = s[0]; t[1] = s[1];

1778 }

1779

1780 void accadd(double s[], double y) {

1781

1782 double u, z = sumx(y, s[1], &u);

1783 s[0] = sumx(z, s[0], &s[1]);

1784 if (s[0] == 0)

1785 s[0] = u;

1786 else

1787 s[1] = s[1] + u;

1788 }

1789

1790 double accsum(const double s[], double y) {

1791

1792 double t[2];

1793 acccopy(s, t);

1794 accadd(t, y);

1795 return t[0];

1796 }

1797

1798 void accneg(double s[]) {

1799

1800 s[0] = -s[0]; s[1] = -s[1];

1801 }

1802

1803 void accrem(double s[], double y) {

1804

1805 s[0] = remainder(s[0], y);

1806 accadd(s, 0.0);

1807 }

1808

1810 p->polyline = (polylinep != 0);

1812 }

1813

1815 p->lat0 = p->lon0 = p->lat = p->lon = NaN;

1816 accini(p->P);

1817 accini(p->A);

1818 p->num = p->crossings = 0;

1819 }

1820

1823 double lat, double lon) {

1824 if (p->num == 0) {

1825 p->lat0 = p->lat = lat;

1826 p->lon0 = p->lon = lon;

1827 } else {

1828 double s12, S12 = 0;

1830 &s12, nullptr, nullptr, nullptr, nullptr, nullptr,

1831 p->polyline ? nullptr : &S12);

1832 accadd(p->P, s12);

1833 if (!p->polyline) {

1834 accadd(p->A, S12);

1835 p->crossings += transit(p->lon, lon);

1836 }

1837 p->lat = lat; p->lon = lon;

1838 }

1839 ++p->num;

1840 }

1841

1844 double azi, double s) {

1845 if (p->num) {

1846

1847

1848 double lat = 0, lon = 0, S12 = 0;

1850 &lat, &lon, nullptr,

1851 nullptr, nullptr, nullptr, nullptr,

1852 p->polyline ? nullptr : &S12);

1853 accadd(p->P, s);

1854 if (!p->polyline) {

1855 accadd(p->A, S12);

1856 p->crossings += transitdirect(p->lon, lon);

1857 }

1858 p->lat = lat; p->lon = lon;

1859 ++p->num;

1860 }

1861 }

1862

1865 boolx reverse, boolx sign,

1866 double* pA, double* pP) {

1867 double s12, S12, t[2];

1868 if (p->num < 2) {

1869 if (pP) *pP = 0;

1870 if (!p->polyline && pA) *pA = 0;

1871 return p->num;

1872 }

1873 if (p->polyline) {

1874 if (pP) *pP = p->P[0];

1875 return p->num;

1876 }

1878 &s12, nullptr, nullptr, nullptr, nullptr, nullptr, &S12);

1879 if (pP) *pP = accsum(p->P, s12);

1880 acccopy(p->A, t);

1881 accadd(t, S12);

1882 if (pA) *pA = areareduceA(t, 4 * pi * g->c2,

1883 p->crossings + transit(p->lon, p->lon0),

1884 reverse, sign);

1885 return p->num;

1886 }

1887

1890 double lat, double lon,

1891 boolx reverse, boolx sign,

1892 double* pA, double* pP) {

1893 double perimeter, tempsum;

1894 int crossings, i;

1895 unsigned num = p->num + 1;

1896 if (num == 1) {

1897 if (pP) *pP = 0;

1898 if (!p->polyline && pA) *pA = 0;

1899 return num;

1900 }

1901 perimeter = p->P[0];

1902 tempsum = p->polyline ? 0 : p->A[0];

1903 crossings = p->crossings;

1904 for (i = 0; i < (p->polyline ? 1 : 2); ++i) {

1905 double s12, S12 = 0;

1907 i == 0 ? p->lat : lat, i == 0 ? p->lon : lon,

1908 i != 0 ? p->lat0 : lat, i != 0 ? p->lon0 : lon,

1909 &s12, nullptr, nullptr, nullptr, nullptr, nullptr,

1910 p->polyline ? nullptr : &S12);

1911 perimeter += s12;

1912 if (!p->polyline) {

1913 tempsum += S12;

1914 crossings += transit(i == 0 ? p->lon : lon,

1915 i != 0 ? p->lon0 : lon);

1916 }

1917 }

1918

1919 if (pP) *pP = perimeter;

1920 if (p->polyline)

1921 return num;

1922

1923 if (pA) *pA = areareduceB(tempsum, 4 * pi * g->c2, crossings, reverse, sign);

1924 return num;

1925 }

1926

1929 double azi, double s,

1930 boolx reverse, boolx sign,

1931 double* pA, double* pP) {

1932 double perimeter, tempsum;

1933 int crossings;

1934 unsigned num = p->num + 1;

1935 if (num == 1) {

1936 if (pP) *pP = NaN;

1937 if (!p->polyline && pA) *pA = NaN;

1938 return 0;

1939 }

1940 perimeter = p->P[0] + s;

1941 if (p->polyline) {

1942 if (pP) *pP = perimeter;

1943 return num;

1944 }

1945

1946 tempsum = p->A[0];

1947 crossings = p->crossings;

1948 {

1949

1950

1951 double lat = 0, lon = 0, s12, S12 = 0;

1953 &lat, &lon, nullptr,

1954 nullptr, nullptr, nullptr, nullptr, &S12);

1955 tempsum += S12;

1956 crossings += transitdirect(p->lon, lon);

1958 &s12, nullptr, nullptr, nullptr, nullptr, nullptr, &S12);

1959 perimeter += s12;

1960 tempsum += S12;

1961 crossings += transit(lon, p->lon0);

1962 }

1963

1964 if (pP) *pP = perimeter;

1965 if (pA) *pA = areareduceB(tempsum, 4 * pi * g->c2, crossings, reverse, sign);

1966 return num;

1967 }

1968

1970 double lats[], double lons[], int n,

1971 double* pA, double* pP) {

1972 int i;

1975 for (i = 0; i < n; ++i)

1978 }

1979

1980 double areareduceA(double area[], double area0,

1981 int crossings, boolx reverse, boolx sign) {

1982 accrem(area, area0);

1983 if (crossings & 1)

1984 accadd(area, (area[0] < 0 ? 1 : -1) * area0/2);

1985

1986

1987 if (!reverse)

1988 accneg(area);

1989

1990 if (sign) {

1991 if (area[0] > area0/2)

1992 accadd(area, -area0);

1993 else if (area[0] <= -area0/2)

1994 accadd(area, +area0);

1995 } else {

1996 if (area[0] >= area0)

1997 accadd(area, -area0);

1998 else if (area[0] < 0)

1999 accadd(area, +area0);

2000 }

2001 return 0 + area[0];

2002 }

2003

2004 double areareduceB(double area, double area0,

2005 int crossings, boolx reverse, boolx sign) {

2006 area = remainder(area, area0);

2007 if (crossings & 1)

2008 area += (area < 0 ? 1 : -1) * area0/2;

2009

2010

2011 if (!reverse)

2012 area *= -1;

2013

2014 if (sign) {

2015 if (area > area0/2)

2016 area -= area0;

2017 else if (area <= -area0/2)

2018 area += area0;

2019 } else {

2020 if (area >= area0)

2021 area -= area0;

2022 else if (area < 0)

2023 area += area0;

2024 }

2025 return 0 + area;

2026 }

2027

2028

API for the geodesic routines in C.

double GEOD_DLL geod_gendirect(const struct geod_geodesic *g, double lat1, double lon1, double azi1, unsigned flags, double s12_a12, double *plat2, double *plon2, double *pazi2, double *ps12, double *pm12, double *pM12, double *pM21, double *pS12)

void GEOD_DLL geod_polygon_addedge(const struct geod_geodesic *g, struct geod_polygon *p, double azi, double s)

void GEOD_DLL geod_polygon_addpoint(const struct geod_geodesic *g, struct geod_polygon *p, double lat, double lon)

void GEOD_DLL geod_setdistance(struct geod_geodesicline *l, double s13)

void GEOD_DLL geod_gensetdistance(struct geod_geodesicline *l, unsigned flags, double s13_a13)

unsigned GEOD_DLL geod_polygon_testpoint(const struct geod_geodesic *g, const struct geod_polygon *p, double lat, double lon, int reverse, int sign, double *pA, double *pP)

void GEOD_DLL geod_polygon_clear(struct geod_polygon *p)

void GEOD_DLL geod_init(struct geod_geodesic *g, double a, double f)

void GEOD_DLL geod_directline(struct geod_geodesicline *l, const struct geod_geodesic *g, double lat1, double lon1, double azi1, double s12, unsigned caps)

void GEOD_DLL geod_inverse(const struct geod_geodesic *g, double lat1, double lon1, double lat2, double lon2, double *ps12, double *pazi1, double *pazi2)

unsigned GEOD_DLL geod_polygon_testedge(const struct geod_geodesic *g, const struct geod_polygon *p, double azi, double s, int reverse, int sign, double *pA, double *pP)

void GEOD_DLL geod_lineinit(struct geod_geodesicline *l, const struct geod_geodesic *g, double lat1, double lon1, double azi1, unsigned caps)

unsigned GEOD_DLL geod_polygon_compute(const struct geod_geodesic *g, const struct geod_polygon *p, int reverse, int sign, double *pA, double *pP)

void GEOD_DLL geod_position(const struct geod_geodesicline *l, double s12, double *plat2, double *plon2, double *pazi2)

double GEOD_DLL geod_geninverse(const struct geod_geodesic *g, double lat1, double lon1, double lat2, double lon2, double *ps12, double *pazi1, double *pazi2, double *pm12, double *pM12, double *pM21, double *pS12)

double GEOD_DLL geod_genposition(const struct geod_geodesicline *l, unsigned flags, double s12_a12, double *plat2, double *plon2, double *pazi2, double *ps12, double *pm12, double *pM12, double *pM21, double *pS12)

void GEOD_DLL geod_direct(const struct geod_geodesic *g, double lat1, double lon1, double azi1, double s12, double *plat2, double *plon2, double *pazi2)

void GEOD_DLL geod_polygonarea(const struct geod_geodesic *g, double lats[], double lons[], int n, double *pA, double *pP)

void GEOD_DLL geod_inverseline(struct geod_geodesicline *l, const struct geod_geodesic *g, double lat1, double lon1, double lat2, double lon2, unsigned caps)

void GEOD_DLL geod_gendirectline(struct geod_geodesicline *l, const struct geod_geodesic *g, double lat1, double lon1, double azi1, unsigned flags, double s12_a12, unsigned caps)

void GEOD_DLL geod_polygon_init(struct geod_polygon *p, int polylinep)