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

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

136

137#if defined(_MSC_VER)

138

139# pragma warning (disable: 4701)

140#endif

141

143

144 using namespace std;

145

146 vectorMath::real& SphericalEngine::sqrttable() {

147 static vector sqrttable(0);

148 return sqrttable;

149 }

150

151 template<bool gradp, SphericalEngine::normalization norm, int L>

153 real x, real y, real z, real a,

154 real& gradx, real& grady, real& gradz)

155 {

156 static_assert(L > 0, "L must be positive");

157 static_assert(norm == FULL || norm == SCHMIDT, "Unknown normalization");

158 int N = c[0].nmx(), M = c[0].mmx();

159

160 real

161 p = hypot(x, y),

162 cl = p != 0 ? x / p : 1,

163 sl = p != 0 ? y / p : 0,

164 r = hypot(z, p),

165 t = r != 0 ? z / r : 0,

166 u = r != 0 ? fmax(p / r, eps()) : 1,

167 q = a / r;

168 real

170 uq = u * q,

172 tu = t / u;

173

174 real vc = 0, vc2 = 0, vs = 0, vs2 = 0;

175

176

177 real vrc = 0, vrc2 = 0, vrs = 0, vrs2 = 0;

178 real vtc = 0, vtc2 = 0, vts = 0, vts2 = 0;

179 real vlc = 0, vlc2 = 0, vls = 0, vls2 = 0;

180 int k[L];

181 const vector& root( sqrttable() );

182 for (int m = M; m >= 0; --m) {

183

184 real

185 wc = 0, wc2 = 0, ws = 0, ws2 = 0,

186 wrc = 0, wrc2 = 0, wrs = 0, wrs2 = 0,

187 wtc = 0, wtc2 = 0, wts = 0, wts2 = 0;

188 for (int l = 0; l < L; ++l)

189 k[l] = c[l].index(N, m) + 1;

190 for (int n = N; n >= m; --n) {

191 real w, A, Ax, B, R;

192 switch (norm) {

194 w = root[2 * n + 1] / (root[n - m + 1] * root[n + m + 1]);

195 Ax = q * w * root[2 * n + 3];

196 A = t * Ax;

197 B = - q2 * root[2 * n + 5] /

198 (w * root[n - m + 2] * root[n + m + 2]);

199 break;

201 w = root[n - m + 1] * root[n + m + 1];

202 Ax = q * (2 * n + 1) / w;

203 A = t * Ax;

204 B = - q2 * w / (root[n - m + 2] * root[n + m + 2]);

205 break;

206 default: break;

207 }

208 R = c[0].Cv(--k[0]);

209 for (int l = 1; l < L; ++l)

210 R += c[l].Cv(--k[l], n, m, f[l]);

211 R *= scale();

212 w = A * wc + B * wc2 + R; wc2 = wc; wc = w;

213 if (gradp) {

214 w = A * wrc + B * wrc2 + (n + 1) * R; wrc2 = wrc; wrc = w;

215 w = A * wtc + B * wtc2 - u*Ax * wc2; wtc2 = wtc; wtc = w;

216 }

217 if (m) {

218 R = c[0].Sv(k[0]);

219 for (int l = 1; l < L; ++l)

220 R += c[l].Sv(k[l], n, m, f[l]);

221 R *= scale();

222 w = A * ws + B * ws2 + R; ws2 = ws; ws = w;

223 if (gradp) {

224 w = A * wrs + B * wrs2 + (n + 1) * R; wrs2 = wrs; wrs = w;

225 w = A * wts + B * wts2 - u*Ax * ws2; wts2 = wts; wts = w;

226 }

227 }

228 }

229

230

231 if (m) {

232 real v, A, B;

233 switch (norm) {

235 v = root[2] * root[2 * m + 3] / root[m + 1];

236 A = cl * v * uq;

237 B = - v * root[2 * m + 5] / (root[8] * root[m + 2]) * uq2;

238 break;

240 v = root[2] * root[2 * m + 1] / root[m + 1];

241 A = cl * v * uq;

242 B = - v * root[2 * m + 3] / (root[8] * root[m + 2]) * uq2;

243 break;

244 default: break;

245 }

246 v = A * vc + B * vc2 + wc ; vc2 = vc ; vc = v;

247 v = A * vs + B * vs2 + ws ; vs2 = vs ; vs = v;

248 if (gradp) {

249

250 wtc += m * tu * wc; wts += m * tu * ws;

251 v = A * vrc + B * vrc2 + wrc; vrc2 = vrc; vrc = v;

252 v = A * vrs + B * vrs2 + wrs; vrs2 = vrs; vrs = v;

253 v = A * vtc + B * vtc2 + wtc; vtc2 = vtc; vtc = v;

254 v = A * vts + B * vts2 + wts; vts2 = vts; vts = v;

255 v = A * vlc + B * vlc2 + m*ws; vlc2 = vlc; vlc = v;

256 v = A * vls + B * vls2 - m*wc; vls2 = vls; vls = v;

257 }

258 } else {

259 real A, B, qs;

260 switch (norm) {

262 A = root[3] * uq;

263 B = - root[15]/2 * uq2;

264 break;

266 A = uq;

267 B = - root[3]/2 * uq2;

268 break;

269 default: break;

270 }

271 qs = q / scale();

272 vc = qs * (wc + A * (cl * vc + sl * vs ) + B * vc2);

273 if (gradp) {

274 qs /= r;

275

276

277

278

279 vrc = - qs * (wrc + A * (cl * vrc + sl * vrs) + B * vrc2);

280 vtc = qs * (wtc + A * (cl * vtc + sl * vts) + B * vtc2);

281 vlc = qs / u * ( A * (cl * vlc + sl * vls) + B * vlc2);

282 }

283 }

284 }

285

286 if (gradp) {

287

288 gradx = cl * (u * vrc + t * vtc) - sl * vlc;

289 grady = sl * (u * vrc + t * vtc) + cl * vlc;

290 gradz = t * vrc - u * vtc ;

291 }

292 return vc;

293 }

294

295 template<bool gradp, SphericalEngine::normalization norm, int L>

297 real p, real z, real a) {

298

299 static_assert(L > 0, "L must be positive");

300 static_assert(norm == FULL || norm == SCHMIDT, "Unknown normalization");

301 int N = c[0].nmx(), M = c[0].mmx();

302

303 real

304 r = hypot(z, p),

305 t = r != 0 ? z / r : 0,

306 u = r != 0 ? fmax(p / r, eps()) : 1,

307 q = a / r;

308 real

310 tu = t / u;

312 int k[L];

313 const vector& root( sqrttable() );

314 for (int m = M; m >= 0; --m) {

315

316 real

317 wc = 0, wc2 = 0, ws = 0, ws2 = 0,

318 wrc = 0, wrc2 = 0, wrs = 0, wrs2 = 0,

319 wtc = 0, wtc2 = 0, wts = 0, wts2 = 0;

320 for (int l = 0; l < L; ++l)

321 k[l] = c[l].index(N, m) + 1;

322 for (int n = N; n >= m; --n) {

323 real w, A, Ax, B, R;

324 switch (norm) {

326 w = root[2 * n + 1] / (root[n - m + 1] * root[n + m + 1]);

327 Ax = q * w * root[2 * n + 3];

328 A = t * Ax;

329 B = - q2 * root[2 * n + 5] /

330 (w * root[n - m + 2] * root[n + m + 2]);

331 break;

333 w = root[n - m + 1] * root[n + m + 1];

334 Ax = q * (2 * n + 1) / w;

335 A = t * Ax;

336 B = - q2 * w / (root[n - m + 2] * root[n + m + 2]);

337 break;

338 default: break;

339 }

340 R = c[0].Cv(--k[0]);

341 for (int l = 1; l < L; ++l)

342 R += c[l].Cv(--k[l], n, m, f[l]);

343 R *= scale();

344 w = A * wc + B * wc2 + R; wc2 = wc; wc = w;

345 if (gradp) {

346 w = A * wrc + B * wrc2 + (n + 1) * R; wrc2 = wrc; wrc = w;

347 w = A * wtc + B * wtc2 - u*Ax * wc2; wtc2 = wtc; wtc = w;

348 }

349 if (m) {

350 R = c[0].Sv(k[0]);

351 for (int l = 1; l < L; ++l)

352 R += c[l].Sv(k[l], n, m, f[l]);

353 R *= scale();

354 w = A * ws + B * ws2 + R; ws2 = ws; ws = w;

355 if (gradp) {

356 w = A * wrs + B * wrs2 + (n + 1) * R; wrs2 = wrs; wrs = w;

357 w = A * wts + B * wts2 - u*Ax * ws2; wts2 = wts; wts = w;

358 }

359 }

360 }

361 if (!gradp)

362 circ.SetCoeff(m, wc, ws);

363 else {

364

365 wtc += m * tu * wc; wts += m * tu * ws;

366 circ.SetCoeff(m, wc, ws, wrc, wrs, wtc, wts);

367 }

368 }

369

370 return circ;

371 }

372

374

375 vector& root( sqrttable() );

376 int L = max(2 * N + 5, 15) + 1, oldL = int(root.size());

377 if (oldL >= L)

378 return;

379 root.resize(L);

380 for (int l = oldL; l < L; ++l)

381 root[l] = sqrt(real(l));

382 }

383

385 vector& C,

386 vector& S,

387 bool truncate) {

388 if (truncate) {

389 if (!((N >= M && M >= 0) || (N == -1 && M == -1)))

390

391 throw GeographicErr("Bad requested degree and order " +

393 }

394 int nm[2];

395 Utility::readarray<int, int, false>(stream, nm, 2);

396 int N0 = nm[0], M0 = nm[1];

397 if (!((N0 >= M0 && M0 >= 0) || (N0 == -1 && M0 == -1)))

398

401 N = truncate ? min(N, N0) : N0;

402 M = truncate ? min(M, M0) : M0;

407 if (N == N0) {

408 Utility::readarray<double, real, false>(stream, C);

409 if (skip) stream.seekg(streamoff(skip), ios::cur);

410 Utility::readarray<double, real, false>(stream, S);

411 if (skip) stream.seekg(streamoff(skip), ios::cur);

412 } else {

413 for (int m = 0, k = 0; m <= M; ++m) {

414 Utility::readarray<double, real, false>(stream, &C[k], N + 1 - m);

415 stream.seekg((N0 - N) * sizeof(double), ios::cur);

416 k += N + 1 - m;

417 }

418 if (skip) stream.seekg(streamoff(skip), ios::cur);

419 for (int m = 1, k = 0; m <= M; ++m) {

420 Utility::readarray<double, real, false>(stream, &S[k], N + 1 - m);

421 stream.seekg((N0 - N) * sizeof(double), ios::cur);

422 k += N + 1 - m;

423 }

424 if (skip) stream.seekg(streamoff(skip), ios::cur);

425 }

426 return;

427 }

428

429

431 SphericalEngine::Value<true, SphericalEngine::FULL, 1>

432 (const coeff[], const real[], real, real, real, real, real&, real&, real&);

434 SphericalEngine::Value<false, SphericalEngine::FULL, 1>

435 (const coeff[], const real[], real, real, real, real, real&, real&, real&);

437 SphericalEngine::Value<true, SphericalEngine::SCHMIDT, 1>

438 (const coeff[], const real[], real, real, real, real, real&, real&, real&);

440 SphericalEngine::Value<false, SphericalEngine::SCHMIDT, 1>

441 (const coeff[], const real[], real, real, real, real, real&, real&, real&);

442

444 SphericalEngine::Value<true, SphericalEngine::FULL, 2>

445 (const coeff[], const real[], real, real, real, real, real&, real&, real&);

447 SphericalEngine::Value<false, SphericalEngine::FULL, 2>

448 (const coeff[], const real[], real, real, real, real, real&, real&, real&);

450 SphericalEngine::Value<true, SphericalEngine::SCHMIDT, 2>

451 (const coeff[], const real[], real, real, real, real, real&, real&, real&);

453 SphericalEngine::Value<false, SphericalEngine::SCHMIDT, 2>

454 (const coeff[], const real[], real, real, real, real, real&, real&, real&);

455

457 SphericalEngine::Value<true, SphericalEngine::FULL, 3>

458 (const coeff[], const real[], real, real, real, real, real&, real&, real&);

460 SphericalEngine::Value<false, SphericalEngine::FULL, 3>

461 (const coeff[], const real[], real, real, real, real, real&, real&, real&);

463 SphericalEngine::Value<true, SphericalEngine::SCHMIDT, 3>

464 (const coeff[], const real[], real, real, real, real, real&, real&, real&);

466 SphericalEngine::Value<false, SphericalEngine::SCHMIDT, 3>

467 (const coeff[], const real[], real, real, real, real, real&, real&, real&);

468

470 SphericalEngine::Circle<true, SphericalEngine::FULL, 1>

471 (const coeff[], const real[], real, real, real);

473 SphericalEngine::Circle<false, SphericalEngine::FULL, 1>

474 (const coeff[], const real[], real, real, real);

476 SphericalEngine::Circle<true, SphericalEngine::SCHMIDT, 1>

477 (const coeff[], const real[], real, real, real);

479 SphericalEngine::Circle<false, SphericalEngine::SCHMIDT, 1>

480 (const coeff[], const real[], real, real, real);

481

483 SphericalEngine::Circle<true, SphericalEngine::FULL, 2>

484 (const coeff[], const real[], real, real, real);

486 SphericalEngine::Circle<false, SphericalEngine::FULL, 2>

487 (const coeff[], const real[], real, real, real);

489 SphericalEngine::Circle<true, SphericalEngine::SCHMIDT, 2>

490 (const coeff[], const real[], real, real, real);

492 SphericalEngine::Circle<false, SphericalEngine::SCHMIDT, 2>

493 (const coeff[], const real[], real, real, real);

494

496 SphericalEngine::Circle<true, SphericalEngine::FULL, 3>

497 (const coeff[], const real[], real, real, real);

499 SphericalEngine::Circle<false, SphericalEngine::FULL, 3>

500 (const coeff[], const real[], real, real, real);

502 SphericalEngine::Circle<true, SphericalEngine::SCHMIDT, 3>

503 (const coeff[], const real[], real, real, real);

505 SphericalEngine::Circle<false, SphericalEngine::SCHMIDT, 3>

506 (const coeff[], const real[], real, real, real);

507

508

509}

Header for GeographicLib::CircularEngine class.

#define GEOGRAPHICLIB_EXPORT

Header for GeographicLib::SphericalEngine class.

Header for GeographicLib::Utility class.

Spherical harmonic sums for a circle.

Exception handling for GeographicLib.

Package up coefficients for SphericalEngine.

Math::real Sv(int k) const

static int Csize(int N, int M)

static int Ssize(int N, int M)

static void readcoeffs(std::istream &stream, int &N, int &M, std::vector< real > &C, std::vector< real > &S, bool truncate=false)

Definition SphericalEngine.cpp:384

Math::real Cv(int k) const

static Math::real Value(const coeff c[], const real f[], real x, real y, real z, real a, real &gradx, real &grady, real &gradz)

Definition SphericalEngine.cpp:152

static void RootTable(int N)

Definition SphericalEngine.cpp:373

static CircularEngine Circle(const coeff c[], const real f[], real p, real z, real a)

Definition SphericalEngine.cpp:296

static std::string str(T x, int p=-1)

Namespace for GeographicLib.