# A Fast Method for Calculation of Marine Gravity Anomaly
## 2 Method
### The Inverse Stokes Formula  
According to [PERSON] et al. [11], the inverse Stokes formula, which uses geoid height undulation to calculate marine gravity anomaly, has the form of,  
\[\Delta g_{P}(\varphi_{p},\lambda_{p})=-\frac{\overline{\gamma}}{R}N_{p}-\frac {\overline{\gamma}}{16\pi R}\iint\limits_{\overline{\sigma}}\frac{N_{q}-N_{p} }{\sin^{3}\frac{\varphi}{2}} \tag{1}\]  
where \(\varphi_{p}\) and \(\lambda_{p}\) are latitude and longitude coordinates of the point \(P\) (the central point of any grid in the study area), respectively; \(\Delta g_{P}\) is the gravity anomaly at point \(P\); \(\overline{\gamma}\) is the average gravity of the Earth; \(R\) is the average radius of the Earth; \(N_{p}\) is the geoid height undulation at the point \(P\); and \(N_{q}\) is the geoid height undulation at a point \(Q\) (the central point of any grid in the integral area) in an integral range. \(\sin^{3}(\psi/2)\) can be computed by:  
\[\sin^{3}\frac{\psi}{2}=\left(\sin^{2}\frac{\varphi_{p}-\varphi_{q}}{2}+\sin^{2} \frac{\lambda_{p}-\lambda_{q}}{2}\cos\varphi_{p}\cos\varphi_{q}\right)^{\frac{ 3}{2}} \tag{2}\]  
where \(\varphi_{q}\) and \(\lambda_{q}\) are the latitude and longitude coordinates of point \(Q\), respectively. When using Equation (1) to compute the marine gravity anomaly, the integral calculation needs to be carried out in the global region. For ease of calculation, Equation (1) can be converted into a cumulative summation form as:  
\[\Delta g_{P}(\varphi_{p},\lambda_{p})=-\frac{\overline{\gamma}}{R}N_{p}-\frac {\Delta\varphi\Delta\lambda}{16\pi R}\overline{\gamma}\sum_{\varphi_{q}= \varphi_{1}}^{\varphi_{n}}\sum_{\lambda_{q}=\lambda_{1}}^{\lambda_{n}}\frac{ \left(N_{q}-N_{p}\right)\cos\varphi_{q}}{\sin^{3}\frac{\psi}{2}} \tag{3}\]  
where \(\Delta\varphi\) and \(\Delta\lambda\) are grid intervals along the latitude and longitude directions, respectively; and \(\varphi_{n}\) and \(\lambda_{n}\) are the numbers of the grid along the latitude and longitude directions, respectively.  
When point \(P\) and point \(Q\) coincide, they have the same latitude and longitude coordinates. The term \(\sin^{3}(\psi/2)\) in Equation (3) is equal to 0, which leads to a singularity problem. However, the term \(N_{q}-N_{p}\) in Equation (3) is also equal to 0, and it can be considered that the grid where point \(Q\) is located does not contribute to Equation (3). Therefore, when point \(P\) and point \(Q\) coincide, the corresponding component term \((N_{q}-N_{p})\cos\varphi_{q}/\sin^{3}(\psi/2)\) is set to 0, which can solve the singularity problem and does not affect the authenticity of the calculation of Equation (3). Judgment statements are generally used to determine whether point \(P\) and point \(Q\) coincide. However, the use of judgment statements in a program can increase the computing time. Moreover, the global geoid height undulation data are divided into lots of sufficiently small grids for obtaining a precise marine gravity anomaly. Therefore, the computing time of Equation (3) increases linearly with the number of the grids. To show the characteristics of the gravity anomaly of the study area in detail, the study area is divided into lots of sufficiently small grids to calculate the marine gravity anomalies, which leads to a further increase in the computing time.