Assessment of Several Interpolation Methods for Precise GPS Orbit

Hamad Yousif; Ahmed El-Rabbany

doi:10.1017/S0373463307004250

Assessment of Several Interpolation Methods for Precise GPS Orbit

Published online by Cambridge University Press: 09 August 2007

Hamad Yousif and

Ahmed El-Rabbany

Show author details

Hamad Yousif: Affiliation:
(Ryerson University, Toronto)
Ahmed El-Rabbany*: Affiliation:
(Ryerson University, Toronto)
*: (Email: rabbany@ryerson.ca)

Article contents

Abstract
INTRODUCTION
PROPOSED INTERPOLATION METHODS
RESULTS AND DISCUSSION
INTERPOLATION USING THE BROADCAST EPHEMERIS
CONCLUSION
References

Rights & Permissions

Abstract

GPS applications such as Precise Point Positioning (PPP) require the availability of precise ephemeris at high rate. To support these applications, several institutions such as the International GNSS Service (IGS) have developed precise orbital service. Unfortunately, however, the data rate of such precise orbits is usually limited to 15 minutes. To overcome this problem, a number of orbital interpolation methods are proposed. This paper examines the performance of four interpolation methods for IGS precise GPS orbits, namely Lagrange, Newton Divided Difference, Cubic Spline and Trigonometric interpolation. In addition, the paper discusses a new approach, which utilizes the residuals between the broadcast and precise ephemeris to generate a high density precise ephemeris. It is shown that the new approach produces better results than previously reported orbital interpolation accuracy.

Keywords

Precise Point Positioning GPS Interpolation

Information

Type: Research Article
Information: The Journal of Navigation , Volume 60 , Issue 3 , September 2007 , pp. 443 - 455

DOI: https://doi.org/10.1017/S0373463307004250 [Opens in a new window]
Copyright: Copyright © The Royal Institute of Navigation 2007

1. INTRODUCTION

Some applications, such as studies of the crustal dynamics of the Earth, require more precise ephemeris data than the broadcast ephemeris. The Precise Point Positioning (PPP) technique has recently evolved, which also requires the availability of precise ephemeris (as well as satellite clock corrections). To support these applications, several institutions (e.g., the International GNSS Service (IGS), formerly International GPS Service) have developed post-mission precise orbital services. Precise ephemeris data is based on GPS data collected at a global GPS network coordinated by the IGS. At present, IGS precise ephemeris data is available to users with some delay, which varies from 3 hours for the observed half of the IGS ultra-rapid orbit, to 17 hours for the rapid orbit, to about 13 days for the most precise IGS final orbit (http://igscb.jpl.nasa.gov/igscb/resource/igssheet.pdf). The three types of precise orbits are referred to the ITRF reference system and have accuracies better than 5 cm. However, they differ in the accuracy of the satellite clock corrections. Users requiring real-time precise orbit data can use the predicted half of the IGS ultra-rapid orbit data, which are accurate to about 10 cm.

IGS precise orbits are available at a typical rate of 15 min. Unfortunately, many GPS applications require precise orbits at higher rates. Several methods of interpolation have been widely used to generate GPS positions and velocities at intermediate points. Previous studies based on the Lagrange method were used by the US National Geodetic Survey (NGS) (Remondi, Reference Remondi1991; Schueler, Reference Schueler1998). Polynomial and trigonometric interpolation were conducted by Schenewerk (Reference Schenewerk2003). Lagrange, Chebyshev and trigonometric were used by Feng and Zheng (Reference Feng and Zheng2005). In this paper, the 24-hour GPS precise orbit is used to generate a higher density orbit using four interpolation methods, namely Lagrange, Newton Divided Difference, Cubic Spline and Trigonometric Methods. Tests to determine the optimal number of terms that can achieve the best accuracy were conducted. In view of these tests, algorithms are designed and implemented. This paper also devises a new approach to obtain a high density precise ephemeris by means of the broadcast ephemeris. Unlike previous work, which interpolates the actual GPS orbit, the new approach interpolates the residuals between the broadcast and precise orbits using a Lagrange interpolation method. It is shown that the new approach produces better results within the two-hour broadcast ephemeris record than previously reported orbital interpolation accuracy.

2. PROPOSED INTERPOLATION METHODS

2.1. Lagrange Interpolation Method

The Lagrange formula is an algebraic expression that can be used to fit a particular data set to a certain polynomial (whose degree is equal to the number of data points) such that it returns the exact value of the function at each data point. Let f ₀, f ₁, f ₂, …, f _n be the values of a given data at times t ₀, t ₁, t ₂, …, t _n, respectively. The approximated value of f, denoted by p(t), at any time t is given by (Spiegel, Reference Spiegel1999):

(1)

$p\lpar t\rpar \equals a_{\setnum{0}} f_{\setnum{0}} \plus a_{\setnum{1}} f_{\setnum{1}\comma } \plus a_{\setnum{2}} f_{\setnum{2}\comma } \plus \cdots \plus a_{n} f_{n} \equals \sum\limits_{i \equals \setnum{0}}^{n} {a_{i} } f_{i}$

where:

(2)

$a_{i} \equals {{\lpar t \minus t_{\setnum{0}} \rpar \lpar t \minus t_{\setnum{1}} \rpar \cdots \lpar t \minus t_{i \minus \setnum{1}} \rpar \lpar t \minus t_{i \plus \setnum{1}} \rpar \cdots \lpar t \minus t_{n} \rpar } \over {\lpar t_{i} \minus t_{\setnum{0}} \rpar \lpar t_{i} \minus t_{\setnum{1}} \rpar \cdots \lpar t_{i} \minus t_{i \minus \setnum{1}} \rpar \lpar t_{i} \minus t_{i \plus \setnum{1}} \rpar \cdots \lpar t_{i} \minus t_{n} \rpar}}$

Since a _i coefficient is a function of t, it can also be referred to as L _i(t) which is the Lagrange operator. Now, substituting t by t ₀, t ₁, t ₂, …. t _n in equation (2) we get:

(3)

$a_{i} \equals L_{i} \lpar t\rpar \equals \left\{ {\matrix{ 1 \tab {for\,t \equals t_{i} } \cr 0 \tab {otherwise} \cr} } \right\}$

Going back to equation (1) and substituting again t by t ₀, t ₁, t ₂, …. t _n we get:

(4)

$p\lpar t_{\setnum{0}} \rpar \equals f_{\setnum{0}} \comma \ p\lpar t_{\setnum{1}} \rpar \equals f_{\setnum{1}} \comma \ p\lpar t_{\setnum{2}} \rpar \equals f_{\setnum{2}} \comma \cdots \comma \ p\lpar t_{n} \rpar \equals f_{n}$

2.2. Newton Divided Difference

The Newton Divided Difference formula is another polynomial expression of a given data set. Similar to Lagrange, Newton Divided Difference returns the exact value of the function at every data point. This formula is termed divided difference because its coefficients involve division of differences. The Newton Divided Difference formula is given, for a given data set (t _i,f _i), by (Spiegel, Reference Spiegel1999):

(5)

$\eqalign{ p\lpar t\rpar \equals \tab a_{\setnum{0}} \plus a_{\setnum{1}} \lpar t\minus t_{\setnum{0}} \rpar \plus a_{\setnum{2}} \lpar t\minus t_{\setnum{0}} \rpar \lpar t\minus t_{\setnum{1}} \rpar \plus a_{\setnum{3}} \lpar t\minus t_{\setnum{0}} \rpar \lpar t\minus t_{\setnum{1}} \rpar \lpar t\minus t_{\setnum{2}} \rpar \cr \tab \plus a_{n} \lpar t\minus t_{\setnum{0}} \ \rpar \lpar t\minus t_{\setnum{1}} \rpar \cdots \lpar t\minus t_{n \minus \setnum{1}} \rpar}$

where: a ₀, a ₁, a ₂, …, a _n are coefficients; and p(t) is the approximated value at any time t. The coefficients a ₀, a ₁, a ₂, …, a _n can be determined recursively by substituting t ₀, t ₁, t ₂, …, t _n as follows:

(6)

$a_{\setnum{0}} \equals p\lpar t_{\setnum{0}} \rpar$

(7)

$a_{\setnum{1}} \equals {{p\lpar t_{\setnum{1}} \rpar \minus a_{\setnum{0}} } \over {\lpar t_{\setnum{1}} \minus t_{\setnum{0}} \rpar}}$

(8)

$a_{\setnum{2}} \equals {{p\lpar t_{\setnum{2}} \rpar \minus a_{\setnum{0}} \minus a_{\setnum{1}} \lpar t_{\setnum{2}} \minus t_{\setnum{0}} \rpar } \over {\lpar t_{\setnum{2}} \minus t_{\setnum{0}} \rpar \lpar t_{\setnum{2}} \minus t_{\setnum{1}} \rpar}}$

⋮

(9)

$a_{n} \equals {{p\lpar t_{n} \rpar \minus a_{\setnum{0}} \minus a_{\setnum{1}} \lpar t_{n} \minus t_{\setnum{0}} \rpar \minus a_{\setnum{2}} \lpar t_{n} \minus t_{\setnum{0}} \rpar \lpar t_{n} \minus t_{\setnum{1}} \rpar \minus \cdots } \over {\lpar t_{n} \minus t_{\setnum{0}} \rpar \lpar t_{n} \minus t_{\setnum{1}} \rpar \lpar t_{n} \minus t_{\setnum{2}} \rpar \cdots \lpar t_{n} \minus t_{n \minus \setnum{1}} \rpar }}$

2.3. Cubic Spline Interpolation

In this method the data set is represented with a 3rd degree polynomial piecewisely (i.e. the data is divided into sections and each section is represented by a cubic polynomial). For example if we have a data set (t _i,f _i), then the Cubic Spline can be defined as follows (Press et al., Reference Press, Teukolosky, Vetterling and Flannery2002):

(10)

$S\lpar t\rpar \equals \left\{ {\matrix{ {s_{\setnum{0}} \lpar t\rpar } \tab {t_{\setnum{0}} \lt t \lt t_{\setnum{1}} } \hfill\cr {s_{\setnum{1}} \lpar t\rpar } \tab {t_{\setnum{1}} \lt t \lt t_{\setnum{2}} } \hfill\cr {s_{\setnum{2}} \lpar t\rpar } \tab {t_{\setnum{2}} \lt t \lt t_{\setnum{3}} } \hfill\cr {} \tab \vdots \cr {s_{i} \lpar t\rpar } \tab {t_{i} \lt t \lt t_{i \plus \setnum{1}} }\hfill \cr {} \tab \vdots \cr {s_{n \minus \setnum{1}} \lpar t\rpar } \tab {t_{n \minus \setnum{1}} \lt t \lt t_{n} } \cr} } \right\}$

where:

(11)

$s_{i} \lpar t\rpar \equals a_{i} t^{\setnum{3}} \plus b_{i} t^{\setnum{2}} \plus c_{i} t \plus d_{i}$

Each section of data has different coefficients a _i, b _i, c _i and d _i. These four coefficients can be determined as follows:

Substitute t=t _i and t=t _i+1

(12)

$s_{i} \lpar t_{i} \rpar \equals a_{i} {t_{i}} ^{\setnum{3}} \plus b_{i} {t_{i}} ^{\setnum{2}} \plus c_{i} t_{i} \plus d_{i}$

(13)

$s_{i} \lpar t_{i \plus \setnum{1}} \rpar \equals a_{i} t_{i \plus \setnum{1}}^{\setnum{3}} \plus b_{i} t_{i \plus \setnum{1}}^{\setnum{2}} \plus c_{i} t_{i \plus \setnum{1}} \plus d_{i}$

Differentiating s _i(t) and substituting t=t _i and t=t _i+1

(14)

$s\prime_{i} \lpar t_{i} \rpar \equals 3a_{i} {t_{i}} ^{\setnum{2}} \plus 2b_{i} t_{i} \plus c_{i}$

(15)

$s\prime_{i} \lpar t_{i \plus \setnum{1}} \rpar \equals 3a_{i} t_{i \plus \setnum{1}}^{\setnum{2}} \plus 2b_{i} t_{i \plus \setnum{1}} \plus c_{i}$

Now using the numerical differentiation formula:

(16)

$s\prime_{i} \lpar t_{i} \rpar \equals {{f_{i \plus \setnum{1}} \minus f_{i} } \over {t_{i \plus \setnum{1}} \minus t_{i}}}$

(17)

${\rm and}\hfill \quad s\prime_{i} \lpar t_{i \plus \setnum{1}} \rpar \equals {{f_{i \plus \setnum{2}} \minus f_{i \plus \setnum{1}} } \over {t_{i \plus \setnum{2}} \minus t_{i \plus \setnum{1}} }}\hfill$

Equations (12) to (17) can be solved simultaneously to determine the coefficients a _i, b _i, c _i and d _i of each section.

2.4. Trigonometric Interpolation

The trigonometric formula is a harmonic representation of a given data set. This representation is mainly convenient for the data that has a periodic nature, which is the case for precise in the inertial frame. The trigonometric formula for a given data set (t _i,f _i), is given by (Spiegel, Reference Spiegel1999):

(18)

$p\lpar t\rpar \equals {{a_{\setnum{0}} } \over 2} \plus \sum\limits_{n \equals \setnum{1}}^{M} {\left[ {a_{n} \cos \left( {{{2\pi n} \over T}t} \right) \plus b_{n} \sin \left( {{{2\pi n} \over T}t} \right)} \right]}$

where:

p(t): is the approximated value at time t
a _n, b _n: are coefficients
T: is the period of f
N: is the number of collocation points
M: is the number of truncated terms (ideally M=∞)

To determine a _n multiply both sides by $\cos \left( {{\textstyle{{2\pi m} \over T}}} \right)t$ and take the integration over one period:

(19)

$\eqalign{ \int_{\setnum{0}}^{T}\!\! {p\lpar t\rpar \cos } \left( {{{2\pi m} \over T}} \right)t\;dt \equals \tab \int_{\setnum{0}}^{T} {{{a_{\setnum{0}} } \over 2}\cos \left( {{{2\pi m} \over T}} \right)} tdt \cr \tab \plus \sum\limits_{n \equals \setnum{1}}^{M} {\int_{\setnum{0}}^{T}\!\!\! {a_{n} \cos \left( {{{2\pi n} \over T}t} \right)} \cos \left( {{{2\pi m} \over T}} \right)tdt} \cr \tab \plus \sum\limits_{n \equals \setnum{1}}^{M} {\int_{\setnum{0}}^{T} \!\!\!{a_{n} \sin \left( {{{2\pi n} \over T}t} \right)\cos \left( {{{2\pi m} \over T}} \right)} tdt} \cr}$

The first integration on the right hand side is equal to 0. Moreover according to the orthogonality theorem:

(20)

$\int_{\setnum{0}}^{T} {\sin \left( {{{2\pi n} \over T}t} \right)\cos \left( {{{2\pi m} \over T}} \right)tdt \equals 0}$

(21)

$\int_{\setnum{0}}^{T} {\cos \left( {{{2\pi n} \over T}t} \right)\cos \left( {{{2\pi m} \over T}} \right)} tdt \equals \left\{ {\matrix{ {0 \quad m \ne n} \cr {{T \over 2} \quad m \equals n} \cr} } \right\}$

making use of equations (19) and (20) we get:

(22)

$a_{n} \equals {2 \over T}\int_{\setnum{0}}^{T} {p\lpar t\rpar \cos \left( {{{2\pi n} \over T}} \right)} dt$

similarly:

(23)

$b_{n} \equals {2 \over T}\int_{\setnum{0}}^{T} {p\lpar t\rpar \sin \left( {{{2\pi n} \over T}} \right)} dt$

Equations (22) and (23) can be reformulated into the following numerical form:

(24)

$a_{n} \equals {2 \over N}\sum\limits_{k \equals \setnum{1}}^{N} {f_{k} \cos \left( {{{2\pi k} \over T}t_{k} } \right)}$

(25)

$b_{n} \equals {2 \over N}\sum\limits_{k \equals \setnum{1}}^{N} {f_{k} \sin \left( {{{2\pi k} \over T}t_{k} } \right)}$

3. RESULTS AND DISCUSSION

To assess the performance of the various interpolation methods, a test data set was acquired from the US National Geodetic Survey (NGS), which is available at http://www.ngs.noaa.gov/gps-toolbox/sp3intrp. The data set consists of the following files:

a. IN_15MIN.200 (A 24-hour precise orbit file spaced at 15 min in the inertial reference frame);
b. ECF_15MIN.200 (same as a but the orbital coordinate refer to the ECEF reference frame);
c. IN_5MIN.200 (same as a but the data interval is 5 minutes);
d. ECF_5MIN.200 (same as b but the data interval is 5 minutes).

The four files (a to d) were generated from the IGS sp3 rapid ephemeris (igr11472.sp3) for January 1, 2002 [http://igscb.jpl.nasa.gov]. The first two files (a and b) were used as input for each interpolation function, while the third and fourth files (c and d) were used as sources to which the interpolation results were compared.

Upon examination of the interpolation results, it was found that the following properties are shared by all interpolation methods:

■ Having too few collocation points produces an unreliable interpolation output. On the other hand, having relatively many points is expected to improve the result but to some limit.
■ The accuracy degrades noticeably near the end points and tends to improve as the interpolator moves towards the centre.

3.1. Lagrange Interpolation Results

As mentioned previously the interpolation accuracy enhances as the number of points increases. The accuracy also assumes its maximum at a specific number of terms after which it continues to degrade and eventually become unstable beyond a certain limit. Trying to investigate the optimum number of terms that gives the best accuracy, we started interpolation by taking the whole number of points for the 24-hour trajectory of 96 points spaced at 15 minutes intervals. We found that 96 points are too many to give an acceptable result as shown in Figure 1.

Figure 1.

Resulting error when taking the whole 24-hour data (i.e., 96 points) to build 96 Lagrange interpolating terms.

Taking smaller portions of data with different number of collocation points (8, 9, 10 and 11 points) we found that the best accuracy was achieved when n=9, as shown in Figures 2 and 3 (notice the difference in vertical scales of both figures). We can also notice that the accuracy degrades slightly near the boundaries for all figures. To overcome this problem, we segmented the 24-hour orbit into 23 overlapping portions; each has 9 terms as shown in Figure 4.

Figure 2.

Resulting error when taking 11 Lagrange terms.

Figure 3.

Resulting error when taking 9 Lagrange terms.

Figure 4.

Segmentation scheme – the 24-hour data is segmented into 23 overlapping segments of 9 terms each.

The segmentation scheme is designed in such a manner that each segment starts at the midpoint of the previous segment. The last segment however starts at point 88 (instead of point 89) to end at point 96 in order to complete the desired 9 terms. Otherwise we would have been left with only 8 terms for the last segment which would lead to a relatively lower accuracy. The idea behind having overlapped segments is to avoid producing spikes near the endpoints, as shown in Figures 2 and 3. Since the accuracy of interpolation assumes it is best at the middle, values near the end points are replaced by those of the mid adjacent segments. A comparison between partitioned and overlapped segments is illustrated in Table 1 and Figures 5 and 6. From these table and figures the accuracy enhancement in the overlapping case is easily observed.

Figure 5.

Interpolation error when the GPS orbit is partitioned into smaller segments without overlapping.

Figure 6.

Interpolation error when the GPS orbit is segmented into overlapping segments.

Table 1.

Comparison between interpolation errors generated by partitioned vs. overlapping segmentation.

Applying Lagrange interpolation using the previous scheme for both inertial and ECEF coordinates, Figures 7 and 8 show the interpolation error obtained for PRN01 (notice the difference in vertical scales of both figures). The interpolation accuracy of inertial orbit is better than that of ECEF, the statistics in Table 2 also depicts this fact numerically.

Figure 7.

Lagrange interpolation error for inertial GPS orbit.

Figure 8.

Lagrange interpolation error for ECEF GPS orbit.

Table 2.

Comparison between inertial vs. ECEF Lagrange interpolation error.

3.2. Newton Divided Difference Interpolation Results

Using the previously discussed algorithm, which can be applied to any other interpolation function, we found that the Newton interpolation result is completely identical to that of Lagrange for both inertial and ECEF orbits as shown in Figures 9 and 10.

Figure 9.

Newton Divided Difference interpolation error for inertial GPS orbit.

Figure 10.

Newton Divided Difference interpolation error for ECEF GPS orbit.

3.3. Cubic Spline Interpolation Results

The Cubic Spline interpolation resulted in a very low accuracy as seen from Figure 11. Perhaps, the reason is due to the fact that the cubic spline is a polynomial of order 3 which is not suitable to represent the periodic nature of the GPS orbit.

Figure 11.

Cubic Spline interpolation error for GPS inertial orbit.

3.4. Trigonometric Interpolation Results

The periodic nature of the GPS orbit in the inertial reference frame makes it reasonable to employ the periodic trigonometric interpolation. The algorithm of applying trigonometric interpolation method for GPS orbit was designed by Mark Schenewerk (Reference Schenewerk2003). Using different values of M in Equation 18 to determine the number of terms that can give the best accuracy, it is found that M=9 terms is the one that satisfies this objective. The statistics shown in Table 3 justifies the choice of M=9.

Table 3.

Statistics of trigonometric interpolation using different number of terms.

The results of interpolating inertial and ECEF orbits using 9 terms are shown in Figures 12 and 13, respectively. Again the transformation error adds up to the interpolation error in case of ECEF orbit.

Figure 12.

Trigonometric interpolation error for GPS inertial orbit.

Figure 13.

Trigonometric interpolation error for GPS ECEF orbit.

4. INTERPOLATION USING THE BROADCAST EPHEMERIS

Previous discussion focused on the interpolation of the actual precise orbit using four different interpolation methods. However, as a result of the high GPS altitude, interpolating the actual orbit is expected to produce a relatively large interpolation error. To overcome this problem, the residuals between the broadcast and precise orbits are formed and interpolated.

The broadcast navigation file, which corresponds to the test data set for January 1, 2002, was downloaded from the Crustal Dynamics Data Information System server (ftp://cddis.gsfc.nasa.gov/pub/gps/gpsdata/brdc/). The orbital parameters of PRN01 were extracted and used to calculate the broadcast ephemeris at 15 minutes intervals. The residuals between the broadcast ephemeris and the precise rapid ephemeris of January 1, 2002 were computed, and then interpolated at 5 minutes intervals using Lagrange interpolation method. The orbital parameters were used again to calculate the broadcast ephemeris but this time at 5 minutes intervals. The broadcast ephemeris at 5 minutes intervals were then added to their respective interpolated residuals to give the ECEF coordinates at 5 minutes intervals. The output was then compared to the reference ECEF coordinates and the errors for the first 2 hours were plotted in Figure 14. As can be seen from Figure 14, the maximum interpolation error is about −1·6 mm, which shows that this approach produces better results than the previous approaches. With the improvement in the GPS broadcast ephemeris as part of the modernization programme, it is expected that the use of more recent navigation files would further reduce the interpolation error. It should be mentioned, however, that as the ephemeris record is valid for two hours only, a larger interpolation error would be expected whenever a new ephemeris record is used.

Figure 14.

Interpolation error using the broadcast ephemeris.

5. CONCLUSION

This paper examined the performance of four different interpolation methods, namely Lagrange, Newton Divided Difference, Cubic Spline and Trigonometric. Apart from Cubic Spline, the other three methods produce a relatively good accuracy. The structure of the Cubic Spline, which is essentially a third degree polynomial, has made it inappropriate to represent the GPS orbit. Lagrange and Newton Divided Difference demonstrate completely identical results in terms of interpolation error. Due to the periodic nature of the GPS orbit, the trigonometric has yielded the best accuracy of all four interpolation methods. One more reason that enhances the accuracy of the trigonometric method is that the 24-hour orbit was well centred among sufficient data by taking several hours from the previous and subsequent days. The interpolation via the new approach, which uses the residuals between the broadcast and precise ephemeris, produced the best results within the two-hour ephemeris record. The maximum interpolation error with the new approach was about 1·6 mm.

References

REFERENCES:

Feng, Y. and Zheng, Y. (2005). “Efficient Interpolation to GPS Orbits for Precise Wide Area Applications.” GPS Solutions, 9, 4, 273–282.CrossRef Google Scholar

Spiegel, M. R. (1999). Mathematical Handbook of Formulas and Tables. McGraw Hill.Google Scholar

Remondi, B. (1991). “NGS Second Generation ASCII and Binary Orbit Formats and Associated Interpolation Studies.” Poster presentation at the 20th General Assembly of the IUGG, Vienna, 11–24 August.Google Scholar

Schenewerk, M. (2003). “A Brief Review Of Basic GPS Orbit Interpolation Strategies.” GPS Solutions, 6, 4, 265–267.CrossRef Google Scholar

Schueler, T (1998). “On the interpolation of SP3 orbit files.” IfEN Technical IFEN-TropAC-TN – 002 – 01. Institute of Geodesy and Navigation, University of Federal Armed Forces, Munich.Google Scholar

Press, W. H., Teukolosky, S. A., Vetterling, W. T., Flannery, B. P. (2002). Numerical Recipes in C++: The Art of Scientific Computing. Cambridge University Press.Google Scholar