1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
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
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
/**
* \file CircularEngine.hpp
* \brief Header for GeographicLib::CircularEngine class
*
* Copyright (c) Charles Karney (2011-2015) <karney@alum.mit.edu> and licensed
* under the MIT/X11 License. For more information, see
* https://geographiclib.sourceforge.io/
**********************************************************************/
#if !defined(GEOGRAPHICLIB_CIRCULARENGINE_HPP)
#define GEOGRAPHICLIB_CIRCULARENGINE_HPP 1
#include <vector>
#include <GeographicLib/Constants.hpp>
#include <GeographicLib/SphericalEngine.hpp>
#if defined(_MSC_VER)
// Squelch warnings about dll vs vector
# pragma warning (push)
# pragma warning (disable: 4251)
#endif
namespace GeographicLib {
/**
* \brief Spherical harmonic sums for a circle
*
* The class is a companion to SphericalEngine. If the results of a
* spherical harmonic sum are needed for several points on a circle of
* constant latitude \e lat and height \e h, then SphericalEngine::Circle can
* compute the inner sum, which is independent of longitude \e lon, and
* produce a CircularEngine object. CircularEngine::operator()() can
* then be used to perform the outer sum for particular vales of \e lon.
* This can lead to substantial improvements in computational speed for high
* degree sum (approximately by a factor of \e N / 2 where \e N is the
* maximum degree).
*
* CircularEngine is tightly linked to the internals of SphericalEngine. For
* that reason, the constructor for this class is private. Use
* SphericalHarmonic::Circle, SphericalHarmonic1::Circle, and
* SphericalHarmonic2::Circle to create instances of this class.
*
* CircularEngine stores the coefficients needed to allow the summation over
* order to be performed in 2 or 6 vectors of length \e M + 1 (depending on
* whether gradients are to be calculated). For this reason the constructor
* may throw a std::bad_alloc exception.
*
* Example of use:
* \include example-CircularEngine.cpp
**********************************************************************/
class GEOGRAPHICLIB_EXPORT CircularEngine {
private:
typedef Math::real real;
enum normalization {
FULL = SphericalEngine::FULL,
SCHMIDT = SphericalEngine::SCHMIDT,
};
int _mM;
bool _gradp;
unsigned _norm;
real _a, _r, _u, _t;
std::vector<real> _wc, _ws, _wrc, _wrs, _wtc, _wts;
real _q, _uq, _uq2;
Math::real Value(bool gradp, real sl, real cl,
real& gradx, real& grady, real& gradz) const;
friend class SphericalEngine;
CircularEngine(int M, bool gradp, unsigned norm,
real a, real r, real u, real t)
: _mM(M)
, _gradp(gradp)
, _norm(norm)
, _a(a)
, _r(r)
, _u(u)
, _t(t)
, _wc(std::vector<real>(_mM + 1, 0))
, _ws(std::vector<real>(_mM + 1, 0))
, _wrc(std::vector<real>(_gradp ? _mM + 1 : 0, 0))
, _wrs(std::vector<real>(_gradp ? _mM + 1 : 0, 0))
, _wtc(std::vector<real>(_gradp ? _mM + 1 : 0, 0))
, _wts(std::vector<real>(_gradp ? _mM + 1 : 0, 0))
{
_q = _a / _r;
_uq = _u * _q;
_uq2 = Math::sq(_uq);
}
void SetCoeff(int m, real wc, real ws)
{ _wc[m] = wc; _ws[m] = ws; }
void SetCoeff(int m, real wc, real ws,
real wrc, real wrs, real wtc, real wts) {
_wc[m] = wc; _ws[m] = ws;
if (_gradp) {
_wrc[m] = wrc; _wrs[m] = wrs;
_wtc[m] = wtc; _wts[m] = wts;
}
}
public:
/**
* A default constructor. CircularEngine::operator()() on the resulting
* object returns zero. The resulting object can be assigned to the result
* of SphericalHarmonic::Circle.
**********************************************************************/
CircularEngine()
: _mM(-1)
, _gradp(true)
, _u(0)
, _t(1)
{}
/**
* Evaluate the sum for a particular longitude given in terms of its
* sine and cosine.
*
* @param[in] sinlon the sine of the longitude.
* @param[in] coslon the cosine of the longitude.
* @return \e V the value of the sum.
*
* The arguments must satisfy <i>sinlon</i><sup>2</sup> +
* <i>coslon</i><sup>2</sup> = 1.
**********************************************************************/
Math::real operator()(real sinlon, real coslon) const {
real dummy;
return Value(false, sinlon, coslon, dummy, dummy, dummy);
}
/**
* Evaluate the sum for a particular longitude.
*
* @param[in] lon the longitude (degrees).
* @return \e V the value of the sum.
**********************************************************************/
Math::real operator()(real lon) const {
real sinlon, coslon;
Math::sincosd(lon, sinlon, coslon);
return (*this)(sinlon, coslon);
}
/**
* Evaluate the sum and its gradient for a particular longitude given in
* terms of its sine and cosine.
*
* @param[in] sinlon the sine of the longitude.
* @param[in] coslon the cosine of the longitude.
* @param[out] gradx \e x component of the gradient.
* @param[out] grady \e y component of the gradient.
* @param[out] gradz \e z component of the gradient.
* @return \e V the value of the sum.
*
* The gradients will only be computed if the CircularEngine object was
* created with this capability (e.g., via \e gradp = true in
* SphericalHarmonic::Circle). If not, \e gradx, etc., will not be
* touched. The arguments must satisfy <i>sinlon</i><sup>2</sup> +
* <i>coslon</i><sup>2</sup> = 1.
**********************************************************************/
Math::real operator()(real sinlon, real coslon,
real& gradx, real& grady, real& gradz) const {
return Value(true, sinlon, coslon, gradx, grady, gradz);
}
/**
* Evaluate the sum and its gradient for a particular longitude.
*
* @param[in] lon the longitude (degrees).
* @param[out] gradx \e x component of the gradient.
* @param[out] grady \e y component of the gradient.
* @param[out] gradz \e z component of the gradient.
* @return \e V the value of the sum.
*
* The gradients will only be computed if the CircularEngine object was
* created with this capability (e.g., via \e gradp = true in
* SphericalHarmonic::Circle). If not, \e gradx, etc., will not be
* touched.
**********************************************************************/
Math::real operator()(real lon,
real& gradx, real& grady, real& gradz) const {
real sinlon, coslon;
Math::sincosd(lon, sinlon, coslon);
return (*this)(sinlon, coslon, gradx, grady, gradz);
}
};
} // namespace GeographicLib
#if defined(_MSC_VER)
# pragma warning (pop)
#endif
#endif // GEOGRAPHICLIB_CIRCULARENGINE_HPP