GeographicLib: TransverseMercator.cpp Source File (original) (raw)

45 bool exact, bool extendp)

46 : _a(a)

47 , _f(f)

48 , _k0(k0)

49 , _exact(exact)

50 , _e2(_f * (2 - _f))

51 , _es((_f < 0 ? -1 : 1) * sqrt(fabs(_e2)))

52 , _e2m(1 - _e2)

55 , _c( sqrt(_e2m) * exp(Math::eatanhe(real(1), _es)) )

56 , _n(_f / (2 - _f))

59 {

60 if (_exact) return;

61 if (!(isfinite(_a) && _a > 0))

62 throw GeographicErr("Equatorial radius is not positive");

63 if (!(isfinite(_f) && _f < 1))

64 throw GeographicErr("Polar semi-axis is not positive");

65 if (!(isfinite(_k0) && _k0 > 0))

67 if (extendp)

68 throw GeographicErr("TransverseMercator extendp not allowed if !exact");

71#if GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 2

72 static const real b1coeff[] = {

74 1, 16, 64, 64,

75 };

76#elif GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 3

77 static const real b1coeff[] = {

79 1, 4, 64, 256, 256,

80 };

81#elif GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 4

82 static const real b1coeff[] = {

84 25, 64, 256, 4096, 16384, 16384,

85 };

86#else

87#error "Bad value for GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER"

88#endif

90#if GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 4

91 static const real alpcoeff[] = {

93 164, 225, -480, 360, 720,

95 557, -864, 390, 1440,

97 -1236, 427, 1680,

99 49561, 161280,

100 };

101#elif GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 5

102 static const real alpcoeff[] = {

103

104 -635, 328, 450, -960, 720, 1440,

105

106 4496, 3899, -6048, 2730, 10080,

107

108 15061, -19776, 6832, 26880,

109

110 -171840, 49561, 161280,

111

112 34729, 80640,

113 };

114#elif GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 6

115 static const real alpcoeff[] = {

116

117 31564, -66675, 34440, 47250, -100800, 75600, 151200,

118

119 -1983433, 863232, 748608, -1161216, 524160, 1935360,

120

121 670412, 406647, -533952, 184464, 725760,

122

123 6601661, -7732800, 2230245, 7257600,

124

125 -13675556, 3438171, 7983360,

126

127 212378941, 319334400,

128 };

129#elif GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 7

130 static const real alpcoeff[] = {

131

132 1804025, 2020096, -4267200, 2204160, 3024000, -6451200, 4838400, 9676800,

133

134 4626384, -9917165, 4316160, 3743040, -5806080, 2620800, 9676800,

135

136 -67102379, 26816480, 16265880, -21358080, 7378560, 29030400,

137

138 155912000, 72618271, -85060800, 24532695, 79833600,

139

140 102508609, -109404448, 27505368, 63866880,

141

142 -12282192400LL, 2760926233LL, 4151347200LL,

143

144 1522256789, 1383782400,

145 };

146#elif GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 8

147 static const real alpcoeff[] = {

148

149 -75900428, 37884525, 42422016, -89611200, 46287360, 63504000, -135475200,

150 101606400, 203212800,

151

152 148003883, 83274912, -178508970, 77690880, 67374720, -104509440,

153 47174400, 174182400,

154

155 318729724, -738126169, 294981280, 178924680, -234938880, 81164160,

156 319334400,

157

158 -40176129013LL, 14967552000LL, 6971354016LL, -8165836800LL, 2355138720LL,

159 7664025600LL,

160

161 10421654396LL, 3997835751LL, -4266773472LL, 1072709352, 2490808320LL,

162

163 175214326799LL, -171950693600LL, 38652967262LL, 58118860800LL,

164

165 -67039739596LL, 13700311101LL, 12454041600LL,

166

167 1424729850961LL, 743921418240LL,

168 };

169#else

170#error "Bad value for GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER"

171#endif

172

173#if GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 4

174 static const real betcoeff[] = {

175

176 -4, 555, -960, 720, 1440,

177

178 -437, 96, 30, 1440,

179

180 -148, 119, 3360,

181

182 4397, 161280,

183 };

184#elif GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 5

185 static const real betcoeff[] = {

186

187 -3645, -64, 8880, -15360, 11520, 23040,

188

189 4416, -3059, 672, 210, 10080,

190

191 -627, -592, 476, 13440,

192

193 -3520, 4397, 161280,

194

195 4583, 161280,

196 };

197#elif GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 6

198 static const real betcoeff[] = {

199

200 384796, -382725, -6720, 932400, -1612800, 1209600, 2419200,

201

202 -1118711, 1695744, -1174656, 258048, 80640, 3870720,

203

204 22276, -16929, -15984, 12852, 362880,

205

206 -830251, -158400, 197865, 7257600,

207

208 -435388, 453717, 15966720,

209

210 20648693, 638668800,

211 };

212#elif GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 7

213 static const real betcoeff[] = {

214

215 -5406467, 6156736, -6123600, -107520, 14918400, -25804800, 19353600,

216 38707200,

217

218 829456, -5593555, 8478720, -5873280, 1290240, 403200, 19353600,

219

220 9261899, 3564160, -2708640, -2557440, 2056320, 58060800,

221

222 14928352, -9132761, -1742400, 2176515, 79833600,

223

224 -8005831, -1741552, 1814868, 63866880,

225

226 -261810608, 268433009, 8302694400LL,

227

228 219941297, 5535129600LL,

229 };

230#elif GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 8

231 static const real betcoeff[] = {

232

233 31777436, -37845269, 43097152, -42865200, -752640, 104428800, -180633600,

234 135475200, 270950400,

235

236 24749483, 14930208, -100683990, 152616960, -105719040, 23224320, 7257600,

237 348364800,

238

239 -232468668, 101880889, 39205760, -29795040, -28131840, 22619520,

240 638668800,

241

242 324154477, 1433121792, -876745056, -167270400, 208945440, 7664025600LL,

243

244 457888660, -312227409, -67920528, 70779852, 2490808320LL,

245

246 -19841813847LL, -3665348512LL, 3758062126LL, 116237721600LL,

247

248 -1989295244, 1979471673, 49816166400LL,

249

250 191773887257LL, 3719607091200LL,

251 };

252#else

253#error "Bad value for GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER"

254#endif

255

256 static_assert(sizeof(b1coeff) / sizeof(real) == maxpow_/2 + 2,

257 "Coefficient array size mismatch for b1");

258 static_assert(sizeof(alpcoeff) / sizeof(real) ==

259 (maxpow_ * (maxpow_ + 3))/2,

260 "Coefficient array size mismatch for alp");

261 static_assert(sizeof(betcoeff) / sizeof(real) ==

262 (maxpow_ * (maxpow_ + 3))/2,

263 "Coefficient array size mismatch for bet");

264 int m = maxpow_/2;

266

267

268 _a1 = _b1 * _a;

269 int o = 0;

270 real d = _n;

271 for (int l = 1; l <= maxpow_; ++l) {

272 m = maxpow_ - l;

273 _alp[l] = d * Math::polyval(m, alpcoeff + o, _n) / alpcoeff[o + m + 1];

274 _bet[l] = d * Math::polyval(m, betcoeff + o, _n) / betcoeff[o + m + 1];

275 o += m + 2;

276 d *= _n;

277 }

278

279

280 }

356 real& x, real& y,

357 real& gamma, real& k) const {

358 if (_exact)

359 return _tmexact.Forward(lon0, lat, lon, x, y, gamma, k);

362

363 int

364 latsign = signbit(lat) ? -1 : 1,

365 lonsign = signbit(lon) ? -1 : 1;

366 lon *= lonsign;

367 lat *= latsign;

368 bool backside = lon > Math::qd;

369 if (backside) {

370 if (lat == 0)

371 latsign = -1;

373 }

374 real sphi, cphi, slam, clam;

377

378

379

380

381

382

383

384

385

386

387

388

389

390

391

392

393

394 real etap, xip;

396 real

397 tau = sphi / cphi,

399 xip = atan2(taup, clam);

400

401

402 etap = asinh(slam / hypot(taup, clam));

403

404

405

406

407 gamma = Math::atan2d(slam * taup, clam * hypot(real(1), taup));

408

409

410

411

412

413

414

415 k = sqrt(_e2m + _e2 * Math::sq(cphi)) * hypot(real(1), tau)

416 / hypot(taup, clam);

417 } else {

419 etap = 0;

420 gamma = lon;

421 k = _c;

422 }

423

424

425

426

427

428

429

430

431

432

433

434

435

436

437

438

439

440

441

442

443

444

445

446

447

448

449

450

451

452

453

454

455

456

457

458

459

460

461

462

463

464

465

466

467

468

469

470

471

472

473

474

475

476

477

478

479

480

481

482

483

484

485

486

487

488 real

489 c0 = cos(2 * xip), ch0 = cosh(2 * etap),

490 s0 = sin(2 * xip), sh0 = sinh(2 * etap);

491 complex a(2 * c0 * ch0, -2 * s0 * sh0);

492 int n = maxpow_;

493 complex

494 y0(n & 1 ? _alp[n] : 0), y1,

495 z0(n & 1 ? 2*n * _alp[n] : 0), z1;

496 if (n & 1) --n;

497 while (n) {

498 y1 = a * y0 - y1 + _alp[n];

499 z1 = a * z0 - z1 + 2*n * _alp[n];

500 --n;

501 y0 = a * y1 - y0 + _alp[n];

502 z0 = a * z1 - z0 + 2*n * _alp[n];

503 --n;

504 }

505 a /= real(2);

506 z1 = real(1) - z1 + a * z0;

507 a = complex(s0 * ch0, c0 * sh0);

508 y1 = complex(xip, etap) + a * y0;

509

510

512 k *= _b1 * abs(z1);

513 real xi = y1.real(), eta = y1.imag();

514 y = _a1 * _k0 * (backside ? Math::pi() - xi : xi) * latsign;

515 x = _a1 * _k0 * eta * lonsign;

516 if (backside)

518 gamma *= latsign * lonsign;

520 k *= _k0;

521 }

524 real& lat, real& lon,

525 real& gamma, real& k) const {

526 if (_exact)

527 return _tmexact.Reverse(lon0, x, y, lat, lon, gamma, k);

528

529

530

531 real

532 xi = y / (_a1 * _k0),

533 eta = x / (_a1 * _k0);

534

535 int

536 xisign = signbit(xi) ? -1 : 1,

537 etasign = signbit(eta) ? -1 : 1;

538 xi *= xisign;

539 eta *= etasign;

540 bool backside = xi > Math::pi()/2;

541 if (backside)

543 real

544 c0 = cos(2 * xi), ch0 = cosh(2 * eta),

545 s0 = sin(2 * xi), sh0 = sinh(2 * eta);

546 complex a(2 * c0 * ch0, -2 * s0 * sh0);

547 int n = maxpow_;

548 complex

549 y0(n & 1 ? -_bet[n] : 0), y1,

550 z0(n & 1 ? -2*n * _bet[n] : 0), z1;

551 if (n & 1) --n;

552 while (n) {

553 y1 = a * y0 - y1 - _bet[n];

554 z1 = a * z0 - z1 - 2*n * _bet[n];

555 --n;

556 y0 = a * y1 - y0 - _bet[n];

557 z0 = a * z1 - z0 - 2*n * _bet[n];

558 --n;

559 }

560 a /= real(2);

561 z1 = real(1) - z1 + a * z0;

562 a = complex(s0 * ch0, c0 * sh0);

563 y1 = complex(xi, eta) + a * y0;

564

566 k = _b1 / abs(z1);

567

568

569

570

571

572 real

573 xip = y1.real(), etap = y1.imag(),

574 s = sinh(etap),

575 c = fmax(real(0), cos(xip)),

576 r = hypot(s, c);

577 if (r != 0) {

578 lon = Math::atan2d(s, c);

579

580 real

581 sxip = sin(xip),

583 gamma += Math::atan2d(sxip * tanh(etap), c);

585

586 k *= sqrt(_e2m + _e2 / (1 + Math::sq(tau))) *

587 hypot(real(1), tau) * r;

588 } else {

590 lon = 0;

591 k *= _c;

592 }

593 lat *= xisign;

594 if (backside)

596 lon *= etasign;

598 if (backside)

600 gamma *= xisign * etasign;

602 k *= _k0;

603 }