GeographicLib  2.6
GeodesicLine3.cpp
Go to the documentation of this file.
1 /**
2  * \file GeodesicLine3.cpp
3  * \brief Implementation for GeographicLib::Triaxial::GeodesicLine3 class
4  *
5  * Copyright (c) Charles Karney (2024-2025) <karney@alum.mit.edu> and licensed
6  * under the MIT/X11 License. For more information, see
7  * https://geographiclib.sourceforge.io/
8  **********************************************************************/
9 
10 #include <iostream>
11 #include <iomanip>
12 // Needed by zsetsdiag
13 #include <sstream>
14 #include <GeographicLib/Triaxial/GeodesicLine3.hpp>
15 
16 namespace GeographicLib {
17  namespace Triaxial {
18 
19  using namespace std;
20 
21  GeodesicLine3::GeodesicLine3(fline f, fline::fics fic,
22  gline g, gline::gics gic)
23  : _tg(f.tg())
24  , _f(f)
25  , _fic(fic)
26  , _g(g)
27  , _gic(gic)
28  {}
29 
30  GeodesicLine3::GeodesicLine3(const Geodesic3& tg)
31  : _f(tg, Geodesic3::gamblk(tg, (tg._umbalt &&
32  tg.kp2() > 0) || tg.k2() == 0))
33  {}
34 
35  GeodesicLine3::GeodesicLine3(const Geodesic3& tg,
36  Angle bet1, Angle omg1, Angle alp1)
37  : _tg(tg)
38  {
39  bet1.round();
40  omg1.round();
41  alp1.round();
42  Geodesic3::gamblk gam = _tg.gamma(bet1, omg1, alp1);
43  _f = fline(tg, gam);
44  _fic = fline::fics(_f, bet1, omg1, alp1);
45  _g = gline(tg, gam);
46  _gic = gline::gics(_g, _fic);
47  }
48 
49  GeodesicLine3::GeodesicLine3(const Geodesic3& tg,
50  real bet1, real omg1, real alp1)
51  : GeodesicLine3(tg, ang(bet1), ang(omg1), ang(alp1))
52  {}
53 
54  void GeodesicLine3::pos1(Angle& bet1, Angle& omg1, Angle& alp1) const {
55  _fic.pos1(_f.transpolar(), bet1, omg1, alp1);
56  }
57 
58  void GeodesicLine3::pos1(real& bet1, real& omg1, real& alp1,
59  bool unroll) const {
60  ang bet1a, omg1a, alp1a;
61  pos1(bet1a, omg1a, alp1a);
62  if (!unroll) {
63  (void) Ellipsoid3::AngNorm(bet1a, omg1a, alp1a);
64  bet1a.setn(); omg1a.setn(); alp1a.setn();
65  }
66  bet1 = real(bet1a);
67  omg1 = real(omg1a);
68  alp1 = real(alp1a);
69  }
70 
71  void GeodesicLine3::Position(real s12,
72  Angle& bet2a, Angle& omg2a, Angle& alp2a) const {
73  // Compute points at distance s12
74  real sig2 = _gic.sig1 + s12/_tg.b();
75  ang &phi2a = bet2a, &tht2a = omg2a;
76  int *countn = nullptr, *countb = nullptr;
77  if (_f.gammax() > 0) {
78  real u2, v2;
79  if constexpr (false && biaxspecial(_tg, _g.gammax())) {
80  u2 = gpsi().inv(sig2, countn, countb);
81  v2 = ftht().inv(fpsi()(u2) - _fic.delta, countn, countb);
82  } else
83  // The general triaxial machinery. This is used for non-meridional
84  // geodesics on biaxial ellipsoids.
85  solve2(_fic.delta, sig2, fpsi(), ftht(), gpsi(), gtht(), u2, v2,
86  countn, countb);
87  tht2a = ang::radians(ftht().rev(v2));
88  ang psi2 = ang::radians(fpsi().rev(u2));
89  // Already normalized
90  phi2a = ang(_f.gm().nup * psi2.s(),
91  _fic.phi0.c() * hypot(psi2.c(), _f.gm().nu * psi2.s()),
92  0, true).rebase(_fic.phi0);
93  alp2a = ang(_fic.Ex * hypot(_f.kx() * _f.gm().nu, _f.kxp() * tht2a.c()),
94  _fic.phi0.c() * _f.kx() * _f.gm().nup * psi2.c());
95  } else if (_f.gammax() == 0) {
96  pair<real, real> sig2n = remx(sig2, 2*_g.s0); // reduce to [-s0, s0)
97  real u2, v2;
98  ang psi2;
99  if (_f.kxp2() == 0) {
100  // meridr
101  // ftht = tht
102  // gtht = 0
103  // meridr
104  // fpsi = 0
105  // gpsi = int(sqrt(1-eps*cos(psi)^2), psi)
106  solve2(_fic.delta, sig2n.first,
107  fpsi(), ftht(), gpsi(), gtht(), u2, v2,
108  countn, countb);
109  phi2a = ang::radians(u2);
110  // u2 in in [-pi/2, pi/2]. With long doubles cos(pi/2) < 0 which leads
111  // to a switch in sheets in ellipsoidal coordinates. Fix by setting
112  // cos(phi2a) = +eps
113  if (signbit(phi2a.c())) // Not triggered with doubles and quads
114  phi2a = ang(copysign(real(1), phi2a.s()),
115  numeric_limits<real>::epsilon()/(1<<11), 0, true);
116  psi2 = phi2a + ang::cardinal(2 * sig2n.second);
117  int parity = fmod(sig2n.second, real(2)) != 0 ? -1 : 1;
118  int Ny = _fic.Nx * parity;
119  phi2a = phi2a.reflect(signbit(_fic.phi0.c() * Ny),
120  signbit(_fic.phi0.c())).rebase(_fic.phi0);
121  tht2a = ang::radians(- _fic.delta) + ang::cardinal(2 * sig2n.second);
122  alp2a = ang(_fic.Ex * real(0), _fic.Nx * parity, 0, true);
123  } else {
124  if (sig2n.first - _g.s0 >= -5 * numeric_limits<real>::epsilon()) {
125  sig2n.first = -_g.s0;
126  ++sig2n.second;
127  }
128  real deltax = bigclamp(_fic.delta + sig2n.second * _f.deltashift(), 1);
129  solve2u(deltax, sig2n.first, fpsi(), ftht(), gpsi(), gtht(), u2, v2,
130  countn, countb);
131  // phi2 = fpsi().rev(u2); tht2 = ftht().rev(v2);
132  phi2a = anglam(u2, _f.kxp());
133  psi2 = phi2a + ang::cardinal(2 * sig2n.second);
134  tht2a = anglam(v2, _f.kx());
135  int parity = fmod(sig2n.second, real(2)) != 0 ? -1 : 1;
136  // if_tg.t().kp2 == 0 then meridional oblate
137  int Ny = _fic.Nx * parity;
138  tht2a += ang::cardinal(2 * sig2n.second);
139  tht2a = tht2a + _fic.tht0;
140  phi2a = phi2a.reflect(signbit(_fic.phi0.c() * Ny),
141  signbit(_fic.phi0.c())).rebase(_fic.phi0);
142  // replace cos(phi)/cos(tht) by sech(u)/sech(v)
143  alp2a = ang(_fic.Ex * _f.kxp() / mcosh(v2, _f.kx()),
144  _f.kx() * Ny / mcosh(u2, _f.kxp()));
145  }
146  } else {
147  // gamma = NaN
148  }
149  phi2a.round();
150  tht2a.round();
151  alp2a.round();
152  tht2a = tht2a.flipsign(_fic.Ex);
153  if (_f.transpolar()) {
154  swap(bet2a, omg2a);
155  alp2a.reflect(false, false, true);
156  }
157  alp2a = alp2a.rebase(_fic.alp0);
158  omg2a += ang::cardinal(1);
159  }
160 
161  void GeodesicLine3::Position(real s12,
162  real& bet2, real& omg2, real& alp2,
163  bool unroll) const {
164  ang bet2a, omg2a, alp2a;
165  Position(s12, bet2a, omg2a, alp2a);
166  if (!unroll) {
167  (void) Ellipsoid3::AngNorm(bet2a, omg2a, alp2a);
168  bet2a.setn(); omg2a.setn(); alp2a.setn();
169  }
170  bet2 = real(bet2a);
171  omg2 = real(omg2a);
172  alp2 = real(alp2a);
173  }
174 
175  void GeodesicLine3::solve2(real f0, real g0,
176  const hfun& fx, const hfun& fy,
177  const hfun& gx, const hfun& gy,
178  real& x, real& y,
179  int* countn, int* countb) {
180  // Return x and y, s.t.
181  // fx(x) - fy(y) = f0
182  // gx(x) + gy(y) = g0
183  //
184  // fx(x) = fxs*x +/- fxm,
185  // fy(y) = fys*y +/- fym,
186  // gx(x) = gxs*x +/- gxm,
187  // gy(y) = gys*y +/- gym;
188  const bool check = false;
189  real fxm = fx.MaxPlus(), fym = fy.MaxPlus(),
190  gxm = gx.MaxPlus(), gym = gy.MaxPlus(),
191  fxs = fx.Slope(), fys = fy.Slope(), gxs = gx.Slope(), gys = gy.Slope(),
192  // solve
193  // x = ( fys*g0 + gys*f0 ) / den +/- Dx
194  // y = ( fxs*g0 - gxs*f0 ) / den +/- Dy
195  // where
196  den = fxs * gys + fys * gxs, // den > 0
197  qf = fxm + fym, qg = gxm + gym,
198  Dx = (qf * gys + qg * fys) / den,
199  Dy = (qf * gxs + qg * fxs) / den,
200  x0 = (fys * g0 + gys * f0) / den, // Initial guess
201  y0 = (fxs * g0 - gxs * f0) / den,
202  xp = x0 + Dx, xm = x0 - Dx,
203  yp = y0 + Dy, ym = y0 - Dy,
204  mm = 0;
205  if constexpr (check) {
206  real DxA = (-qf * gys + qg * fys) / den,
207  DyA = (-qf * gxs + qg * fxs) / den;
208  real
209  fA = (fx(x0 - Dx) - fy(y0 - DyA)) - f0,
210  gA = (gx(x0 - Dx) + gy(y0 - DyA)) - g0,
211  fB = (fx(x0 - DxA) - fy(y0 - Dy)) - f0,
212  gB = (gx(x0 - DxA) + gy(y0 - Dy)) - g0,
213  fC = (fx(x0 + Dx) - fy(y0 + DyA)) - f0,
214  gC = (gx(x0 + Dx) + gy(y0 + DyA)) - g0,
215  fD = (fx(x0 + DxA) - fy(y0 + Dy)) - f0,
216  gD = (gx(x0 + DxA) + gy(y0 + Dy)) - g0;
217  if (!( fabs(DxA) <= Dx && fabs(DyA) <= Dy ))
218  throw GeographicLib::GeographicErr("Bad Dx/Dy");
219  if (!( fA <= 0 && gA <= 0 )) {
220  cout << scientific << "midA " << fA << " " << gA << "\n";
221  cout << "DA " << Dx << " " << DxA << " " << Dy << " " << DyA << "\n";
222  throw GeographicLib::GeographicErr
223  ("Bad initial midpoints A GeodesicLine3::newt2");
224  }
225  if (!( fB >= 0 && gB <= 0 )) {
226  cout << scientific << "midB " << fB << " " << gB << "\n";
227  throw GeographicLib::GeographicErr
228  ("Bad initial midpoints B GeodesicLine3::newt2");
229  }
230  if (!( fC >= 0 && gC >= 0 )) {
231  cout << scientific << "midC " << fC << " " << gC << "\n";
232  throw GeographicLib::GeographicErr
233  ("Bad initial midpoints C GeodesicLine3::newt2");
234  }
235  if (!( fD <= 0 && gD >= 0 )) {
236  cout << scientific << "midD " << fD << " " << gD << "\n";
237  throw GeographicLib::GeographicErr
238  ("Bad initial midpoints D GeodesicLine3::newt2");
239  }
240  }
241  newt2(f0, g0, fx, fy, gx, gy,
242  xm-mm*Dx, xp+mm*Dx, fx.HalfPeriod(),
243  ym-mm*Dy, yp+mm*Dy, fy.HalfPeriod(),
244  (fx.HalfPeriod() * fxs + fy.HalfPeriod() * fys) / 2,
245  (gx.HalfPeriod() * gxs + gy.HalfPeriod() * gys) / 2,
246  x, y, countn, countb);
247  }
248 
249  void GeodesicLine3::solve2u(real d0, real s0,
250  const hfun& fx, const hfun& fy,
251  const hfun& gx, const hfun& gy,
252  real& u, real& v,
253  int* countn, int* countb) {
254  // Return u and v, s.t.
255  // fx(u) - fy(v) = d0
256  // gx(u) + gy(v) = s0
257  // specialized for umbilics
258  //
259  // fx, fy, gx, gy are increasing functions defined in [-1, 1]*pi2
260  // Assume fx(0) = fy(0) = gx(0) = gy(0) = 0
261  real pi2 = Geodesic3::BigValue(),
262  sbet = gx.Max(), somg = gy.Max(), stot = sbet + somg,
263  dbet = fx.Max(), domg = fy.Max(), del = dbet - domg;
264  if (fabs(s0) - stot >= -5 * numeric_limits<real>::epsilon()) {
265  // close to umbilic points we have
266  // fx(u) = u -/+ dbet
267  // fy(v) = v -/+ domg
268  // fx(u) - fy(v) = d0
269  // constrain u+v = +/- 2 * pi2
270  // u = +/- pi2 + (d0 +/- (dbet - domg))/2
271  // v = +/- pi2 - (d0 +/- (dbet - domg))/2
272  real t0 = copysign(pi2, s0),
273  t1 = (d0 + (1 - 2 * signbit(s0)) * del) / 2;
274  u = t0 + t1; v = t0 - t1;
275  } else if (fabs(d0) > 2*pi2/3 &&
276  fabs((1 - 2 * signbit(d0)) * s0 - (sbet - somg)) <=
277  7 * numeric_limits<real>::epsilon()) {
278  if (d0 > 0) {
279  u = 2*d0/3; v = -1*d0/3;
280  } else {
281  u = 1*d0/3; v = -2*d0/3;
282  }
283  } else {
284  real mm = 2;
285  newt2(d0, s0, fx, fy, gx, gy,
286  -mm * pi2, mm * pi2, pi2,
287  -mm * pi2, mm * pi2, pi2,
288  pi2, pi2, u, v, countn, countb);
289  }
290  }
291 
292  void GeodesicLine3::newt2(real f0, real g0,
293  const hfun& fx, const hfun& fy,
294  const hfun& gx, const hfun& gy,
295  real xa, real xb, real xscale,
296  real ya, real yb, real yscale,
297  real fscale, real gscale,
298  real& x, real& y,
299  int* countn, int* countb) {
300  // solve
301  // f = fx(x) - fy(y) - f0 = 0
302  // g = gx(x) + gy(y) - g0 = 0
303  // for x, y. Assume that
304  // fx and fy are increasing functions
305  // gx and gy are non-decreasing functions
306  // The solution is bracketed by x in [xa, xb], y in [ya, yb]
307  static const real tol = numeric_limits<real>::epsilon();
308  const bool debug = Geodesic3::debug_, check = false;
309  // Relax [fg]tol to /10 instead of /100. Otherwise solution resorts to
310  // too much bisection. Example:
311  // echo 63 -173 -61.97921997838416712 -4.64409746197940890408 |
312  // ./Geod3Solve $PROX
313  const real
314  ftol = tol * fscale/10,
315  gtol = tol * gscale/10,
316  xtol = pow(tol, real(0.75)) * xscale,
317  ytol = pow(tol, real(0.75)) * yscale,
318  // If tolmult == 0, umbilical case
319  // echo 70 0 180 0.46138515584816834054 | ./Geod3Solve $SET
320  // take countn/countb = 55/54 iterations to converge. With tolmult == 1,
321  // it needs 5/3 iterations.
322  tolmult = 1;
323  int cntn = 0, cntb = 0;
324  real oldf = Math::infinity(), oldg = oldf, olddx = oldf, olddy = oldf;
325  zset xset(zvals(xa, fx(xa), gx(xa)),
326  zvals(xb, fx(xb), gx(xb))),
327  yset(zvals(ya, fy(ya), gy(ya)),
328  zvals(yb, fy(yb), gy(yb)));
329  // x = xset.bisect(), y = yset.bisect();
330  auto p = zsetsbisect(xset, yset, f0, g0, false);
331  x = p.first; y = p.second;
332  if constexpr (check) {
333  // A necessary condition for a root is
334  // f01 <= 0 <= f10
335  // g00 <= 0 <= g11
336  real
337  f01 = (xset.min().fz - yset.max().fz) - f0,
338  f10 = (xset.max().fz - yset.min().fz) - f0,
339  g00 = (xset.min().gz + yset.min().gz) - g0,
340  g11 = (xset.max().gz + yset.max().gz) - g0;
341  // Allow equality on the initial points
342  if (!( f01 <= 0 && 0 <= f10 && g00 <= 0 && 0 <= g11 ))
343  throw GeographicLib::GeographicErr
344  ("Bad initial points GeodesicLine3::newt2");
345  zvals xv(x, fx(x), gx(x)), yv(y, fy(y), gy(y));
346  }
347  // degen is a flag detected degeneracy to a 1d system. Only the case gy(y)
348  // = const is treated. Then g = gx(x) + gy(y) - g0 = 0 converges as a 1d
349  // Newton system which fixes the value of x. Once x is found (detected by
350  // xset.min().z and xset.max().z being consecutive numbers), g0 is adjusted
351  // to that the g equation is satisfied exactly. This allows f = fx(x) -
352  // fy(y) - f0 = 0 to be solved reliably for y.
353  bool bis = false, degen = false;
354  int ibis = -1, i = 0;
355  for (; i < maxit_ ||
356  (Geodesic3::throw_ && (throw GeographicLib::GeographicErr
357  ("Convergence failure GeodesicLine3::newt2"),
358  false));
359  ++i) {
360  ++cntn;
361  if (!degen && nextafter(xset.min().z, xset.max().z) == xset.max().z &&
362  yset.min().gz == yset.max().gz) {
363  degen = true;
364  real ga = (xset.min().gz + yset.min().gz) - g0,
365  gb = (xset.max().gz + yset.max().gz) - g0;
366  x = gb < -ga ? xset.max().z : xset.min().z;
367  // Adjust g0 so that g = 0 and dx = 0
368  if (gb < -ga) {
369  g0 = xset.max().gz + yset.max().gz;
370  zvals xx = xset.max();
371  xset.insert(xx, -1);
372  } else {
373  g0 = xset.min().gz + yset.min().gz;
374  zvals xx = xset.min();
375  xset.insert(xx, +1);
376  }
377  }
378  zvals xv(x, fx(x), gx(x)), yv(y, fy(y), gy(y));
379  // zsetsinsert updates xv and yv to enforce monotonicity of f and g
380  zsetsinsert(xset, yset, xv, yv, f0, g0);
381  real f = (xv.fz - yv.fz) - f0, g = (xv.gz + yv.gz) - g0;
382  if ((fabs(f) <= ftol && fabs(g) <= gtol) || isnan(f) || isnan(g)) {
383  if constexpr (debug)
384  cout << "break0 " << scientific << f << " " << g << "\n";
385  break;
386  }
387  real
388  fxp = fx.deriv(x), fyp = fy.deriv(y),
389  gxp = gx.deriv(x), gyp = gy.deriv(y),
390  den = fxp * gyp + fyp * gxp,
391  dx = -( gyp * f + fyp * g) / den, // if degen, dx = 0
392  dy = -(-gxp * f + fxp * g) / den, // if degen dy = f / fyp
393  xn = x + dx, yn = y + dy,
394  dxa = 0, dya = 0;
395  xa = xset.min().z; xb = xset.max().z;
396  ya = yset.min().z; yb = yset.max().z;
397  if constexpr (debug) {
398  bool bb = gyp == 0 &&
399  nextafter(xset.min().z, xset.max().z) == xset.max().z;
400  cout << "DERIV " << i << " " << fxp << " " << fyp << " " << gxp << " " << gyp << " " << bb << "\n";
401  cout << "DY " << i << " " << -gxp * f << " " << fxp * g << "\n";
402  cout << "FG " << i << " " << f << " " << g << "\n";
403  cout << "DXY " << i << " " << degen << " " << dx << " " << dy << "\n";
404  }
405  if constexpr (check) {
406  if (!( fxp >= 0 && fyp >= 0 && gxp >= 0 && gyp >= 0 && den > 0 )) {
407  cout << "DERIVS " << x << " " << y << " "
408  << fxp << " " << fyp << " "
409  << gxp << " " << gyp << " " << den << "\n";
410  throw GeographicLib::GeographicErr
411  ("Bad derivatives GeodesicLine3::newt2");
412  }
413  }
414  bool cond1 = den > 0 &&
415  (i < ibis + 2 ||
416  ((2*fabs(f) < oldf || 2*fabs(g) < oldg) ||
417  (2*fabs(dx) < olddx || 2*fabs(dy) < olddy))),
418  cond2 = xn >= xa-xtol*tolmult && xn <= xb+xtol*tolmult &&
419  yn >= ya-ytol*tolmult && yn <= yb+ytol*tolmult;
420  if (cond1 && cond2) {
421  oldf = fabs(f); oldg = fabs(g); olddx = fabs(dx); olddy = fabs(dy);
422  x = xn; y = yn;
423  bis = false;
424  if (!(fabs(dx) > xtol || fabs(dy) > ytol)) {
425  if constexpr (debug)
426  cout << "break1 " << scientific << dx << " " << dy << " "
427  << f << " " << g << "\n";
428  break;
429  }
430  } else {
431  // xn = xset.bisect(); yn = yset.bisect();
432  p = zsetsbisect(xset, yset, f0, g0, false);
433  xn = p.first; yn = p.second;
434  ++cntb;
435  if (x == xn && y == yn) {
436  if constexpr (debug)
437  cout << "break2 " << f << " " << g << "\n";
438  break;
439  }
440  dxa = xn - x; dya = yn - y;
441  x = xn; y = yn;
442  bis = true;
443  ibis = i;
444  }
445  (void) bis;
446  if constexpr (debug)
447  cout << "AA " << scientific << setprecision(4)
448  << x-xa << " " << xb-x << " "
449  << y-ya << " " << yb-y << "\n";
450  if constexpr (debug)
451  cout << "CC " << i << " "
452  << bis << " " << cond1 << " " << cond2 << " "
453  << scientific << setprecision(2) << f << " " << g << " "
454  << dx << " " << dy << " " << dxa << " " << dya << " "
455  << xb-xa << " " << yb-ya << "\n";
456  if constexpr (debug)
457  cout << "BOX " << i << " " << bis << " "
458  << xset.num() << " " << yset.num() << " "
459  << scientific << setprecision(3)
460  << xset.max().z - xset.min().z << " "
461  << yset.max().z - yset.min().z << "\n";
462  }
463  if (countn)
464  *countn += cntn;
465  if (countb)
466  *countb += cntb;
467  if constexpr (debug) {
468  cout << "CNT " << cntn << " " << cntb << "\n";
469  cout << "XY " << setprecision(18) << x << " " << y << "\n";
470  }
471  }
472 
473  int GeodesicLine3::zset::insert(zvals& t, int flag) {
474  // Inset t into list. flag = -/+ 1 indicates new min/max.
475  // Return -1 if t was already present; othersize return index of newly
476  // inserted value.
477  // If monotonic, value of t.fz and t.gz is adjusted to ensuire monotonicity:
478  // if a.z == b.z then a.fz == b.gz && a.fz == b.gz
479  // if a.z < b.z then a.fz <= b.gz && a.fz <= b.gz
480  using std::isnan;
481  // Need this flag to be set otherwise conditions for setting the bounds on
482  // the allowed results are not met and newt2 fails to converge.
483  const bool monotonic = true;
484  int ind = -1;
485  if (isnan(t.z)) return ind;
486  if (t < min()) {
487  if constexpr (monotonic) {
488  t.fz = fmin(t.fz, min().fz);
489  t.gz = fmin(t.gz, min().gz);
490  }
491  } else if (t == min()) {
492  // Check if t is "other" endpoint and collapse bracket to zero
493  if (flag > 0) _s.resize(1);
494  t = min();
495  } else if (t == max()) {
496  // Check if t is "other" endpoint and collapse bracket to zero
497  if (flag < 0) { _s[0] = _s.back(); _s.resize(1); }
498  t = max();
499  } else if (max() < t) {
500  if constexpr (monotonic) {
501  t.fz = fmax(t.fz, max().fz);
502  t.gz = fmax(t.gz, max().gz);
503  }
504  }
505  if (!(min() < t && t < max())) // Not in range
506  return ind;
507  // Now min() < t < max()
508  auto p = std::lower_bound(_s.begin(), _s.end(), t);
509  if (p == _s.end()) return ind; // Can't happen
510  // Fix components of t
511  bool ins = !(*p == t);
512  if constexpr (monotonic) {
513  if (ins) {
514  t.fz = Math::clamp(t.fz, (p-1)->fz, p->fz);
515  t.gz = Math::clamp(t.gz, (p-1)->gz, p->gz);
516  } else // z components match
517  t = *p; // set fz and gz values
518  }
519  if (flag < 0) {
520  _s.erase(_s.begin(), p);
521  if (ins) {
522  _s.insert(_s.begin(), t);
523  ind = 0;
524  }
525  } else if (flag > 0) {
526  if (ins) {
527  _s.erase(p, _s.end());
528  _s.push_back(t);
529  ind = int(_s.size()) - 1;
530  } else
531  _s.erase(p+1, _s.end());
532  } else if (ins) {
533  ind = int(p - _s.begin());
534  _s.insert(p, t);
535  }
536  // else it's a duplicate and not a new end value
537  return ind;
538  }
539 
540  void GeodesicLine3::zsetsinsert(zset& xset, zset& yset,
541  zvals& xfg, zvals& yfg,
542  real f0, real g0) {
543  const bool debug = Geodesic3::debug_;
544  real x0 = 0, y0 = 0;
545  int xind = xset.insert(xfg), yind = yset.insert(yfg);
546  if constexpr (debug) {
547  cout << "BOXA " << xset.num() << " " << yset.num() << " "
548  << xind << " " << yind << " "
549  << xset.min().z - x0 << " " << xset.max().z - x0 << " "
550  << yset.min().z - y0 << " " << yset.max().z - y0 << " "
551  << xset.max().z - xset.min().z << " "
552  << yset.max().z - yset.min().z << "\n";
553  zsetsdiag(xset, yset, f0, g0);
554  }
555  if (xind < 0 && yind < 0) return;
556  // We update two sets of bounds xa[0], xb[0], etc. use the values at the
557  // newly inserted row and column. xa[1], xb[1], etc. use just the value at
558  // the center point.
559  zvals xa[2] = {xset.min(), xset.min()},
560  xb[2] = {xset.max(), xset.max()},
561  ya[2] = {yset.min(), yset.min()},
562  yb[2] = {yset.max(), yset.max()};
563  for (int i = 0; i < xset.num(); ++i) {
564  const zvals& x = xset.val(i);
565  for (int j = 0; j < yset.num(); ++j) {
566  if (i == xind || j == yind) {
567  const zvals& y = yset.val(j);
568  real f = (x.fz - y.fz) - f0, g = (x.gz + y.gz) - g0;
569  // Update bounds as follows
570  //
571  // y=-x y=x
572  // g=0 f=0
573  // \ ind=-2 /
574  // \ f<0 /
575  // \ g>0 /
576  // \ yb /
577  // \ D /
578  // ind=-4 \ / ind=4
579  // f<0 \ / f>0
580  // g<0 X g>0
581  // xa / \ xb
582  // A / \ C
583  // / B \.
584  // / f>0 \.
585  // / g<0 \.
586  // / ya \.
587  // / ind=2 \.
588 
589  // If the pattern is
590  //
591  // BOXF --+
592  // BOXF -++
593  // BOXF -++
594  //
595  // BOXG ..+
596  // BOXG ..+
597  // BOXG ---
598  //
599  // The, allowing equality (tight bounds) collapses the set to 1x1
600  // (the middle element). This is clearly wrong since now we can't
601  // get f = 0. We detect this by checking the values of f and g at
602  // the corners:
603  //
604  // f <= 0 at NW corner f >= 0 at SE corner
605  // g <= 0 at SW corner g >= 0 at NE corner
606  //
607  // If theses requirements fail, use bounds updated with just the
608  // center point i == xind && j == yind.
609  if (f <= 0) {
610  if (g <= 0 && xa[0] < x) xa[0] = x;
611  if (g >= 0 && y < yb[0]) yb[0] = y;
612  if (i == xind && j == yind) {
613  if (g <= 0 && xa[1] < x) xa[1] = x;
614  if (g >= 0 && y < yb[1]) yb[1] = y;
615  }
616  }
617  if (f >= 0) {
618  if (g <= 0 && ya[0] < y) ya[0] = y;
619  if (g >= 0 && x < xb[0]) xb[0] = x;
620  if (i == xind && j == yind) {
621  if (g <= 0 && ya[1] < y) ya[1] = y;
622  if (g >= 0 && x < xb[1]) xb[1] = x;
623  }
624  }
625  }
626  }
627  }
628  int k = 0;
629  for (; k < 2; ++k) {
630  // Check signs at corners
631  if ((xa[k].fz - yb[k].fz) - f0 <= 0 &&
632  (xb[k].fz - ya[k].fz) - f0 >= 0 &&
633  (xa[k].gz + ya[k].gz) - g0 <= 0 &&
634  (xb[k].gz + yb[k].gz) - g0 >= 0)
635  // corner constraints met
636  break;
637  }
638  int kk = min(k, 1);
639  if constexpr (debug)
640  cout << "BOXK " << k << " " << xind << " " << yind << "\n";
641  xset.insert(xa[kk], -1);
642  xset.insert(xb[kk], +1);
643  yset.insert(ya[kk], -1);
644  yset.insert(yb[kk], +1);
645  if constexpr (debug) {
646  cout << "BOXB " << xset.num() << " " << yset.num() << " "
647  << xset.min().z - x0 << " " << xset.max().z - x0 << " "
648  << yset.min().z - y0 << " " << yset.max().z - y0 << " "
649  << xa[kk].z - x0 << " " << xb[kk].z - x0 << " "
650  << ya[kk].z - y0 << " " << yb[kk].z - y0 << "\n";
651  zsetsdiag(xset, yset, f0, g0);
652  }
653  if (k == 2)
654  throw GeographicLib::GeographicErr
655  ("Bad corner points GeodesicLine3::zsetsinsert");
656  }
657 
658  void GeodesicLine3::zsetsdiag(const zset& xset, const zset& yset,
659  real f0, real g0) {
660  ostringstream fs, gs;
661  for (int j = yset.num() - 1; j >= 0; --j) {
662  const zvals& y = yset.val(j);
663  fs << "BOXF ";
664  gs << "BOXG ";
665  for (int i = 0; i < xset.num(); ++i) {
666  const zvals& x = xset.val(i);
667  real f = (x.fz - y.fz) - f0, g = (x.gz + y.gz) - g0;
668  fs << (f == 0 ? '.' : f < 0 ? '-' : '+');
669  gs << (g == 0 ? '.' : g < 0 ? '-' : '+');
670  }
671  fs << "\n";
672  gs << "\n";
673  }
674  cout << fs.str() << gs.str();
675  }
676 
677  pair<Math::real, Math::real>
678  GeodesicLine3::zsetsbisect(const zset& xset, const zset& yset,
679  real f0, real g0, bool secant) {
680  if constexpr (true)
681  return pair<real, real>(xset.bisect(), yset.bisect());
682  else if (secant && xset.num() <= 2 && yset.num() <= 2) {
683  // Use secant solution
684  real
685  dx = xset.max().z - xset.min().z,
686  dy = yset.max().z - yset.min().z,
687  fx1 = (xset.max().fz - xset.min().fz) / (dx == 0 ? 1 : dx),
688  fy1 = (yset.max().fz - yset.min().fz) / (dy == 0 ? 1 : dy),
689  gx1 = (xset.max().gz - xset.min().gz) / (dx == 0 ? 1 : dx),
690  gy1 = (yset.max().gz - yset.min().gz) / (dy == 0 ? 1 : dy),
691  den = fx1*gy1 + fy1*gx1,
692  f00 = (xset.min().fz + yset.min().fz) - f0,
693  g00 = (xset.min().gz + yset.min().gz) - g0,
694  x = -gy1*f00 - fy1*g00,
695  y = gx1*f00 - fx1*g00;
696  if (den <= 0) den = 1;
697  x /= den;
698  y /= den;
699  if (x <= 0 || x >= 1) x = 1/real(2);
700  if (y <= 0 || y >= 1) y = 1/real(2);
701  x = xset.min().z + x * dx;
702  y = yset.min().z + y * dy;
703  return pair<real, real>(x, y);
704  } else {
705  // A necessary condition for a root is
706  // f01 <= 0 <= f10
707  // g00 <= 0 <= g11
708  int cnt = 0;
709  real xgap = -1, ygap = -1, x = Math::NaN(), y = Math::NaN();
710  for (int i = 0; i < max(1, xset.num() - 1); ++i) {
711  int i1 = min(i + 1, xset.num() - 1);
712  real xgap1 = xset.val(i1).z - xset.val(i).z,
713  xmean = (xset.val(i1).z + xset.val(i).z) / 2;
714  for (int j = 0; j < max(1, yset.num() - 1); ++j) {
715  int j1 = min(j + 1, yset.num() - 1);
716  real ygap1 = yset.val(j1).z - yset.val(j).z,
717  ymean = (yset.val(j1).z + yset.val(j).z) / 2;
718  real
719  f01 = (xset.val(i ).fz - yset.val(j1).fz) - f0,
720  f10 = (xset.val(i1).fz - yset.val(j ).fz) - f0,
721  g00 = (xset.val(i ).gz + yset.val(j ).gz) - g0,
722  g11 = (xset.val(i1).gz + yset.val(j1).gz) - g0;
723  if (f01 <= 0 && 0 <= f10 && g00 <= 0 && 0 <= g11) {
724  ++cnt;
725  if (xgap1 > xgap) { xgap = xgap1; x = xmean; }
726  if (ygap1 > ygap) { ygap = ygap1; y = ymean; }
727  }
728  }
729  }
730  if (cnt == 0)
731  throw GeographicLib::GeographicErr
732  ("No legal box GeodesicLine3::zsetsbisect");
733  return pair<real, real>(x, y);
734  }
735  }
736 
737  void GeodesicLine3::Hybrid(ang betomg2,
738  ang& bet2a, ang& omg2a, ang& alp2a,
739  real& s12, bool betp) const {
740  fline::disttx d = _f.Hybrid(_fic, betomg2, bet2a, omg2a, alp2a, betp);
741  s12 = _g.dist(_gic, d);
742  }
743 
744  GeodesicLine3::fline::disttx
745  GeodesicLine3::fline::Hybrid(const fics& fic, ang betomg2,
746  ang& bet2a, ang& omg2a, ang& alp2a,
747  bool betp) const {
748  // Is the control variable psi or tht?
749  bool psip = !transpolar() ? betp : !betp;
750  if (!betp) betomg2 -= ang::cardinal(1);
751  ang tau12;
752  if (psip) {
753  ang phi2 = betomg2;
754  if (gammax() > 0) {
755  real spsi = phi2.s(),
756  // In evaluating equivalent expressions, choose the one with minimum
757  // cancelation. Need the 0 + x to convert -0 to +0. (Note sqrt(-0) =
758  // -0 and fmax(+0, -0) may be -0.)
759  cpsi = nu() < nup() ?
760  (phi2.c() - nu()) * (phi2.c() + nu()) :
761  (nup() - phi2.s()) * (nup() + phi2.s());
762  // Return nan if geodesic can't reach phi2 -- given by sqrt(neg) -- but
763  // allow a little slop.
764  cpsi = !(cpsi > -numeric_limits<real>::epsilon()) ? Math::NaN() :
765  (signbit(cpsi) ? 0 : sqrt(cpsi));
766  // Need Angle(0, 0) to be treated like Angle(0, 1) here.
767  ang psi12 = (ang(spsi, cpsi) - fic.psi1).base();
768  // convert -180deg to 180deg
769  if (signbit(psi12.s()))
770  psi12 = ang(0, copysign(real(1), psi12.c()), 0, true);
771  tau12 = psi12;
772  } else if (gammax() == 0) {
773  tau12 = (fic.Nx > 0 ? phi2 - fic.phi1 :
774  phi2 + fic.phi1 + ang::cardinal(2)).base();
775  } else
776  tau12 = ang::NaN();
777  } else {
778  ang tht2 = betomg2.flipsign(fic.Ex);
779  if (gammax() >= 0) {
780  // FIX THIS!! betp shouldn't be appearing here.
781  // Test case
782  // echo -88 21 88 -111 | ./Geod3Solve -i $SET --hybridalt
783  ang tht2b = tht2; tht2b.reflect(false, betp && fic.Ex < 0);
784  ang tht12 = tht2b - fic.tht1;
785  // convert -180deg to 180deg
786  if (signbit(tht12.s()))
787  tht12 = ang(0, copysign(real(1), tht12.c()), 0, true);
788  tau12 = tht12;
789  } else
790  tau12 = ang::NaN();
791  }
792  disttx ret = ArcPos0(fic, tau12.base(), bet2a, omg2a, alp2a, betp);
793  return ret;
794  }
795 
796  GeodesicLine3::fline::fline(const Geodesic3& tg, bool neg)
797  : _tg(tg)
798  , _gm(tg, neg)
799  {}
800 
801  GeodesicLine3::fline::fline(const Geodesic3& tg, Geodesic3::gamblk gam)
802  : _tg(tg)
803  , _gm(gam)
804  , _fpsi(false, _gm.kx2 , _gm.kxp2,
805  +(_gm.transpolar ? -1 : 1) * _tg.e2(),
806  -(_gm.transpolar ? -1 : 1) * _gm.gamma, _tg)
807  , _ftht(false, _gm.kxp2 , _gm.kx2,
808  -(_gm.transpolar ? -1 : 1) * _tg.e2(),
809  +(_gm.transpolar ? -1 : 1) * _gm.gamma, _tg)
810  {
811  // Only needed for umbilical lines
812  _deltashift = _gm.gamma == 0 ?
813  (_tg.k2() > 0 && _tg.kp2() > 0 ? 2 * (_fpsi.Max() - _ftht.Max()) : 0) :
814  Math::NaN();
815  }
816 
817  GeodesicLine3::gline::gline(const Geodesic3& tg, bool neg)
818  : _tg(tg)
819  , _gm(tg, neg)
820  {}
821 
822  GeodesicLine3::gline::gline(const Geodesic3& tg, const Geodesic3::gamblk& gm)
823  : _tg(tg)
824  , _gm(gm)
825  , _gpsi(true, _gm.kx2 , _gm.kxp2,
826  +(_gm.transpolar ? -1 : 1) * _tg.e2(),
827  -(_gm.transpolar ? -1 : 1) * _gm.gamma, _tg)
828  , _gtht(true, _gm.kxp2 , _gm.kx2,
829  -(_gm.transpolar ? -1 : 1) * _tg.e2(),
830  +(_gm.transpolar ? -1 : 1) * _gm.gamma, _tg)
831  , s0(_gm.gammax == 0 ? _gpsi.Max() + _gtht.Max() : 0)
832  {}
833 
834  Math::real GeodesicLine3::fline::Hybrid0(const fics& fic, ang bet2, ang omg2,
835  bool betp) const {
836  ang bet2a, omg2a, alp2a;
837  (void) Hybrid(fic, betp ? bet2 : omg2, bet2a, omg2a, alp2a, betp);
838  bool angnorm = true || betp;
839  if (angnorm)
840  (void) Ellipsoid3::AngNorm(bet2a, omg2a, alp2a, !betp);
841  if (betp) {
842  omg2a -= omg2;
843  return omg2a.radians0();
844  } else {
845  bet2a -= bet2;
846  return angnorm ? bet2a.radians0() : bet2.radians();
847  }
848  }
849 
850  GeodesicLine3::fline::disttx
851  GeodesicLine3::fline::ArcPos0(const fics& fic, ang tau12,
852  ang& bet2a, ang& omg2a, ang& alp2a,
853  bool betp) const {
854  // XXX fix for biaxial
855  disttx ret{Math::NaN(), Math::NaN(), 0};
856  bool psip = transpolar() ? !betp : betp;
857  ang &phi2a = bet2a, &tht2a = omg2a;
858  if (gammax() > 0) {
859  ang psi2;
860  real u2, v2, u2x = 0;
861  if (psip) {
862  psi2 = tau12 + fic.psi1;
863  v2 = fpsi().fwd(psi2.radians());
864  u2 = ftht().inv(fpsi()(v2) - fic.delta);
865  tht2a = ang::radians(ftht().rev(u2));
866  } else {
867  tht2a = fic.tht1 + tau12;
868  u2 = ftht().fwd(tht2a.radians());
869  u2x = ftht()(u2) + fic.delta;
870  v2 = fpsi().inv(u2x);
871  psi2 = ang::radians(fpsi().rev(v2));
872  }
873  // Already normalized
874  phi2a = ang(nup() * psi2.s(),
875  fic.phi0.c() * hypot(psi2.c(), nu() * psi2.s()),
876  0, true).rebase(fic.phi0);
877  real s = fic.Ex * hypot(kx() * nu(), kxp() * tht2a.c()),
878  c = fic.phi0.c() * kx() * nup() * psi2.c();
879  if (s == 0 && c == 0)
880  (transpolar() ? s : c) = 1;
881  alp2a = ang(s, c);
882  ret.phiw2 = v2;
883  ret.thtw2 = u2;
884  } else if (gammax() == 0) {
885  real u2, v2;
886  int ii;
887  if (psip) {
888  phi2a = fic.phi1 + tau12.flipsign(fic.Nx);
889  // remainder with result in [-pi/2, pi/2)
890  pair<real, real> phi2n =
891  // remx(fic.Nx * (phi2a - fic.phi0).radians(), Math::pi());
892  remx((phi2a - fic.phi0).flipsign(fic.Nx));
893  u2 = fpsi().fwd(phi2n.first);
894  int parity = fmod(phi2n.second, real(2)) != 0 ? -1 : 1;
895  if (kxp() == 0) {
896  // v2 is independent on u2
897  v2 = 0;
898  tht2a = ang::radians(-fic.Ex * fic.delta) +
899  ang::cardinal(2*fic.Ex*phi2n.second);
900  alp2a = fic.alp1.nearest(2U) + ang::cardinal(parity < 0 ? 2 : 0);
901  } else {
902  real deltax = bigclamp(fic.delta + phi2n.second * _deltashift, 2);
903  v2 = ftht().inv(fpsi()(u2) - deltax);
904  tht2a = ang::radians(parity * ftht().rev(v2))
905  .rebase(fic.tht0);
906  // Conflict XXX
907  // testset -50 180 20 0 want fic.Nx multiplying s()
908  // testspha -20 90 20 -90 wants fic.Nx multiplying c()
909  alp2a = ang(kxp() * fic.Ex * parity / mcosh(v2, kx()),
910  fic.Nx * kx() / mcosh(u2, kxp()));
911  // Move forward from umbilical point
912  phi2a += ang::eps().flipsign(fic.Nx);
913  }
914  ii = int(phi2n.second);
915  } else {
916  tht2a = fic.tht1 + tau12;
917  // remainder with result in [-pi/2, pi/2)
918  pair<real, real> tht2n =
919  // remx(fic.Ex * (tht2a - fic.tht0).radians(), Math::pi());
920  remx(tht2a - fic.tht0);
921  v2 = ftht().fwd(tht2n.first);
922  u2 = 0;
923  int parity = fmod(tht2n.second, real(2)) != 0 ? -1 : 1;
924  if (kxp() == 0) {
925  if (fic.phi1.c() != 0 && tau12 == tau12.nearest(2U)) {
926  // STILL TO DO...
927  // Special case for finding conjugate points on meridional
928  // geodesics, gamma = 0, tau12 = multiple of pi. This is
929  // specialized for prolate ellipsoids for now.
930  real npi = (tau12.ncardinal() + fic.phi1.nearest(2U).ncardinal())
931  * Math::pi()/2;
932  // Don't worry about the case where we start at a pole -- this is
933  // already handled in Geodesic3::Inverse.
934  //
935  // Solve F(tpsi2) = tpsi2 - Deltaf(n*pi + atan(tpsi2))
936  // = (tan(psi1) - Deltaf(psi1)) = c
937  //
938  // for tpsi2, starting guess tpsi2 = tan(psi1).
939  // Test case prolate 1 3 0 1
940  // p1 = [-90, -1]
941  // [sx,a1x,a2x] = t.distance(p1,[90,178.9293750483])
942  // -> a1x = 90
943  // [sx,a1x,a2x] = t.distance(p1,[90,178.9293750484])
944  // -> a1x = 90.0023
945  // omg1 = -91 -> omg2 = 178.9293750483-90 = 88.9293750483
946  real c = fic.Nx * (fic.phi1.t() - fpsi().df(fic.phi1.radians())),
947  l = exp(Geodesic3::BigValue()),
948  tpsi2 = Trigfun::root(Trigfun::ARCPOS0,
949  [this, npi] (real tpsi) -> pair<real, real>
950  {
951  real psi = atan(tpsi);
952  return pair<real, real>
953  (tpsi - fpsi().df(npi + psi),
954  1 - fpsi().dfp(psi) /
955  (1 + Math::sq(tpsi)));
956  },
957  c, fic.psi1.t(), -l, l);
958  v2 = atan(tpsi2);
959  phi2a = (ang(tpsi2, 1) + fic.phi1.nearest(2U))
960  .flipsign(parity*fic.Nx).rebase(fic.phi0);
961  alp2a = ang::cardinal(fic.Nx * parity == 1 ? 0 : 2);
962  } else {
963  u2 = v2 == 0 ? 0 : copysign(Math::pi()/2, tht2n.first);
964  phi2a = ang::cardinal(fabs(v2) == 0
965  ? 0 : copysign(real(1), v2 * fic.Nx))
966  .rebase(fic.phi0);
967  alp2a = fabs(v2) == 0 ? ang::cardinal(2) :
968  fic.alp1.nearest(2U) +
969  ang::cardinal(parity == 1 ? 0 : 2) +
970  ang::radians(v2).flipsign(parity * phi2a.s());
971  }
972  } else {
973  real deltax = bigclamp(fic.delta + tht2n.second * _deltashift, 2);
974  u2 = fpsi().inv(ftht()(v2) + deltax);
975  real phi2 = fic.Nx * parity * fpsi().rev(u2);
976  phi2a = ang::radians(phi2);
977  alp2a = ang(fic.Ex * kxp() / mcosh(v2, kx()),
978  kx() * fic.Nx * parity / mcosh(u2, kxp()));
979  tht2a += ang::eps();
980  }
981  ii = int(tht2n.second);
982  // Move forward from umbilical point
983  }
984  ret.phiw2 = u2;
985  ret.thtw2 = v2;
986  ret.ind2 = ii;
987  } else {
988  // gamma == NaN
989  }
990 
991  tht2a = tht2a.flipsign(fic.Ex);
992  if (transpolar()) {
993  swap(bet2a, omg2a);
994  alp2a.reflect(false, false, true);
995  }
996  omg2a += ang::cardinal(1);
997  alp2a = alp2a.rebase(fic.alp0);
998  return ret;
999  }
1000 
1001  GeodesicLine3::fline::fics::fics(const fline& f, ang bet10, ang omg10,
1002  ang alp10)
1003  : tht1(omg10 - ang::cardinal(1))
1004  , phi1(bet10)
1005  , alp1(alp10)
1006  {
1007  static const real eps = numeric_limits<real>::epsilon();
1008  alp0 = alp1.nearest(f.transpolar() ? 2U : 1U);
1009  if (f.transpolar()) {
1010  swap(tht1, phi1);
1011  alp1.reflect(false, false, true);
1012  }
1013  const Geodesic3& tg = f.tg();
1014  if (!f.transpolar() && phi1.s() == 0 && fabs(alp1.c()) <= Math::sq(eps))
1015  alp1 = ang(alp1.s(), - Math::sq(eps), alp1.n(), true);
1016  Ex = signbit(alp1.s()) ? -1 : 1;
1017  Nx = signbit(alp1.c()) ? -1 : 1;
1018  tht1 = tht1.flipsign(Ex);
1019  if (f.gammax() > 0) {
1020  phi0 = phi1.nearest(2U);
1021  psi1 = ang(f.kx() * phi1.s(),
1022  phi0.c() * alp1.c() *
1023  hypot(f.kx() * phi1.c(), f.kxp() * tht1.c()));
1024  // k = 1, kp = 0, sqrt(-mu) = bet1.c() * fabs(alp1.s())
1025  // modang(psi1, sqrt(-mu)) = atan2(bet1.s() * fabs(alp1.s()),
1026  // bet0.c() * alp1.c());
1027  // assume fbet().fwd(x) = x in this case
1028  u0 = f.fpsi().fwd(psi1.radians());
1029  v0 = f.ftht().fwd(tht1.radians());
1030  delta = (biaxspecial(tg, f.gammax()) ?
1031  atan2(phi1.s() * fabs(alp1.s()), phi0.c() * alp1.c())
1032  - sqrt(f.gammax()) * f.fpsi().df(u0)
1033  : f.fpsi()(u0)) - f.ftht()(v0);
1034  } else if (f.gammax() == 0) {
1035  if (f.kxp2() == 0) {
1036  // meridonal geodesic on biaxial ellipsoid
1037  // N.B. factor of sqrt(f.gammax()) omitted
1038  // Implicitly assume phi1 in [-90, 90] for now
1039 
1040  // phi1, tht1, alp1 are starting conditions.
1041 
1042  // If alp1 = 0/180, alp1.s() == 0, and phi1 != +/-90, phi1.c() != 0,
1043  // then, at baseline equator crossing, phi = phi1.nearest(2U), tht =
1044  // tht0 = tht1, alp = alp1
1045  phi0 = phi1.nearest(2U); // phi0 = 0
1046  // psi1 not used, perhaps it should be?
1047  // psi1 = ang(phi1.s(), phi0.c() * alp1.c() * fabs(phi1.c())));
1048  tht0 = tht1;
1049 
1050  // Othersise at pole, phi1.c() == 0. Define baseline equator crossing
1051  // by alp = alp1.nearest(2U),
1052  // phi = phi1 + cardinal(signbit(alp1.c()) ? -1 : 1)
1053  // tht0 = tht1 +
1054  // (alp1.nearest(2U)-alp1).flipsign(phi1.s())
1055  if (phi1.c() == 0) {
1056  // phi0 = phi1 + ang::cardinal(signbit(alp1.c()) ? -1 : 1);
1057  tht0 += (alp1.nearest(2U)-alp1).flipsign(phi1.s() * Ex);
1058  }
1059  u0 = f.fpsi().fwd((phi1-phi0).radians());
1060  v0 = 0;
1061  delta = -f.ftht()(tht0.radians());
1062  } else {
1063  // N.B. factor of k*kp omitted
1064  // phi0, tht0 are the middle of the initial umbilical segment
1065  if (fabs(phi1.c()) < 8*eps && fabs(tht1.c()) < 8*eps) {
1066  phi0 = phi1.nearest(1U) + ang::cardinal(Nx);
1067  tht0 = tht1.nearest(1U) + ang::cardinal(1);
1068  delta = f.deltashift()/2 - log(fabs(alp1.t()));
1069  } else {
1070  phi0 = phi1.nearest(2U);
1071  tht0 = tht1.nearest(2U);
1072  delta = Nx * f.fpsi()(lamang(phi1 - phi0, tg.kp())) -
1073  f.ftht()(lamang(tht1 - tht0, tg.k()));
1074  }
1075  }
1076  } else {
1077  // gamma = NaN
1078  }
1079  }
1080 
1081  void GeodesicLine3::fline::fics::pos1(bool transpolar,
1082  ang& bet10, ang& omg10, ang& alp10)
1083  const {
1084  bet10 = phi1; omg10 = tht1.flipsign(Ex);
1085  alp10 = alp1;
1086  if (transpolar) {
1087  swap(bet10, omg10);
1088  alp10.reflect(false, false, true);
1089  }
1090  omg10 += ang::cardinal(1);
1091  }
1092 
1093  void GeodesicLine3::fline::fics::setquadrant(const fline& f, unsigned q) {
1094  ang bet1, omg1, alp1x; // TODO: Fix so fics uses tau1
1095  pos1(f.transpolar(), bet1, omg1, alp1x);
1096  alp1x.setquadrant(q);
1097  *this = fics(f, bet1, omg1, alp1x);
1098  }
1099 
1100  GeodesicLine3::gline::gics::gics(const gline& g, const fline::fics& fic) {
1101  if (g.gammax() > 0) {
1102  sig1 = g.gpsi()(fic.u0) + g.gtht()(fic.v0);
1103  } else if (g.gammax() == 0) {
1104  sig1 = g.kxp2() == 0 ? fic.Nx * g.gpsi()(fic.u0) :
1105  fic.Nx * g.gpsi()(lamang(fic.phi1 - fic.phi0, g.kxp())) +
1106  g.gtht()(lamang(fic.tht1 - fic.tht0, g.kx()));
1107  } else {
1108  // gamma = NaN
1109  }
1110  }
1111 
1112  Math::real GeodesicLine3::gline::dist(gics ic, fline::disttx d) const {
1113  real sig2 = gpsi()(d.phiw2) + gtht()(d.thtw2) + d.ind2 * 2*s0;
1114  return (sig2 - ic.sig1) * tg().b();
1115  }
1116 
1117  GeodesicLine3::hfun::hfun(bool distp, real kap, real kapp, real eps, real mu,
1118  const Geodesic3& tg)
1119  : _kap(kap)
1120  , _kapp(kapp)
1121  , _eps(eps)
1122  , _mu(mu)
1123  , _sqrtmu(sqrt(fabs(_mu)))
1124  , _sqrtkap(sqrt(_kap))
1125  , _sqrtkapp(sqrt(_kapp))
1126  , _distp(distp)
1127  , _umb(!tg.biaxial() && _mu == 0)
1128  , _meridr(_kap == 0 && _mu == 0)
1129  , _meridl(_kapp == 0 && _mu == 0)
1130  , _biaxr(biaxspecial(tg, _mu) && _kap == 0)
1131  , _biaxl(biaxspecial(tg, _mu) && _kapp == 0)
1132  {
1133  // mu in [-kap, kapp], eps in (-inf, 1/kap)
1134  if (!_distp) {
1135  if (_meridr || _biaxr) {
1136  // biaxial rotating coordinate
1137  // _kapp == 1, mu < 0 not allowed
1138  _tx = false;
1139  // f multiplied by sqrt(mu)
1140  _fun = TrigfunExt(
1141  [eps = _eps, mu = _mu]
1142  (real tht) -> real
1143  // This is a trivial case f' = 1
1144  { return fthtbiax(tht, eps, mu); },
1145  Math::pi()/2, false);
1146  } else if (_meridl || _biaxl) {
1147  // biaxial librating coordinate
1148  // _kap == 1, mu > 0 not allowed
1149  // DON'T USE tx: _tx = _mu < 0 && -_mu < tg._ellipthresh;
1150  _tx = false;
1151  // f explicitly multiplied by sqrt(-mu) in operator()() for _biaxl
1152  // For _meridl operator() ignores this
1153  _fun = TrigfunExt(
1154  [eps = _eps, mu = _mu]
1155  (real psi) -> real
1156  { return dfpsibiax(sin(psi), cos(psi), eps, mu); },
1157  Math::pi()/2, false);
1158  } else if (_mu > 0) {
1159  _tx = _mu / (_kap + _mu) < tg._ellipthresh;
1160  if (_tx) {
1161  _ell = EllipticFunction(_kap / (_kap + _mu), 0,
1162  _mu / (_kap + _mu), 1);
1163  _fun = TrigfunExt(
1164  [kap = _kap, kapp = _kapp,
1165  eps = _eps, mu = _mu, ell = _ell]
1166  (real u) -> real
1167  { real sn, cn, dn; (void) ell.am(u, sn, cn, dn);
1168  return fup(cn, kap, kapp, eps, mu); },
1169  _ell.K(), false);
1170  } else
1171  _fun = TrigfunExt(
1172  [kap = _kap, kapp = _kapp, eps = _eps, mu = _mu]
1173  (real tht) -> real
1174  { return fthtp(cos(tht), kap, kapp, eps, mu); },
1175  Math::pi()/2, false);
1176  } else if (_mu < 0) {
1177  _tx = -_mu / _kap < tg._ellipthresh;
1178  if (_tx) {
1179  _ell = EllipticFunction((_kap + _mu) / _kap, 0, -_mu / _kap, 1);
1180  _fun = TrigfunExt(
1181  [kap = _kap, kapp = _kapp,
1182  eps = _eps, mu = _mu, ell = _ell]
1183  (real v) -> real
1184  { real sn, cn, dn; (void) ell.am(v, sn, cn, dn);
1185  return fvp(dn, kap, kapp, eps, mu); },
1186  _ell.K(), false);
1187  } else
1188  _fun = TrigfunExt(
1189  [kap = _kap, kapp = _kapp, eps = _eps, mu = _mu]
1190  (real psi) -> real
1191  { return fpsip(sin(psi), cos(psi),
1192  kap, kapp, eps, mu); },
1193  Math::pi()/2, false);
1194  } else if (_umb) {
1195  _tx = _kapp < tg._ellipthresh;
1196  // f multiplied by sqrt(kap*kapp)
1197  // Include scale = 1 in TrigfunExt constructor because this function
1198  // gets added to u.
1199  if (_tx) {
1200  _ell = EllipticFunction(_kap, 0, _kapp, 1);
1201  _fun = TrigfunExt(
1202  [kap = _kap, kapp = _kapp,
1203  eps = _eps, ell = _ell]
1204  (real v) -> real
1205  { real sn, cn, dn; (void) ell.am(v, sn, cn, dn);
1206  return dfvp(cn, dn, kap, kapp, eps); },
1207  2 * _ell.K(), true, 1);
1208  } else
1209  _fun = TrigfunExt(
1210  [
1211 #if defined(_MSC_VER)
1212  // Visual Studio requires the capture of this.
1213  // I can't see why -- probably a bug?
1214  this,
1215 #endif
1216  kap = _kap, kapp = _kapp, eps = _eps]
1217  (real tht) -> real
1218  { return dfp(cos(tht), kap, kapp, eps); },
1219  Math::pi(), true, 1);
1220  } else {
1221  // _mu == NaN
1222  _tx = false;
1223  }
1224  } else { // _distp
1225  if (_meridr || _biaxr) {
1226  // biaxial symmetry coordinate
1227  // _kapp == 1, mu < 0 not allowed
1228  _tx = false;
1229  _fun = TrigfunExt(
1230  [eps = _eps, mu = _mu]
1231  (real tht) -> real
1232  // degenerate f' = 0
1233  { return gthtbiax(tht, eps, mu); },
1234  Math::pi()/2, false);
1235  } else if (_meridl || _biaxl) {
1236  // biaxial non-symmetry coordinate
1237  // _kap == 1, mu > 0 not allowed
1238  // DON'T USE tx: _tx = _mu < 0 && -_mu < tg._ellipthresh;
1239  _tx = false;
1240  _fun = TrigfunExt(
1241  [eps = _eps, mu = _mu]
1242  (real psi) -> real
1243  { return gpsibiax(sin(psi), cos(psi), eps, mu); },
1244  Math::pi()/2, false);
1245  } else if (_mu > 0) {
1246  _tx = _mu / (_kap + _mu) < tg._ellipthresh;
1247  if (_tx) {
1248  _ell = EllipticFunction(_kap / (_kap + _mu), 0,
1249  _mu / (_kap + _mu), 1);
1250  _fun = TrigfunExt(
1251  [kap = _kap, kapp = _kapp,
1252  eps = _eps, mu = _mu, ell = _ell]
1253  (real u) -> real
1254  { real sn, cn, dn; (void) ell.am(u, sn, cn, dn);
1255  return gup(cn, dn, kap, kapp, eps, mu); },
1256  _ell.K());
1257  } else
1258  _fun = TrigfunExt(
1259  [kap = _kap, kapp = _kapp, eps = _eps, mu = _mu]
1260  (real tht) -> real
1261  { return gthtp(cos(tht), kap, kapp, eps, mu); },
1262  Math::pi()/2);
1263  } else if (_mu < 0) {
1264  _tx = -_mu / _kap < tg._ellipthresh;
1265  if (_tx) {
1266  _ell = EllipticFunction((_kap + _mu) / _kap, 0, -_mu / _kap, 1);
1267  _fun = TrigfunExt(
1268  [kap = _kap, kapp = _kapp,
1269  eps = _eps, mu = _mu, ell = _ell]
1270  (real v) -> real
1271  { real sn, cn, dn; (void) ell.am(v, sn, cn, dn);
1272  return gvp(cn, dn, kap, kapp, eps, mu); },
1273  _ell.K());
1274  } else
1275  _fun = TrigfunExt(
1276  [kap = _kap, kapp = _kapp, eps = _eps, mu = _mu]
1277  (real psi) -> real
1278  { return gpsip(sin(psi), cos(psi),
1279  kap, kapp, eps, mu); },
1280  Math::pi()/2);
1281  } else if (_umb) {
1282  _tx = _kapp < tg._ellipthresh;
1283  if (_tx) {
1284  _ell = EllipticFunction(_kap, 0, _kapp, 1);
1285  _fun = TrigfunExt(
1286  [kap = _kap, kapp = _kapp, eps = _eps, ell = _ell]
1287  (real v) -> real
1288  { real sn, cn, dn; (void) ell.am(v, sn, cn, dn);
1289  return g0vp(cn, kap, kapp, eps); },
1290  2*_ell.K(), true);
1291  } else
1292  _fun = TrigfunExt(
1293  [kap = _kap, kapp = _kapp, eps = _eps]
1294  (real tht) -> real
1295  { return g0p(cos(tht), kap, kapp, eps); },
1296  Math::pi(), true);
1297  } else {
1298  // _mu == NaN
1299  _tx = false;
1300  }
1301  }
1302  // N.B. _max < 0 for _umb && eps < 0
1303  _max = _umb ? _fun(_tx ? _ell.K() : Math::pi()/2) :
1304  !_distp ? ( _biaxl ? Math::pi()/2 + _sqrtmu * _fun.Max() :
1305  _fun.Max() ) :
1306  _meridl ? _fun(Math::pi()/2) : _fun.Max();
1307  }
1308 
1309  Math::real GeodesicLine3::hfun::operator()(real u) const {
1310  if (!_distp) {
1311  if (_biaxl)
1312  // This is sqrt(-mu) * f(u)
1313  return modang(u, _sqrtmu) - _sqrtmu * _fun(u);
1314  else if (_meridl)
1315  return 0;
1316  else if (_umb) {
1317  // This is sqrt(kap * kapp) * f(u)
1318  real phi = gd(u, _sqrtkapp);
1319  return u - _fun(_tx ? _ell.F(phi) : phi);
1320  } else
1321  return _fun(u);
1322  } else {
1323  if (_umb) {
1324  real phi = gd(u, _sqrtkapp);
1325  return _fun(_tx ? _ell.F(phi) : phi);
1326  } else
1327  return _fun(u);
1328  }
1329  }
1330 
1331  Math::real GeodesicLine3::hfun::deriv(real u) const {
1332  if (!_distp) {
1333  if (_biaxl)
1334  // This is sqrt(-mu) * f'(u)
1335  return _sqrtmu / (Math::sq(cos(u)) - _mu * Math::sq(sin(u)))
1336  - _sqrtmu * _fun.deriv(u);
1337  else if (_meridl)
1338  return 0;
1339  else if (_umb) {
1340  // This is sqrt(kap * kapp) * f'(u)
1341  real phi = gd(u, _sqrtkapp),
1342  t = _kapp + Math::sq(sinh(u));
1343  // dphi/du = _sqrtkapp * cosh(u) / t;
1344  // for tx w = F(phi, sqrtkap)
1345  // dw/dphi = 1/sqrt(kapp + kap*cos(phi)^2)
1346  // = sqrt(t)/(sqrtkapp *cosh(u))
1347  // N.B. cos(phi)^2 = kapp/t
1348  // dw/du = dw/dphi*dphi/du = 1/sqrt(t)
1349  return 1 - _fun.deriv(_tx ? _ell.F(phi) : phi) /
1350  ( _tx ? sqrt(t) : t / (_sqrtkapp * cosh(u)) );
1351  } else
1352  return _fun.deriv(u);
1353  } else {
1354  if (_umb) {
1355  real phi = gd(u, _sqrtkapp),
1356  t = _kapp + Math::sq(sinh(u));
1357  // See comments in ffun::deriv
1358  return _fun.deriv(_tx ? _ell.F(phi) : phi) /
1359  ( _tx ? sqrt(t) : t / (_sqrtkapp * cosh(u)) );
1360  } else
1361  return _fun.deriv(u);
1362  }
1363  }
1364 
1365  Math::real GeodesicLine3::hfun::root(real z, real u0,
1366  int* countn, int* countb, real tol)
1367  const {
1368  if (!_distp) {
1369  if (!isfinite(z)) return z; // Deals with +/-inf and nan
1370  if (_biaxl) {
1371  // z = 1/sqrt(-mu)
1372  // here f(u) = atan(sqrt(-mu)*tan(u))/sqrt(-mu)-Deltaf(u)
1373  // let fx(u) = sqrt(-mu) * f(u); fun(u) = sqrt(-mu)*Deltaf(u)
1374  // z = fx(u) = modang(u, sqrt(-mu)) - fun(u)
1375  // fun(u) is monotonically increasing/decreasing quasilinear function
1376  // period pi. Combined function is quasilinear
1377  // z = s * u +/- m; period = pi
1378  // Inverting:
1379  // u = z/s +/- m/s; period s*pi
1380  // u = fxinv(z) = modang(z/x, 1/sqrt(-mu)) + funinv(z)
1381  // funinv(z) periodic function of z period (1-s)*pi
1382  real d = Max(),
1383  ua = (z - d) / Slope(),
1384  ub = (z + d) / Slope();
1385  u0 = fmin(ub, fmax(ua, u0));
1386  return Trigfun::root(Trigfun::FFUNROOT,
1387  [this]
1388  (real u) -> pair<real, real>
1389  { return pair<real, real>((*this)(u), deriv(u)); },
1390  z,
1391  u0, ua, ub,
1392  HalfPeriod(), HalfPeriod()/Slope(), 1,
1393  countn, countb, tol);
1394  } else if (_umb) {
1395  real d = fabs(Max())
1396  + 2 * numeric_limits<real>::epsilon() * fmax(real(1), fabs(z)),
1397  ua = z - d,
1398  ub = z + d;
1399  u0 = fmin(ub, fmax(ua, u0));
1400  return Trigfun::root(Trigfun::FFUNROOT,
1401  [this]
1402  (real u) -> pair<real, real>
1403  { return pair<real, real>((*this)(u), deriv(u)); },
1404  z,
1405  u0, ua, ub,
1406  Math::pi()/2, Math::pi()/2, 1,
1407  countn, countb, tol);
1408  } else
1409  return Math::NaN();
1410  } else {
1411  // This function isn't neeed. General inversion mechanisms in Trigfun
1412  // suffice. NO, the trigfun for _umb is not invertible.
1413  if (!(isfinite(z) && _umb))
1414  return Math::NaN(); // Deals with +/-inf and nan
1415  // Now we're dealing with _umb.
1416  if (fabs(z) >= Max())
1417  return copysign(Geodesic3::BigValue(), z);
1418  real ua = -Geodesic3::BigValue(), ub = -ua;
1419  u0 = fmin(ub, fmax(ua, u0));
1420  // Solve z = _fun(_tx ? _ell.F(gd(u)) : gd(u)) for u
1421  return Trigfun::root(Trigfun::GFUNROOT,
1422  [this]
1423  (real u) -> pair<real, real>
1424  { return pair<real, real>((*this)(u), deriv(u)); },
1425  z,
1426  u0, ua, ub,
1427  Math::pi()/2, Math::pi()/2, 1, countn, countb, tol);
1428  }
1429  }
1430 
1431  // Accurate inverse by direct Newton (not using _finv)
1432  Math::real GeodesicLine3::hfun::inv(real z, int* countn, int* countb)
1433  const {
1434  if (!_distp)
1435  return _umb ? root(z, z, countn, countb) :
1436  _biaxl ? root(z, modang(z/Slope(), 1/_sqrtmu), countn, countb) :
1437  _fun.inv1(z, countn, countb);
1438  else { // distp
1439  if (_biaxr) return Math::NaN();
1440  else if (_umb) {
1441  // In limit _eps -> 0
1442  // g(u) = atan(_sqrtkap/_sqrtkapp * tanh(u))
1443  // at u = 0, dg/du = _sqrtkap/_sqrtkapp
1444  // u = inf, g = atan(_sqrtkap/_sqrtkapp)
1445  //
1446  // For _eps finite
1447  // at u = 0, dg/du = _sqrtkap/_sqrtkapp * sqrt(1 - _eps*_kap)
1448  // u = inf, g = _max
1449  // Note: at u = 0
1450  // du/dphi = _sqrtkapp, so
1451  // dg/dphi = _sqrtkap * sqrt(1 - _eps*_kap) -- OK
1452  //
1453  // Approximate g(u) for _eps finite by
1454  // g(u) = _max / atan(_sqrtkap/_sqrtkapp) *
1455  // atan(_sqrtkap/_sqrtkapp *
1456  // tanh(sqrt(1 - _eps*_kap) * atan(_sqrtkap/_sqrtkapp) / _max
1457  // * u))
1458  // Values at +/- inf and slope at origin match.
1459  //
1460  // Solve z = g(u) gives u = u0:
1461  real u0 = atanh(_sqrtkapp/_sqrtkap *
1462  tan(atan(_sqrtkap/_sqrtkapp) / _max * z)) /
1463  (sqrt(1 - _eps*_kap) * atan(_sqrtkap/_sqrtkapp) / _max);
1464  return root(z, u0, countn, countb);
1465  } else
1466  return _fun.inv1(z, countn, countb);
1467  }
1468  }
1469 
1470  // _mu > 0 && !_tx
1471  Math::real GeodesicLine3::hfun::fthtp(real c, real kap, real kapp,
1472  real eps, real mu) {
1473  real c2 = kap * Math::sq(c);
1474  return sqrt((1 - eps * c2) / ((kapp + c2) * (c2 + mu)) );
1475  }
1476  // This is non-negative
1477  Math::real GeodesicLine3::hfun::gthtp(real c, real kap, real kapp,
1478  real eps, real mu) {
1479  real c2 = kap * Math::sq(c);
1480  return c2 * sqrt((1 - eps * c2) / ((kapp + c2) * (c2 + mu)) );
1481  }
1482 
1483  // _mu > 0 && _tx
1484  Math::real GeodesicLine3::hfun::fup(real cn, real kap, real kapp,
1485  real eps, real mu) {
1486  real c2 = kap * Math::sq(cn);
1487  return sqrt( (1 - eps * c2) / ((kapp + c2) * (kap + mu)) );
1488  }
1489  // This is non-negative
1490  Math::real GeodesicLine3::hfun::gup(real cn, real /* dn */,
1491  real kap, real kapp, real eps, real mu) {
1492  real c2 = kap * Math::sq(cn);
1493  return c2 * sqrt( (1 - eps * c2) / ((kapp + c2) * (kap + mu)) );
1494  }
1495 
1496  // _mu == 0 && !_tx
1497  Math::real GeodesicLine3::hfun::dfp(real c, real kap, real kapp, real eps) {
1498  // function dfp = dfpf(phi, kappa, epsilon)
1499  // return derivative of sqrt(kap * kapp) * Delta f
1500  // s = sqrt(1 - kap * sin(phi)^2)
1501  real c2 = kap * Math::sq(c), s = sqrt(kapp + c2);
1502  return eps*kap * sqrt(kapp) * c / (s * (1 + sqrt(1 - eps*c2)));
1503  }
1504  Math::real GeodesicLine3::hfun::g0p(real c, real kap, real kapp, real eps) {
1505  real c2 = kap * Math::sq(c);
1506  return sqrt( kap * (1 - eps * c2) / (kapp + c2) ) * c;
1507  }
1508 
1509  // _mu == 0 && _tx
1510  Math::real GeodesicLine3::hfun::dfvp(real cn, real /* dn */,
1511  real kap, real kapp, real eps) {
1512  // function dfvp = dfvpf(v, kap, eps)
1513  // return derivative of sqrt(kap * kapp) * Delta f_v
1514  return eps*kap * sqrt(kapp) * cn /
1515  (1 + sqrt(1 - eps*kap * Math::sq(cn)));
1516  }
1517  Math::real GeodesicLine3::hfun::g0vp(real cn, real kap, real /* kapp */,
1518  real eps) {
1519  real c2 = kap * Math::sq(cn);
1520  return sqrt( kap * (1 - eps * c2) ) * cn;
1521  }
1522 
1523  // _mu < 0 && !_tx
1524  Math::real GeodesicLine3::hfun::fpsip(real s, real c, real kap, real kapp,
1525  real eps, real mu) {
1526  real c2 = kap * Math::sq(c) - mu * Math::sq(s);
1527  return sqrt( (1 - eps * c2) / ((kapp + c2) * c2) ) ;
1528  }
1529  // This is positive
1530  Math::real GeodesicLine3::hfun::gpsip(real s, real c, real kap, real kapp,
1531  real eps, real mu) {
1532  real c2 = kap * Math::sq(c) - mu * Math::sq(s);
1533  return sqrt(c2 * (1 - eps * c2) / (kapp + c2)) ;
1534  }
1535 
1536  // _mu < 0 && _tx
1537  Math::real GeodesicLine3::hfun::fvp(real dn, real kap, real kapp,
1538  real eps, real /* mu */) {
1539  real c2 = kap * Math::sq(dn);
1540  return sqrt( (1 - eps * c2) / ((kapp + c2) * kap) );
1541  }
1542  // This is positive
1543  Math::real GeodesicLine3::hfun::gvp(real /* cn */, real dn,
1544  real kap, real kapp,
1545  real eps, real /* mu */) {
1546  real dn2 = Math::sq(dn), c2 = kap * dn2;
1547  return dn2 * sqrt( kap * (1 - eps * c2) / (kapp + c2) );
1548  }
1549 
1550  // biaxial variants for kap = 0, kapp = 1, mu > 0
1551  Math::real
1552  GeodesicLine3::hfun::fthtbiax(real /* tht */, real /* eps */,
1553  real /* mu */) {
1554  // Multiply by f functions by sqrt(abs(mu))
1555  // return 1 / sqrt(mu);
1556  return 1;
1557  }
1558  Math::real
1559  GeodesicLine3::hfun::gthtbiax(real /* tht */, real /* eps */,
1560  real /* mu */) {
1561  return 0;
1562  }
1563 
1564  // biaxial variants for kap = 1, kapp = 0, mu <= 0, !_tx
1565  Math::real
1566  GeodesicLine3::hfun::dfpsibiax(real s, real c, real eps, real mu) {
1567  real c2 = Math::sq(c) - mu * Math::sq(s);
1568  // f functions are multiplied by sqrt(abs(mu)) but don't include this
1569  // factor here; instead include it in operator()(). etc. This was we can
1570  // still use this function in the limit mu -> 0 to determine the conjugate
1571  // point on a meridian.
1572  return eps / (1 + sqrt(1 - eps * c2));
1573  }
1574  Math::real
1575  GeodesicLine3::hfun::gpsibiax(real s, real c, real eps, real mu) {
1576  real c2 = Math::sq(c) - mu * Math::sq(s);
1577  // return gpsip(s, c, 1, 0, eps, mu) but with the factor c2 canceled
1578  return sqrt(1 - eps * c2);
1579  }
1580 
1581  } // namespace Triaxial
1582 } // namespace GeographicLib
GeographicLib::TrigfunExt
A function defined by its derivative and its inverse.
Definition: Trigfun.hpp:498
GeographicLib::Math::pi
static T pi()
Definition: Math.hpp:187
GeographicLib::AngleT::radians0
T radians0() const
Definition: Angle.hpp:549
GeographicLib::Math::infinity
static T infinity()
Definition: Math.cpp:308
GeographicLib::AngleT::radians
static AngleT radians(T rad)
Definition: Angle.hpp:534
GeographicLib::Triaxial::GeodesicLine3::fline::fics::pos1
void pos1(bool transpolar, ang &bet10, ang &omg10, ang &alp10) const
Definition: GeodesicLine3.cpp:1081
GeodesicLine3.hpp
Header for GeographicLib::Triaxial::GeodesicLine3 class.
GeographicLib::AngleT::reflect
AngleT & reflect(bool flips, bool flipc=false, bool swapp=false)
Definition: Angle.hpp:696
GeographicLib::AngleT::nearest
AngleT nearest(unsigned ind=0U) const
Definition: Angle.hpp:743
GeographicLib::AngleT::ncardinal
T ncardinal() const
Definition: Angle.hpp:572
GeographicLib::Math::clamp
static T clamp(T x, T a, T b)
Definition: Math.cpp:296
GeographicLib::Triaxial::GeodesicLine3::pos1
void pos1(Angle &bet1, Angle &omg1, Angle &alp1) const
Definition: GeodesicLine3.cpp:54
GeographicLib::AngleT< Math::real >::cardinal
static AngleT cardinal(Math::real q)
Definition: Angle.hpp:721
GeographicLib::AngleT::setn
AngleT & setn(T n=0)
Definition: Angle.hpp:662
GeographicLib::AngleT< Math::real >::radians
Math::real radians() const
Definition: Angle.hpp:542
GeographicLib::AngleT::setquadrant
AngleT & setquadrant(unsigned q)
Definition: Angle.hpp:682
GeographicLib::Math::sq
static T sq(T x)
Definition: Math.hpp:209
GeographicLib::Triaxial::GeodesicLine3
The direct geodesic problem for a triaxial ellipsoid.
Definition: GeodesicLine3.hpp:54
GeographicLib::Triaxial::Ellipsoid3::AngNorm
static bool AngNorm(Angle &bet, Angle &omg, Angle &alp, bool alt=false)
Definition: Ellipsoid3.hpp:247
GeographicLib
Namespace for GeographicLib.
Definition: Accumulator.cpp:12
GeographicLib::Math::NaN
static T NaN()
Definition: Math.cpp:301
GeographicLib::AngleT::flipsign
AngleT flipsign(T mult) const
Definition: Angle.hpp:705
real
GeographicLib::Math::real real
Definition: Geod3Solve.cpp:25
std::swap
void swap(GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &a, GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &b)
Definition: NearestNeighbor.hpp:819
GeographicLib::Triaxial::Geodesic3
The solution of the geodesic problem for a triaxial ellipsoid.
Definition: Geodesic3.hpp:65
GeographicLib::GeographicErr
Exception handling for GeographicLib.
Definition: Constants.hpp:344
GeographicLib::AngleT::base
AngleT base() const
Definition: Angle.hpp:642
ang
GeographicLib::Angle ang
Definition: Geod3Solve.cpp:26
GeographicLib::AngleT::rebase
AngleT rebase(const AngleT &c) const
Definition: Angle.hpp:647
GeographicLib::AngleT::round
AngleT & round()
Definition: Angle.hpp:636
GeographicLib::Triaxial::GeodesicLine3::Position
void Position(real s12, Angle &bet2, Angle &omg2, Angle &alp2) const
Definition: GeodesicLine3.cpp:71