Discrete evaluation of Stokes’s integral by means of Voronoi and Delaunay structures

Two alternative approaches are investigated to compute the discrete Stokes integral for gravimetric geoid determination so that geographical grid subdivision and gridding is not required. The techniques are based on Voronoi and Delaunay structures, in which the target area is partitioned into polygons and triangles, respectively, and the computation is carried out by point-wise integration. In the Voronoi scheme, polygonal areas just contain the observed gravity anomalies, instead of the interpolated ones; thus no gridding process or data interpolation is necessary, and only the original data are used. In the Delaunay scheme, gridding is also not required, but observed anomalies are interpolated into triangular compartments whose vertices hold the gravity stations. Geoidal undulations are thus computed at the barycenters (centroids) of the triangles. Both schemes were applied to the local gravimetric geoid determination in two distinct areas of Brazil (municipality of Rio de Janeiro, and Resende). The gravity observations are almost uniformly distributed spatially at both sites, and their topographies are very rugged. The Stokes component was also computed by means of classical numerical integration (space-domain), and compared with the Voronoi and Delaunay schemes to give root-mean-square (RMS) differences of 0.022 and 0.024 m, respectively, at the Rio de Janeiro site. In Resende, the comparisons yielded RMS differences of 0.040 and 0.053 m. The largest difference did not exceed 0.100 m for both methods and datasets. The one-dimensional (1-D) FFT (spectral domain) technique was also used on the Rio de Janeiro dataset, which gave RMS differences of 0.031 m for the classical method, 0.039 m for the Voronoi scheme, and 0.047 m for the Delaunay scheme. Relative comparisons with 465 GPS-leveling baselines in the Rio de Janeiro site gave RMS differences of 0.069, 0.061, 0.071, and 0.071 m, for the Voronoi, Delaunay, classical, and 1-D FFT methods, respectively. Since the Voronoi and Delaunay schemes avoid the gridding step, the pre-processing time and labor are reduced. As with other methods, the dependence upon data quality and distribution is the main drawback of both schemes. Finally, the Voronoi and Delaunay schemes proved to be computationally as efficient as the 1-D FFT method for only the geoid height computation.     

Updating and computing the geoid using two dimensional fast Hartley transform and fast T transform

In the analyses of 2D real arrays, fast Hartley (FHT), fast T (FTT) and real-valued fast Fourier transforms are generally preferred in lieu of a complex fast Fourier transform due to the advantages of the former with respect to disk storage and computation time. Although the FHT and the FTT in one dimension are identical, they are different in two or more dimensions. Therefore, first, definitions and some properties of both transforms and the related 2D FHT and FTT algorithms are stated. After reviewing the 2D FHT and FTT solutions of Stokes' formula in planar approximation, 2D FHT and FTT methods are developed for geoid updating to incorporate additional gravity anomalies. The methods are applied for a test area which includes a 64×64 grid of 3^′×3^′ point gravity anomalies and geoid heights calculated from point masses. The geoids computed by 2D FHT and FTT are found to be identical. However, the RMS value of the differences between the computed and test geoid is ±15 mm. The numerical simulations indicate that the new methods of geoid updating are practical and accurate with considerable savings on storage requirements. Received: 15 February 1996; Accepted: 22 January 1997     

Minimization and estimation of geoid undulation errors

The objective of this paper is to minimize the geoid undulation errors by focusing on the contribution of the global geopotential model and regional gravity anomalies, and to estimate the accuracy of the predicted gravimetric geoid.The geopotential model's contribution is improved by (a) tailoring it using the regional gravity anomalies and (b) introducing a weighting function to the geopotential coefficients. The tailoring and the weighting function reduced the difference (1) between the geopotential model and the GPS/levelling-derived geoid undulations in British Columbia by about 55% and more than 10%, respectively.Geoid undulations computed in an area of 40° by 120° by Stokes' integral with different kernel functions are analyzed. The use of the approximated kernels results in about 25 cm () and 190 cm (maximum) geoid errors. As compared with the geoid derived by GPS/levelling, the gravimetric geoid gives relative differences of about 0.3 to 1.4 ppm in flat areas, and 1 to 2.5 ppm in mountainous areas for distances of 30 to 200 km, while the absolute difference (1) is about 5 cm and 20 cm, respectively.A optimal Wiener filter is introduced for filtering of the gravity anomaly noise, and the performance is investigated by numerical examples. The internal accuracy of the gravimetric geoid is studied by propagating the errors of the gravity anomalies and the geopotential coefficients into the geoid undulations. Numerical computations indicate that the propagated geoid errors can reasonably reflect the differences between the gravimetric and GPS/levelling-derived geoid undulations in flat areas, such as Alberta, and is over optimistic in the Rocky Mountains of British Columbia.Paper presented at the IAG General Meeting, Beijing, China, August 8–13, 1993.     

Gravity field convolutions without windowing and edge effects

A new set of formulas has been developed for the computation of geoid undulations and terrain corrections by FFT when the input gravity anomalies and heights are mean gridded values. The effects of the analytical and the discrete spectra of kernel functions and that of zero-padding on the computation of geoid undulations and terrain corrections are studied in detail.Numerical examples show that the discrete spectrum is superior to the analytically-defined one. By using the discrete spectrum and 100% zero-padding, the RMS differences are 0.000 m for the FFT geoid undulations and 0.200 to 0.000 mGal for the FFT terrain corrections compared with results obtained by numerical integration.     

A numerical investigation on height anomaly prediction in mountainous areas

This paper provides numerical examples for the prediction of height anomalies by the solution of Molodensky's boundary value problem. Computations are done within two areas in the Canadian Rockies. The data used are on a grid with various grid spacings from 100 m to 5 arc-minutes. Numerical results indicate that the Bouguer or the topographicisostatic gravity anomalies should be used in gravity interpolation. It is feasible to predict height anomalies in mountainous areas with an accuracy of 10 cm (1) if sufficiently dense data grids are used. After removing the systematic bias, the differences between the geoid undulations converted from height anomalies and those derived from GPS/levelling on 50 benchmarks is 12 cm (1) when the grid spacing is 1km, and 50 cm (1) when the grid spacing is 5. It is not necessary, in most cases, to require a grid spacing finer than 1 km, because the height anomaly changes only by 3 cm (1) when the grid spacing is increased from 100 m to 1000 m. Numerical results also indicate that, only the first two terms of the Molodensky series have to be evaluated in all but the extreme cases, since the contributions of the higher order terms are negligible compared to the objective accuracy.     

Combined regional geoid determination for the ERS-1 radar altimeter calibration

Summary A local model of the geoid in NE Italy and its section along the Venice ground track of the ERS-1 satellite of the European Space Agency is presented. The observational data consist of geoid undulations determined with a network of 25 stations of known orthometric (by spirit leveling) and ellipsoidal (by GPS differential survey) and of 13 deflections of the vertical measured at sites of the network for which, besides the ellipsoidal (WGS84) coordinates, also astronomic coordinates were known. The network covers an area of 1×1 degrees and is tied to a vertical and horizontal datum: one vertex of the network is the tide gauge of Punta Salute, in Venice, providing a tie to a mean sea level; a second vertex is the site for mobile laser systems at Monte Venda, on the Euganei Hills, for which geocentric coordinates resulted from the analysis of several LAGEOS passes.The interpolation algorithm used to map sparse and heterogeneous data to a regular grid of geoid undulations is based on least squares collocation and the autocorrelation function of the geoid undulations is modeled by a third order Markov process on flat earth. The algorithm has been applied to the observed undulations and deflections of the vertical after subtraction of the corresponding predictions made on the basis of the OSU91A global geoid model of the Ohio State University, complete to degree and order 360. The locally improved geoid results by adding back, at the nodes of a regular grid, the predictions of the global field to the least squares interpolated values. Comparison of the model values with the raw data at the observing stations indicates that the mean discrepancy is virtually zero with a root mean square dispersion of 8 cm, assuming that the ellipsoidal heights and vertical deflections data are affected by a random error of 3 cm and 0.5 respectively. The corrections resulting from the local data and added to the background 360×360 global model are described by a smooth surface with excursions from the reference surface not larger than ±30 cm.     

高原山区局部大地水准面精化方法研究

局部大地水准面精化的实质是精确计算出大地水准面的起伏变化情况。一般情况下,需要密度足够的重力数据,依重力异常密集计算大地水准面差距或高程异常。但是在大陆西部高原山区重力点密度是不够的,无法达到大地水准面精化的目的。本文从理论上证实了用地形和岩石密度数据进行局部大地水准面精化的可行性。     

Inverse Vening Meinesz formula and deflection-geoid formula: applications to the predictions of gravity and geoid over the South China Sea

Using the spherical harmonic representations of the earth's disturbing potential and its functionals, we derive the inverse Vening Meinesz formula, which converts deflection of the vertical to gravity anomaly using the gradient of the H function. The deflection-geoid formula is also derived that converts deflection to geoidal undulation using the gradient of the C function. The two formulae are implemented by the 1D FFT and the 2D FFT methods. The innermost zone effect is derived. The inverse Vening Meinesz formula is employed to compute gravity anomalies and geoidal undulations over the South China Sea using deflections from Seasat, Geosat, ERS-1 and TOPEX//POSEIDON satellite altimetry. The 1D FFT yields the best result of 9.9-mgal rms difference with the shipborne gravity anomalies. Using the simulated deflections from EGM96, the deflection-geoid formula yields a 4-cm rms difference with the EGM96-generated geoid. The predicted gravity anomalies and geoidal undulations can be used to study the tectonic structure and the ocean circulations of the South China Sea. Received: 7 April 1997 / Accepted: 7 January 1998     

Effect of the SRTM global DEM on the determination of a high-resolution geoid model: a case study in Iran

Any errors in digital elevation models (DEMs) will introduce errors directly in gravity anomalies and geoid models when used in interpolating Bouguer gravity anomalies. Errors are also propagated into the geoid model by the topographic and downward continuation (DWC) corrections in the application of Stokes’s formula. The effects of these errors are assessed by the evaluation of the absolute accuracy of nine independent DEMs for the Iran region. It is shown that the improvement in using the high-resolution Shuttle Radar Topography Mission (SRTM) data versus previously available DEMs in gridding of gravity anomalies, terrain corrections and DWC effects for the geoid model are significant. Based on the Iranian GPS/levelling network data, we estimate the absolute vertical accuracy of the SRTM in Iran to be 6.5 m, which is much better than the estimated global accuracy of the SRTM (say 16 m). Hence, this DEM has a comparable accuracy to a current photogrammetric high-resolution DEM of Iran under development. We also found very large differences between the GLOBE and SRTM models on the range of −750 to 550 m. This difference causes an error in the range of −160 to 140 mGal in interpolating surface gravity anomalies and −60 to 60 mGal in simple Bouguer anomaly correction terms. In the view of geoid heights, we found large differences between the use of GLOBE and SRTM DEMs, in the range of −1.1 to 1 m for the study area. The terrain correction of the geoid model at selected GPS/levelling points only differs by 3 cm for these two DEMs.     

A new isostatic model of the lithosphere and gravity field

Based on the analysis of various factors controlling isostatic gravity anomalies and geoid undulations, it is concluded that it is essential to model the lithospheric density structure as accurately as possible. Otherwise, if computed in the classical way (i.e. based on the surface topography and the simple Airy compensation scheme), isostatic anomalies mostly reflect differences of the real lithosphere structure from the simplified compensation model, and not necessarily the deviations from isostatic equilibrium. Starting with global gravity, topography and crustal density models, isostatic gravity anomalies and geoid undulations have been determined. The initial crust and upper-mantle density structure has been corrected in a least squares adjustment using gravity. To model the long-wavelength (>2000 km) features in the gravity field, the isostatic condition (i.e. equal mass for all columns above the compensation level) is applied in the adjustment to uncover the signals from the deep-Earth interior, including dynamic deformations of the Earths surface. The isostatic gravity anomalies and geoid undulations, rather than the observed fields, then represent the signals from mantle convection and deep density inhomogeneities including remnants of subducted slabs. The long-wavelength non-isostatic (i.e. the dynamic) topography was estimated to range from –0.4 to 0.5 km. For shorter wavelengths (<2000 km), the isostatic condition is not applied in the adjustment in order to obtain the non-isostatic topography due to regional deviations from classical Airy isostasy. The maximum deviations from Airy isostasy (–1.5 to 1 km) occur at currently active plate boundaries. As another result, a new global model of the lithosphere density distribution is generated. The most pronounced negative density anomalies in the upper mantle are found near large plume provinces, such as Iceland and East Africa, and in the vicinity of the mid-ocean ridge axes. Positive density anomalies in the upper mantle under the continents are not correlated with the cold and thick lithosphere of cratons, indicating a compensation mechanism due to thermal and compositional density.     

A new,high-resolution geoid for Canada and part of the U.S. by the 1D-FFT method

A new, high-resolution and high-precision geoid has been computed for the whole of Canada and part of the U.S., ranging from 35°N to about 90°N in latitude and 210°E to 320°E in longitude. The OSU91A geopotential model complete to degree and order 360 was combined with a 5 × 5 mean gravity anomaly grid and 1km × 1km topographical information to generate the geoid file. The remove-restore technique was adopted for the computation of terrain effects by Helmert's condensation reduction. The contribution of the local gravity data to the geoid was computed strictly by the 1D-FFT technique, which allows for the evaluation of the discrete spherical Stokes integral without any approximation, parallel by parallel. The indirect effects of up to second order were considered. The internal precision of the geoid, i.e. the contribution of the gravity data and the model coefficients noise, was also evaluated through error propagation by FFT. In a relative sense, these errors seem to agree quite well with the external errors and show clearly the weak areas of the geoid which are mostly due to insufficient gravity data coverage. Comparison of the gravimetric geoid with the GPS/levelling-derived geoidal heights of eight local GPS networks with a total of about 900 stations shows that the absolute agreement with respect to the GPS/levelling datum is generally better than 10 cm RMS and the relative agreement ranges, in most cases, from 4 to 1 ppm over short distances of about 20 to 100km, 1 to 0.5 ppm over distances of about 100 to 200 km, and 0.5 to 0.1 ppm for baselines of 200 to over 1000 km. Other existing geoids, such as UNB90, GEOID90 and GSD91, were also included in the comparison, showing that the new geoid achieves the best agreement with the GPS/levelling data.Presented at theIAG General Meeting, Beijing, P.R. China, Aug. 6–13, 1993     

Comparisons of spectral techniques for geoid computations over large regions

The Stokes formula is efficiently evaluated by the one-and two- dimensional (1D, 2D) fast Fourier transform (FFT) technique in the plane and on the sphere in order to obtain precise geoid determinatiover a large area such as Europe. Using a high-pass filtered spherical harmonic reference model (OSU91A truncated to different degrees), gridded gravity anomalies and geoid heights were produced and the anomalies were used as input in the FFT software. Various tests were performed with respect to the different kernel functions used, to the spherical computations in bands, as well as to windowing, edge effects and extent of the area. It is thus demonstrated that, in geoid computations over large regions, the 1D spherical FFT and the 2D multiband spherical FFT in combination with discrete spectra for the kernel functions and 100% zero-padding give better results than those obtained by the other transform techniques. Additionally, numerical tests were carried out at the same test area using the planar fast Hartley transform (FHT) instead of the FFT and the results obtained by the two attractive alternatives were compared regarding the requirements in both computer time and computer memory needed in geoid height computations.A slightly modified version of the paper has been presented at the XX EGS General Assembly, Hamburg, 3–7 April, 1995     

Harmonic continuation and gridding effects on geoid height prediction

Least-squares collocation and Stokes integral formula, as implemented using the Fast Fourier Technique, handle the harmonic downward continuation problem quite differently. FFT furthermore requires gridded data, amplifying the difference of methods.We have in this paper studied numerically the effects of downward continuation and gridding in a mountainous area in central Norway. Topographically smoothed data were used in order to reduce these effects. Despite the smoothing, it was found that the vertical gravity gradient had values up to -11 mgal/km. The corresponding differences between geoid heights and the height anomalies at altitude reached 12 cm.The differences between geoid heights obtained using collocation or FFT with gravity data at terrain level or sea level showed differences between the values of up to 10 cm r.m.s. A part of this difference was a consequence of different data areas used in the FFT and collocation solution, though.Major discrepancies between the solutions were found in areas where the topographic smoothing could not be applied (deep fjords with no depth information in the used DTM) or where there seemed to be gross errors in the data.We conclude that proper handling of harmonic continuation is important, even when we as here have used a 1 km resolution DTM for the calculation of topographic effects. The effect of data gridding, required for the FFT method, seems not to be as serious as the need to limit the data distribution area, required when least squares collocation is used with randomly distributed data.     

计算Stokes公式的快速Hartley变换(FHT)技术

本文提出了利用快速Hartley变换(FHT)计算Stokes公式的方法,这一算法最适合于用来计算实序列的积分变换,而快速Fourier变换(FFT)较适合于用来计算复序列的积分变换。计算Stokes公式只涉及实序列问题,用FHT计算Stokes公式比用FFT算法更有效。本文详细地描述了用FHΥ计算Stokes公式的算法,进行了数值计算,与相应的FFT计算结果作了比较。结果表明,两种算法可以得到相同的精度,但是,FHT的计算速度比FFT的计算速度快一倍以上,且所需要的内存空间只是后者的一半。     

Gravity anomalies from satellite altimetry: comparison between computation via geoid heights and via deflections of the vertical

The accumulation of good quality satellite altimetry missions allows us to have a precise geoid with fair resolution and to compute free air gravity anomalies easily by fast Fourier transform (FFT) techniques.In this study we are comparing two methods to get gravity anomalies. The first one is to establish a geoid grid and transform it into anomalies using inverse Stokes formula in the spectral domain via FFT. The second one computes deflection of the vertical grids and transforms them into anomalies.The comparison is made using different data sets: Geosat, ERS-1 and Topex-Poseidon exact repeat misions (ERMs) north of 30°S and Geosat geodetic mission (GM) south of 30°S. The second method which transforms the geoid gradients converted into deflection of the vertical values is much better and the results have been favourably evaluated by comparison with marine gravity data.     

A new strategy for the atmospheric gravity effect in gravimetric geoid determination

Prior to Stokes integration, the gravitational effect of atmospheric masses must be removed from the gravity anomaly g. One theory for the atmospheric gravity effect on the geoid is the well-known International Association of Geodesy approach in connection with Stokes integral formula. Another strategy is the use of a spherical harmonic representation of the topography, i.e. the use of a global topography computed from a set of spherical harmonics. The latter strategy is improved to account for local information. A new formula is derived by combining the local contribution of the atmospheric effect computed from a detailed digital terrain model and the global contribution computed from a spherical harmonic model of the topography. The new formula is tested over Iran and the results are compared with corresponding results from the old formula which only uses the global information. The results show significant differences. The differences between the two formulas reach 17 cm in a test area in Iran.     

Confirming regional 1 cm differential geoid accuracy from airborne gravimetry: the Geoid Slope Validation Survey of 2011

A terrestrial survey, called the Geoid Slope Validation Survey of 2011 (GSVS11), encompassing leveling, GPS, astrogeodetic deflections of the vertical (DOV) and surface gravity was performed in the United States. The general purpose of that survey was to evaluate the current accuracy of gravimetric geoid models, and also to determine the impact of introducing new airborne gravity data from the ‘Gravity for the Redefinition of the American Vertical Datum’ (GRAV-D) project. More specifically, the GSVS11 survey was performed to determine whether or not the GRAV-D airborne gravimetry, flown at 11 km altitude, can reduce differential geoid error to below 1 cm in a low, flat gravimetrically uncomplicated region. GSVS11 comprises a 325 km traverse from Austin to Rockport in Southern Texas, and includes 218 GPS stations ( $\sigma _{\Delta h }= 0.4$ cm over any distance from 0.4 to 325 km) co-located with first-order spirit leveled orthometric heights ( $\sigma _{\Delta H }= 1.3$ cm end-to-end), including new surface gravimetry, and 216 astronomically determined vertical deflections $(\sigma _{\mathrm{DOV}}= 0.1^{\prime \prime })$ . The terrestrial survey data were compared in various ways to specific geoid models, including analysis of RMS residuals between all pairs of points on the line, direct comparison of DOVs to geoid slopes, and a harmonic analysis of the differences between the terrestrial data and various geoid models. These comparisons of the terrestrial survey data with specific geoid models showed conclusively that, in this type of region (low, flat) the geoid models computed using existing terrestrial gravity, combined with digital elevation models (DEMs) and GRACE and GOCE data, differential geoid accuracy of 1 to 3 cm (1 $\sigma )$ over distances from 0.4 to 325 km were currently being achieved. However, the addition of a contemporaneous airborne gravity data set, flown at 11 km altitude, brought the estimated differential geoid accuracy down to 1 cm over nearly all distances from 0.4 to 325 km.     

APPLICATION OF THE HARTLEY TRANSFORM IN THE COMPUTATION OF THE GEOID IN CHINA

This paper presents a method for the computation of the Stokes for-mula using the Fast Hartley Transform（FHT）techniques.The algorithm is mostsuitable for the computation of real sequence transform,while the Fast FourierTransform（FFT）techniques are more suitable for the computaton of complex se-quence transform.A method of spherical coordinate transformation is presented inthis paper.By this method the errors,which are due to the approximate term inthe convolution of Stokes formula,can be effectively eliminated.Some numericaltests are given.By a comparison with both FFT techniques and numerical integra-tion method,the results show that the resulting values of geoidal undulations byFHT techniques are almost the same as by FFT techniques,and the computation-al speed of FHT techniques is about two times faster than that of FFT techniques.     

FFT-Evaluation and applications of gravity-field convolution integrals with mean and point data

Due to the fact that the spectrum of a convolution is the product of the spectra of the two convolved functions, the convolution integrals of physical geodesy can be evaluated very efficiently by the use of the fast Fourier transform (FFT) provided that gravity and/or terrain data are available on a regular grid. All Fourier transform-based methods usually treat the gridded data as point values despite the fact that these discrete values may have been obtained by averaging and they represent mean values over the whole area of a grid element. In the frequency domain, this fact can be taken into account very easily, because the spectra of mean and point data are related via a two-dimensional (2D) sinc function. The paper shows explicitly this relationship using the convolution integrals for the evaluation of geoid undulations, deflections of the vertical, and gravity and gradiometry terrain effects. Numerical tests are presented, indicating that the differences in the two approaches may become significant when highly accurate results are wanted. The application of the2D sinc function in the evaluation, update, and inversion of other convolution integrals is briefly discussed as well.     

High-resolution regional geoid computation without applying Stokes’s formula: a case study of the Iranian geoid

A new theory for high-resolution regional geoid computation without applying Stokess formula is presented. Operationally, it uses various types of gravity functionals, namely data of type gravity potential (gravimetric leveling), vertical derivatives of the gravity potential (modulus of gravity intensity from gravimetric surveys), horizontal derivatives of the gravity potential (vertical deflections from astrogeodetic observations) or higher-order derivatives such as gravity gradients. Its algorithmic version can be described as follows: (1) Remove the effect of a very high degree/order potential reference field at the point of measurement (POM), in particular GPS positioned, either on the Earths surface or in its external space. (2) Remove the centrifugal potential and its higher-order derivatives at the POM. (3) Remove the gravitational field of topographic masses (terrain effect) in a zone of influence of radius r. A proper choice of such a radius of influence is 2r=4×10⁴ km/n, where n is the highest degree of the harmonic expansion. (cf. Nyquist frequency). This third remove step aims at generating a harmonic gravitational field outside a reference ellipsoid, which is an equipotential surface of a reference potential field. (4) The residual gravitational functionals are downward continued to the reference ellipsoid by means of the inverse solution of the ellipsoidal Dirichlet boundary-value problem based upon the ellipsoidal Abel–Poisson kernel. As a discretized integral equation of the first kind, downward continuation is Phillips–Tikhonov regularized by an optimal choice of the regularization factor. (5) Restore the effect of a very high degree/order potential reference field at the corresponding point to the POM on the reference ellipsoid. (6) Restore the centrifugal potential and its higher-order derivatives at the ellipsoidal corresponding point to the POM. (7) Restore the gravitational field of topographic masses ( terrain effect) at the ellipsoidal corresponding point to the POM. (8) Convert the gravitational potential on the reference ellipsoid to geoidal undulations by means of the ellipsoidal Bruns formula. A large-scale application of the new concept of geoid computation is made for the Iran geoid. According to the numerical investigations based on the applied methodology, a new geoid solution for Iran with an accuracy of a few centimeters is achieved.Acknowledgments. The project of high-resolution geoid computation of Iran has been support by National Cartographic Center (NCC) of Iran. The University of Tehran, via grant number 621/3/602, supported the computation of a global geoid solution for Iran. Their support is gratefully acknowledged. A. Ardalan would like to thank Mr. Y. Hatam, and Mr. K. Ghazavi from NCC and Mr. M. Sharifi, Mr. A. Safari, and Mr. M. Motagh from the University of Tehran for their support in data gathering and computations. The authors would like to thank the comments and corrections made by the four reviewers and the editor of the paper, Professor Will Featherstone. Their comments helped us to correct the mistakes and improve the paper.