Small changes suggestion in Documentation/ Docstring 1

[thom-le]

(see German below)

Hi,

while working through the geoid computation I compared the documentation of the geoidHeight kernel with the code, and the two do not quite match. The implementation is fine — it is only the docstring that is misleading. Three small points, all in source/classes/kernel/kernelGeoid.h, block docstringKernelGeoid:

1. The quadratic height term is missing. The docstring stops after the linear term:

\gamma  = \gamma_0 - 0.30877\cdot 10^{-5}/s^2(1-0.00142\sin^2(B))h

whereas Planets::normalGravity in source/base/planets.cpp includes a second-order term:

return gamma0 - 2*ga/a*(1+1/f+m+(-3/f+5*m/2)*sin2)*h + 3*ga/a/a*h*h;

3*ga/a/a evaluates to 7.212e-13, i.e. the classical 0.72e-12. At h = 3000 m the omitted term is about 0.65 mGal, so a reader who trusts the docstring may draw the wrong conclusion about the accuracy at high altitudes.

2. \gamma_0 is given as a series expansion, while the code uses the closed Somigliana formula. Since that formula introduces $a$, $b$, $\gamma_a$ and $\gamma_b$, the explanatory line below should name them as constants of GRS80.

3. B is called "the longitude" in the explanatory line; it is the ellipsoidal latitude, since normalGravity converts via ellipsoid(p, L, B, h) beforehand.

Suggested replacement:

 \begin{equation}\label{normalgravity}
-\gamma  = \gamma_0 - 0.30877\cdot 10^{-5}/s^2(1-0.00142\sin^2(B))h
+\gamma  = \gamma_0 - 0.30877\cdot 10^{-5}/s^2(1-0.00142\sin^2(B))h + 0.72\cdot 10^{-12}/(m\,s^2)h^2
 \end{equation}
 and
 \begin{equation}
-\gamma_0 = 9.780327\,m/s^2(1+0.0053024\sin^2(B)-0.0000058\sin^2(2B))
+\gamma_0 = \frac{a\gamma_a\cos^2(B)+b\gamma_b\sin^2(B)}{\sqrt{a^2\cos^2(B)+b^2\sin^2(B)}}
 \end{equation}
-where $h$ is the ellipsoidal height in meter and $B$ the longitude.
+where $h$ is the ellipsoidal height in meter, $B$ the ellipsoidal latitude,
+$a$, $b$ the semi-axes of the GRS80 ellipsoid and $\gamma_a=9.7803267715\,m/s^2$,
+$\gamma_b=9.8321863685\,m/s^2$ the normal gravity at the equator and at the
+pole, as given in GRS80 (Moritz 1980).

Thanks and best regards
Thomas

---

Hallo,

bei der Arbeit an der Geoidberechnung habe ich die Dokumentation des geoidHeight-Kernels mit dem Code verglichen, und beide passen nicht ganz zusammen. Die Implementierung ist in Ordnung — irreführend ist allein der Docstring. Drei kleine Punkte, alle in source/classes/kernel/kernelGeoid.h, Block docstringKernelGeoid:

1. Das quadratische Höhenglied fehlt. Der Docstring bricht nach dem linearen Term ab:

\gamma  = \gamma_0 - 0.30877\cdot 10^{-5}/s^2(1-0.00142\sin^2(B))h

Planets::normalGravity in source/base/planets.cpp enthält dagegen einen Term zweiter Ordnung:

return gamma0 - 2*ga/a*(1+1/f+m+(-3/f+5*m/2)*sin2)*h + 3*ga/a/a*h*h;

3*ga/a/a ergibt 7.212e-13, also das klassische 0.72e-12. Bei h = 3000 m macht das weggelassene Glied rund 0,65 mGal aus — wer sich auf den Docstring verlässt, schätzt die Genauigkeit im Hochgebirge falsch ein.

2. \gamma_0 steht als Reihenentwicklung, im Code ist es die geschlossene Formel von Somigliana. Da darin $a$, $b$, $\gamma_a$ und $\gamma_b$ auftreten, sollte die Erklärzeile darunter diese Größen als Konstanten des GRS80 benennen.

3. B wird als „the longitude" bezeichnet; es ist die ellipsoidische Breite, da normalGravity zuvor über ellipsoid(p, L, B, h) umrechnet.

Der Änderungsvorschlag steht als Diff im englischen Teil.

Danke und LG
Thomas

[tmayerguerr]

Hi Thomas,

Thank you for your comment. I have updated the documentation in line with your suggestion.

Best regards
Torsten

[thom-le]

Hi Torsten,

thank you very much. Just a small thing, I know, I just happened to come across it.

Best regards
Thomas