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