This file is indexed.

/usr/include/GeographicLib/SphericalHarmonic.hpp is in libgeographiclib-dev 1.21-1ubuntu1.

This file is owned by root:root, with mode 0o644.

The actual contents of the file can be viewed below.

  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
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
/**
 * \file SphericalHarmonic.hpp
 * \brief Header for GeographicLib::SphericalHarmonic class
 *
 * Copyright (c) Charles Karney (2011) <charles@karney.com> and licensed under
 * the MIT/X11 License.  For more information, see
 * http://geographiclib.sourceforge.net/
 **********************************************************************/

#if !defined(GEOGRAPHICLIB_SPHERICALHARMONIC_HPP)
#define GEOGRAPHICLIB_SPHERICALHARMONIC_HPP \
  "$Id: 6fa804c46efd01670cfb7835dd022791b60d2942 $"

#include <vector>
#include <GeographicLib/Constants.hpp>
#include <GeographicLib/SphericalEngine.hpp>
#include <GeographicLib/CircularEngine.hpp>
#include <GeographicLib/Geocentric.hpp>

namespace GeographicLib {

  /**
   * \brief Spherical Harmonic series
   *
   * This class evaluates the spherical harmonic sum \verbatim
 V(x, y, z) = sum(n = 0..N)[ q^(n+1) * sum(m = 0..n)[
   (C[n,m] * cos(m*lambda) + S[n,m] * sin(m*lambda)) *
   P[n,m](cos(theta)) ] ]
\endverbatim
   * where
   * - <i>p</i><sup>2</sup> = <i>x</i><sup>2</sup> + <i>y</i><sup>2</sup>,
   * - <i>r</i><sup>2</sup> = <i>p</i><sup>2</sup> + <i>z</i><sup>2</sup>,
   * - \e q = <i>a</i>/<i>r</i>,
   * - \e theta = atan2(\e p, \e z) = the spherical \e colatitude,
   * - \e lambda = atan2(\e y, \e x) = the longitude.
   * - P<sub>\e nm</sub>(\e t) is the associated Legendre polynomial of degree
   *   \e n and order \e m.
   *
   * Two normalizations are supported for P<sub>\e nm</sub>
   * - fully normalized denoted by SphericalHarmonic::FULL.
   * - Schmidt semi-normalized denoted by SphericalHarmonic::SCHMIDT.
   *
   * Clenshaw summation is used for the sums over both \e n and \e m.  This
   * allows the computation to be carried out without the need for any
   * temporary arrays.  See SphericalEngine.cpp for more information on the
   * implementation.
   *
   * References:
   * - C. W. Clenshaw, A note on the summation of Chebyshev series,
   *   %Math. Tables Aids Comput. 9(51), 118-120 (1955).
   * - R. E. Deakin, Derivatives of the earth's potentials, Geomatics
   *   Research Australasia 68, 31-60, (June 1998).
   * - W. A. Heiskanen and H. Moritz, Physical Geodesy, (Freeman, San
   *   Francisco, 1967).  (See Sec. 1-14, for a definition of Pbar.)
   * - S. A. Holmes and W. E. Featherstone, A unified approach to the
   *   Clenshaw summation and the recursive computation of very high degree
   *   and order normalised associated Legendre functions, J. Geod. 76(5),
   *   279-299 (2002).
   * - C. C. Tscherning and K. Poder, Some geodetic applications of Clenshaw
   *   summation, Boll. Geod. Sci. Aff. 41(4), 349-375 (1982).
   *
   * Example of use:
   * \include example-SphericalHarmonic.cpp
   **********************************************************************/

  class GEOGRAPHIC_EXPORT SphericalHarmonic {
  public:
    /**
     * Supported normalizations for the associated Legendre polynomials.
     **********************************************************************/
    enum normalization {
      /**
       * Fully normalized associated Legendre polynomials.
       *
       * These are defined by <i>P</i><sub><i>nm</i></sub><sup>full</sup>(\e z)
       * = (-1)<sup><i>m</i></sup> sqrt(\e k (2\e n + 1) (\e n - \e m)! / (\e n
       * + \e m)!) <b>P</b><sub><i>n</i></sub><sup><i>m</i></sup>(\e z), where
       * <b>P</b><sub><i>n</i></sub><sup><i>m</i></sup>(\e z) is Ferrers
       * function (also known as the Legendre function on the cut or the
       * associated Legendre polynomial) http://dlmf.nist.gov/14.7.E10 and \e k
       * = 1 for \e m = 0 and \e k = 2 otherwise.
       *
       * The mean squared value of
       * <i>P</i><sub><i>nm</i></sub><sup>full</sup>(cos \e theta) cos(\e m \e
       * lambda) and <i>P</i><sub><i>nm</i></sub><sup>full</sup>(cos \e theta)
       * sin(\e m \e lambda) over the sphere is 1.
       *
       * @hideinitializer
       **********************************************************************/
      FULL = SphericalEngine::FULL,
      /**
       * Schmidt semi-normalized associated Legendre polynomials.
       *
       * These are defined by <i>P</i><sub><i>nm</i></sub><sup>schmidt</sup>(\e
       * z) = (-1)<sup><i>m</i></sup> sqrt(\e k (\e n - \e m)! / (\e n + \e
       * m)!)  <b>P</b><sub><i>n</i></sub><sup><i>m</i></sup>(\e z), where
       * <b>P</b><sub><i>n</i></sub><sup><i>m</i></sup>(\e z) is Ferrers
       * function (also known as the Legendre function on the cut or the
       * associated Legendre polynomial) http://dlmf.nist.gov/14.7.E10 and \e k
       * = 1 for \e m = 0 and \e k = 2 otherwise.
       *
       * The mean squared value of
       * <i>P</i><sub><i>nm</i></sub><sup>schmidt</sup>(cos \e theta) cos(\e m
       * \e lambda) and <i>P</i><sub><i>nm</i></sub><sup>schmidt</sup>(cos \e
       * theta) sin(\e m \e lambda) over the sphere is 1/(2\e n + 1).
       *
       * @hideinitializer
       **********************************************************************/
      SCHMIDT = SphericalEngine::SCHMIDT,
      /// \cond SKIP
      // These are deprecated...
      full = FULL,
      schmidt = SCHMIDT,
      /// \endcond
    };

  private:
    typedef Math::real real;
    SphericalEngine::coeff _c[1];
    real _a;
    unsigned _norm;

  public:
    /**
     * Constructor with a full set of coefficients specified.
     *
     * @param[in] C the coefficients \e C<sub>\e nm</sub>.
     * @param[in] S the coefficients \e S<sub>\e nm</sub>.
     * @param[in] N the maximum degree and order of the sum
     * @param[in] a the reference radius appearing in the definition of the
     *   sum.
     * @param[in] norm the normalization for the associated Legendre
     *   polynomials, either SphericalHarmonic::full (the default) or
     *   SphericalHarmonic::schmidt.
     *
     * The coefficients \e C<sub>\e nm</sub> and \e S<sub>\e nm</sub> are
     * stored in the one-dimensional vectors \e C and \e S which must contain
     * (\e N + 1)(\e N + 2)/2 and N (\e N + 1)/2 elements, respectively, stored
     * in "column-major" order.  Thus for \e N = 3, the order would be:
     * <i>C</i><sub>00</sub>,
     * <i>C</i><sub>10</sub>,
     * <i>C</i><sub>20</sub>,
     * <i>C</i><sub>30</sub>,
     * <i>C</i><sub>11</sub>,
     * <i>C</i><sub>21</sub>,
     * <i>C</i><sub>31</sub>,
     * <i>C</i><sub>22</sub>,
     * <i>C</i><sub>32</sub>,
     * <i>C</i><sub>33</sub>.
     * In general the (\e n,\e m) element is at index \e m*\e N - \e m*(\e m -
     * 1)/2 + \e n.  The layout of \e S is the same except that the first
     * column is omitted (since the \e m = 0 terms never contribute to the sum)
     * and the 0th element is <i>S</i><sub>11</sub>
     *
     * The class stores <i>pointers</i> to the first elements of \e C and \e S.
     * These arrays should not be altered or destroyed during the lifetime of a
     * SphericalHarmonic object.
     **********************************************************************/
    SphericalHarmonic(const std::vector<real>& C,
                      const std::vector<real>& S,
                      int N, real a, unsigned norm = FULL)
      : _a(a)
      , _norm(norm)
    { _c[0] = SphericalEngine::coeff(C, S, N); }

    /**
     * Constructor with a subset of coefficients specified.
     *
     * @param[in] C the coefficients \e C<sub>\e nm</sub>.
     * @param[in] S the coefficients \e S<sub>\e nm</sub>.
     * @param[in] N the degree used to determine the layout of \e C and \e S.
     * @param[in] nmx the maximum degree used in the sum.  The sum over \e n is
     *   from 0 thru \e nmx.
     * @param[in] mmx the maximum order used in the sum.  The sum over \e m is
     *   from 0 thru min(\e n, \e mmx).
     * @param[in] a the reference radius appearing in the definition of the
     *   sum.
     * @param[in] norm the normalization for the associated Legendre
     *   polynomials, either SphericalHarmonic::FULL (the default) or
     *   SphericalHarmonic::SCHMIDT.
     *
     * The class stores <i>pointers</i> to the first elements of \e C and \e S.
     * These arrays should not be altered or destroyed during the lifetime of a
     * SphericalHarmonic object.
     **********************************************************************/
    SphericalHarmonic(const std::vector<real>& C,
                      const std::vector<real>& S,
                      int N, int nmx, int mmx,
                      real a, unsigned norm = FULL)
      : _a(a)
      , _norm(norm)
    { _c[0] = SphericalEngine::coeff(C, S, N, nmx, mmx); }

    /**
     * A default constructor so that the object can be created when the
     * constructor for another object is initialized.  This default object can
     * then be reset with the default copy assignment operator.
     **********************************************************************/
    SphericalHarmonic() {}

    /**
     * Compute the spherical harmonic sum.
     *
     * @param[in] x cartesian coordinate.
     * @param[in] y cartesian coordinate.
     * @param[in] z cartesian coordinate.
     * @return \e V the spherical harmonic sum.
     *
     * This routine requires constant memory and thus never throws an
     * exception.
     **********************************************************************/
    Math::real operator()(real x, real y, real z) const throw() {
      real f[] = {1};
      real v = 0;
      real dummy;
      switch (_norm) {
      case FULL:
        v = SphericalEngine::Value<false, SphericalEngine::FULL, 1>
          (_c, f, x, y, z, _a, dummy, dummy, dummy);
        break;
      case SCHMIDT:
        v = SphericalEngine::Value<false, SphericalEngine::SCHMIDT, 1>
          (_c, f, x, y, z, _a, dummy, dummy, dummy);
        break;
      }
      return v;
    }

    /**
     * Compute a spherical harmonic sum and its gradient.
     *
     * @param[in] x cartesian coordinate.
     * @param[in] y cartesian coordinate.
     * @param[in] z cartesian coordinate.
     * @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 spherical harmonic sum.
     *
     * This is the same as the previous function, except that the components of
     * the gradients of the sum in the \e x, \e y, and \e z directions are
     * computed.  This routine requires constant memory and thus never throws
     * an exception.
     **********************************************************************/
    Math::real operator()(real x, real y, real z,
                          real& gradx, real& grady, real& gradz) const throw() {
      real f[] = {1};
      real v = 0;
      switch (_norm) {
      case FULL:
        v = SphericalEngine::Value<true, SphericalEngine::FULL, 1>
          (_c, f, x, y, z, _a, gradx, grady, gradz);
        break;
      case SCHMIDT:
        v = SphericalEngine::Value<true, SphericalEngine::SCHMIDT, 1>
          (_c, f, x, y, z, _a, gradx, grady, gradz);
        break;
      }
      return v;
    }

    /**
     * Create a CircularEngine to allow the efficient evaluation of several
     * points on a circle of latitude.
     *
     * @param[in] p the radius of the circle.
     * @param[in] z the height of the circle above the equatorial plane.
     * @param[in] gradp if true the returned object will be able to compute the
     *   gradient of the sum.
     * @return the CircularEngine object.
     *
     * SphericalHarmonic::operator()() exchanges the order of the sums in the
     * definition, i.e., sum(n = 0..N)[sum(m = 0..n)[...]] becomes sum(m =
     * 0..N)[sum(n = m..N)[...]].  SphericalHarmonic::Circle performs the inner
     * sum over degree \e n (which entails about <i>N</i><sup>2</sup>
     * operations).  Calling CircularEngine::operator()() on the returned
     * object performs the outer sum over the order \e m (about \e N
     * operations).  This routine may throw a bad_alloc exception in the
     * CircularEngine constructor.
     *
     * Here's an example of computing the spherical sum at a sequence of
     * longitudes without using a CircularEngine object
     \code
  SphericalHarmonic h(...);     // Create the SphericalHarmonic object
  double r = 2, lat = 33, lon0 = 44, dlon = 0.01;
  double
    phi = lat * Math::degree<double>(),
    z = r * sin(phi), p = r * cos(phi);
  for (int i = 0; i <= 100; ++i) {
    real
      lon = lon0 + i * dlon,
      lam = lon * Math::degree<double>();
    std::cout << lon << " " << h(p * cos(lam), p * sin(lam), z) << "\n";
  }
     \endcode
     * Here is the same calculation done using a CircularEngine object.  This
     * will be about <i>N</i>/2 times faster.
     \code
  SphericalHarmonic h(...);     // Create the SphericalHarmonic object
  double r = 2, lat = 33, lon0 = 44, dlon = 0.01;
  double
    phi = lat * Math::degree<double>(),
    z = r * sin(phi), p = r * cos(phi);
  CircularEngine c(h(p, z, false)); // Create the CircularEngine object
  for (int i = 0; i <= 100; ++i) {
    real
      lon = lon0 + i * dlon;
    std::cout << lon << " " << c(lon) << "\n";
  }
     \endcode
     **********************************************************************/
    CircularEngine Circle(real p, real z, bool gradp) const {
      real f[] = {1};
      switch (_norm) {
      case FULL:
        return gradp ?
          SphericalEngine::Circle<true, SphericalEngine::FULL, 1>
          (_c, f, p, z, _a) :
          SphericalEngine::Circle<false, SphericalEngine::FULL, 1>
          (_c, f, p, z, _a);
        break;
      case SCHMIDT:
      default:                  // To avoid compiler warnings
        return gradp ?
          SphericalEngine::Circle<true, SphericalEngine::SCHMIDT, 1>
          (_c, f, p, z, _a) :
          SphericalEngine::Circle<false, SphericalEngine::SCHMIDT, 1>
          (_c, f, p, z, _a);
        break;
      }
    }

    /**
     * @return the zeroth SphericalEngine::coeff object.
     **********************************************************************/
    const SphericalEngine::coeff& Coefficients() const throw()
    { return _c[0]; }
  };

} // namespace GeographicLib

#endif  // GEOGRAPHICLIB_SPHERICALHARMONIC_HPP