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