GeographicLib  2.1.2
JacobiConformal.hpp
Go to the documentation of this file.
1 /**
2  * \file JacobiConformal.hpp
3  * \brief Header for GeographicLib::JacobiConformal class
4  *
5  * \note This is just sample code. It is not part of GeographicLib
6  * itself.
7  *
8  * Copyright (c) Charles Karney (2014-2020) <charles@karney.com> and licensed
9  * under the MIT/X11 License. For more information, see
10  * https://geographiclib.sourceforge.io/
11  **********************************************************************/
12 
13 #if !defined(GEOGRAPHICLIB_JACOBICONFORMAL_HPP)
14 #define GEOGRAPHICLIB_JACOBICONFORMAL_HPP 1
15 
17 
18 namespace GeographicLib {
19  /**
20  * \brief Jacobi's conformal projection of a triaxial ellipsoid
21  *
22  * <b>NOTE:</b> This is just sample code. It is not part of GeographicLib
23  * itself.
24  *
25  * This is a conformal projection of the ellipsoid to a plane in which
26  * the grid lines are straight; see Jacobi,
27  * <a href="https://books.google.com/books?id=ryEOAAAAQAAJ&pg=PA212">
28  * Vorlesungen &uuml;ber Dynamik, &sect;28</a>. The constructor takes the
29  * semi-axes of the ellipsoid (which must be in order). Member functions map
30  * the ellipsoidal coordinates &omega; and &beta; separately to \e x and \e
31  * y. Jacobi's coordinates have been multiplied by
32  * (<i>a</i><sup>2</sup>&minus;<i>c</i><sup>2</sup>)<sup>1/2</sup> /
33  * (2<i>b</i>) so that the customary results are returned in the cases of
34  * a sphere or an ellipsoid of revolution.
35  *
36  * The ellipsoid is oriented so that the large principal ellipse, \f$Z=0\f$,
37  * is the equator, \f$\beta=0\f$, while the small principal ellipse,
38  * \f$Y=0\f$, is the prime meridian, \f$\omega=0\f$. The four umbilic
39  * points, \f$\left|\omega\right| = \left|\beta\right| = \frac12\pi\f$, lie
40  * on middle principal ellipse in the plane \f$X=0\f$.
41  *
42  * For more information on this projection, see \ref jacobi.
43  **********************************************************************/
45  typedef Math::real real;
46  real _a, _b, _c, _ab2, _bc2, _ac2;
47  EllipticFunction _ex, _ey;
48  static void norm(real& x, real& y) {
49  using std::hypot;
50  real z = hypot(x, y); x /= z; y /= z;
51  }
52  public:
53  /**
54  * Constructor for a trixial ellipsoid with semi-axes.
55  *
56  * @param[in] a the largest semi-axis.
57  * @param[in] b the middle semi-axis.
58  * @param[in] c the smallest semi-axis.
59  *
60  * The semi-axes must satisfy \e a &ge; \e b &ge; \e c > 0 and \e a >
61  * \e c. This form of the constructor cannot be used to specify a
62  * sphere (use the next constructor).
63  **********************************************************************/
64  JacobiConformal(real a, real b, real c)
65  : _a(a), _b(b), _c(c)
66  , _ab2((_a - _b) * (_a + _b))
67  , _bc2((_b - _c) * (_b + _c))
68  , _ac2((_a - _c) * (_a + _c))
69  , _ex(_ab2 / _ac2 * Math::sq(_c / _b), -_ab2 / Math::sq(_b),
70  _bc2 / _ac2 * Math::sq(_a / _b), Math::sq(_a / _b))
71  , _ey(_bc2 / _ac2 * Math::sq(_a / _b), +_bc2 / Math::sq(_b),
72  _ab2 / _ac2 * Math::sq(_c / _b), Math::sq(_c / _b))
73  {
74  using std::isfinite;
75  if (!(isfinite(_a) && _a >= _b && _b >= _c && _c > 0))
76  throw GeographicErr("JacobiConformal: axes are not in order");
77  if (!(_a > _c))
78  throw GeographicErr
79  ("JacobiConformal: use alternate constructor for sphere");
80  }
81  /**
82  * Alternate constructor for a triaxial ellipsoid.
83  *
84  * @param[in] a the largest semi-axis.
85  * @param[in] b the middle semi-axis.
86  * @param[in] c the smallest semi-axis.
87  * @param[in] ab the relative magnitude of \e a &minus; \e b.
88  * @param[in] bc the relative magnitude of \e b &minus; \e c.
89  *
90  * This form can be used to specify a sphere. The semi-axes must
91  * satisfy \e a &ge; \e b &ge; c > 0. The ratio \e ab : \e bc must equal
92  * (<i>a</i>&minus;<i>b</i>) : (<i>b</i>&minus;<i>c</i>) with \e ab
93  * &ge; 0, \e bc &ge; 0, and \e ab + \e bc > 0.
94  **********************************************************************/
95  JacobiConformal(real a, real b, real c, real ab, real bc)
96  : _a(a), _b(b), _c(c)
97  , _ab2(ab * (_a + _b))
98  , _bc2(bc * (_b + _c))
99  , _ac2(_ab2 + _bc2)
100  , _ex(_ab2 / _ac2 * Math::sq(_c / _b),
101  -(_a - _b) * (_a + _b) / Math::sq(_b),
102  _bc2 / _ac2 * Math::sq(_a / _b), Math::sq(_a / _b))
103  , _ey(_bc2 / _ac2 * Math::sq(_a / _b),
104  +(_b - _c) * (_b + _c) / Math::sq(_b),
105  _ab2 / _ac2 * Math::sq(_c / _b), Math::sq(_c / _b))
106  {
107  using std::isfinite;
108  if (!(isfinite(_a) && _a >= _b && _b >= _c && _c > 0 &&
109  ab >= 0 && bc >= 0))
110  throw GeographicErr("JacobiConformal: axes are not in order");
111  if (!(ab + bc > 0 && isfinite(_ac2)))
112  throw GeographicErr("JacobiConformal: ab + bc must be positive");
113  }
114  /**
115  * @return the quadrant length in the \e x direction.
116  **********************************************************************/
117  Math::real x() const { return Math::sq(_a / _b) * _ex.Pi(); }
118  /**
119  * The \e x projection.
120  *
121  * @param[in] somg sin(&omega;).
122  * @param[in] comg cos(&omega;).
123  * @return \e x.
124  **********************************************************************/
125  Math::real x(real somg, real comg) const {
126  real somg1 = _b * somg, comg1 = _a * comg; norm(somg1, comg1);
127  return Math::sq(_a / _b)
128  * _ex.Pi(somg1, comg1, _ex.Delta(somg1, comg1));
129  }
130  /**
131  * The \e x projection.
132  *
133  * @param[in] omg &omega; (in degrees).
134  * @return \e x (in degrees).
135  *
136  * &omega; must be in [&minus;180&deg;, 180&deg;].
137  **********************************************************************/
138  Math::real x(real omg) const {
139  real somg, comg;
140  Math::sincosd(omg, somg, comg);
141  return x(somg, comg) / Math::degree();
142  }
143  /**
144  * @return the quadrant length in the \e y direction.
145  **********************************************************************/
146  Math::real y() const { return Math::sq(_c / _b) * _ey.Pi(); }
147  /**
148  * The \e y projection.
149  *
150  * @param[in] sbet sin(&beta;).
151  * @param[in] cbet cos(&beta;).
152  * @return \e y.
153  **********************************************************************/
154  Math::real y(real sbet, real cbet) const {
155  real sbet1 = _b * sbet, cbet1 = _c * cbet; norm(sbet1, cbet1);
156  return Math::sq(_c / _b)
157  * _ey.Pi(sbet1, cbet1, _ey.Delta(sbet1, cbet1));
158  }
159  /**
160  * The \e y projection.
161  *
162  * @param[in] bet &beta; (in degrees).
163  * @return \e y (in degrees).
164  *
165  * &beta; must be in (&minus;180&deg;, 180&deg;].
166  **********************************************************************/
167  Math::real y(real bet) const {
168  real sbet, cbet;
169  Math::sincosd(bet, sbet, cbet);
170  return y(sbet, cbet) / Math::degree();
171  }
172  };
173 
174 } // namespace GeographicLib
175 
176 #endif // GEOGRAPHICLIB_JACOBICONFORMAL_HPP
Header for GeographicLib::EllipticFunction class.
Elliptic integrals and functions.
Math::real Delta(real sn, real cn) const
Exception handling for GeographicLib.
Definition: Constants.hpp:316
Jacobi's conformal projection of a triaxial ellipsoid.
Math::real y(real bet) const
Math::real x(real omg) const
JacobiConformal(real a, real b, real c)
Math::real y(real sbet, real cbet) const
Math::real x(real somg, real comg) const
JacobiConformal(real a, real b, real c, real ab, real bc)
Mathematical functions needed by GeographicLib.
Definition: Math.hpp:76
static T degree()
Definition: Math.hpp:200
static void sincosd(T x, T &sinx, T &cosx)
Definition: Math.cpp:106
static T sq(T x)
Definition: Math.hpp:212
Namespace for GeographicLib.
Definition: Accumulator.cpp:12