Blažej Bucha

Lecturer at the Department of Theoretical Geodesy and Geoinformatics, Slovak University of Technology in Bratislava.

This page offers outputs from some of the research papers that I co-authored, all of them being related to gravity field modelling. The computer codes are written in C, Fortran and/or MATLAB. Anyone is encouraged to tailor the routines as freely as needed.

In case you have any comments, suggestions, questions or criticism, please feel free to contact me at If you too prefer secure communication, you may want to sign your email and/or encrypt it with my public_key using GNU Privacy Guard. The key fingerprint is EF9A 8C56 E625 8A23 A844 EFB3 24DD 5C83 468E 5E30.




Name: GrafLab
Description: Package for high-degree spherical harmonic synthesis
Programming Language: MATLAB
Version: 2.1.4
Updated: Mar 24, 2020

GrafLab (GRAvity Field LABoratory) is a MATLAB-based routine to compute functionals of the geopotential up to high degrees (tens of thousands and beyond).


A user can interact with GrafLab either through a GUI (the figure below) or from the command line.

GrafLab was published in

Bucha, B., Janák, J., 2013. A MATLAB-based graphical user interface program for computing functionals of the geopotential up to ultra-high degrees and orders. Computers & Geosciences 56, 186-196,

Please use this reference in your works when using GrafLab.

Some sample outputs from GrafLab

Computation time

Below are shown computation times required by GrafLab to synthesize the Earth's topography as represented by Earth2014.TBI2014.degree10800 up to harmonic degrees nmax = 90, 180, 360, 720, 1080, 2160, 4320, 6480, 8640 and 10800. The position of evaluation points is defined by the Gauss—Legendre quadrature. The dimensions of the grids are (nmax+1)*(2*nmax+2) in terms of the latitude and the longitude, respectively, and Legendre's functions were evaluated after Fukushima (2012). The tests were conducted on a PC with Intel® CoreTM i7-6800K CPU.

Go to top


Name: isGrafLab
Description: Package for efficient high-degree spherical harmonic synthesis on irregular surfaces
Programming Language: MATLAB
Version: 2.1.3
Updated: Nov 8, 2019

isGrafLab (Irregular Surface GRAvity Field LABoratory) is a modified version of GrafLab that allows accurate and fast computation of functionals of the geopotential at dense grids residing on irregular surfaces such as the Earth's surface. It employs the efficient lumped coefficients approach for surface spherical harmonic synthesis of the functionals and their radial derivatives at regular surfaces (sphere or ellipsoid of revolution), and the Taylor series expansions to continue the functionals to the irregular surface.

Further details can be found in the paper:

Bucha, B., Janák, J., 2014. A MATLAB-based graphical user interface program for computing functionals of the geopotential up to ultra-high degrees and orders: Efficient computation at irregular surfaces. Computers & Geosciences 66, 219-227,

Please use this reference in your works when using isGrafLab.

Go to top

GOCE-only global gravity field models from the kinematic orbit

Description: GOCE-only spherical harmonic models of the Earth's gravitational field

The models were estimated from the kinematic orbit of the GOCE satellite (01 November 2009 to 11 January 2010, the R1 period). Further details:

The effect caused by the stabilization of the near-zonal coefficients is shown in the figure below (Bucha et al., 2015).

When using the models, please use the following reference:

Bucha, B., Bezděk, A., Sebera, J., Janák, J., 2015. Global and regional gravity field determination from GOCE kinematic orbit by means of spherical radial basis functions. Surveys in Geophysics 36, 773-801,

Go to top

Divergence-free spherical harmonic gravity models of the Moon's topographic masses

Description: Spherical harmonic models of the gravitational field implied by the Moon's topographic masses

These models approximate the gravitational field induced by topographic masses of the Moon and mitigate the divergence effect of spherical harmonics in the vicinity of the Moon's topography. The models rely on the Runge—Krarup theorem and enable generally a more accurate evaluation of gravity field quantities close to the lunar topography as compared to the models from spectral gravity forward modelling. This is because the latter ones may suffer from the divergence effect when evaluating the spherical harmonic series on or below the limit sphere encompassing all the gravitating masses (that is, also on the topography). In the figure below, the two approaches are validated with respect to a third (reference) one (Bucha et al., 2019). The upper map shows the reference signal, the middle one depicts the differences related to the Runge—Krarup solution and the bottom one validates the spectral gravity forward modelling technique; note the different colour bars.

We provide 12 models that approximate the gravity field implied by the lunar topographic masses. The topographic masses are expanded up to degrees 90, 180, 360 and 720, and the maximum degree of the gravity field models varies from 360 up to 2160. All models are available in the gfc format as defined by ICGEM. One of the models, "STU_Moon_topography_to720_gravity_to2160", can also be accessed from ICGEM, where it can be find under a shortened name "STU_MoonTopo720".

The corresponding models that may suffer from the divergence effect were derived by Hirt and Kuhn (2017) and can be downloaded here, including grids of reference gravity disturbances from spatial-domain Newtonian integration.

Please use the following reference when using the models in your works:

Bucha, B., Hirt, C., Kuhn, M., 2019. Divergence-free spherical harmonic gravity field modelling based on the Runge—Krarup theorem: a case study for the Moon. Journal of Geodesy 93, 489-513,

Go to top

SRTM2gravity — Near-global 3 arcsec maps of topography-implied gravity

Name: SRTM2gravity
Description: Near-global 3 arcsec maps of topography-implied gravity
Authors: Christian Hirt et al.
Download: or
Academic Torrent:

For details on SRTM2gravity, visit

Go to top

CHarm — C/Python library for spherical harmonic transforms

Name: CHarm
Description: C/Python library for high-degree spherical harmonic transforms
Programming Language: C, Python
Source Code:

CHarm is a C library to work with spherical harmonics up to almost arbitrarily high degrees. The library is accompanied by a Python wrapper called PyHarm.


The predecessor of CHarm, the SHA_SHS package written in C, Fortran and MATLAB, is archived here. The package works correctly, has its own brief documentation, but is no longer maintained.

Go to top

Molodensky's truncation coefficients for cap-modified spectral gravity forward modelling

Description: Molodensky's truncation coefficients for cap integration in spectral gravity forward modelling
Download coefficients:
Download source code:

Cap-modified spectral gravity forward modelling is a spectral technique to deliver near- and/or far-zone gravity effects implied by topographic masses, the shape of which is given as a surface spherical harmonic series. From the numerical point of view, one of its most difficult parts is the accurate computation of the associated Molodensky's truncation coefficients. This is because this evaluation may necessitate to extend the number of significant digits from 16 (double precision) to, say, 64 or even 256. Here, this was achieved using the MATLAB's ADVANPIX toolbox.

MATLAB routines to compute the truncation coefficients and their radial and horizontal derivatives are also available from the link above. The ADVANPIX toolbox is required to run the codes.

Available are also some pre-computed ready-to-use near-zone and far-zone truncation coefficients from the Bucha et al. (2019) study. The coefficients are evaluated for

256 significant digits were used to compute the truncation coefficients, ensuring 24-digit accuracy or better. After the evaluation, the coefficients were converted to double precision with 16 significant digits. Importantly, in some cases, the loss of significance errors may be encountered during the spherical harmonic synthesis when using the coefficients (see the reference below).

Please use the following reference when using either the codes or the coefficients:

Bucha, B., Hirt, C., Kuhn, M., 2019. Cap integration in spectral gravity forward modelling up to the full gravity tensor. Journal of Geodesy,

Go to top


Name: SRBF_package
Description: Routines for gravity field modelling with spherical radial basis functions
Programming Languages: Fortran (+ Python wrapper), MATLAB

SRBF_package is a collection of a few routines for gravity field modelling with band-limited spherical radial basis functions (SRBFs), allowing both analysis and synthesis of the gravitational field. Harmonically analysed can be the gravitational potential (or any other continuous scalar function on a 2D sphere), while synthesis is possible for the gravitational potential, the gravitational vector in LNOF and the gravitational tensor in LNOF. To our knowledge, the equations to synthesize the gravitational tensor are novel and are more efficient than other formulae (e.g., that from Appendix of Bucha et al., 2016). This is because only three radial basis functions are needed to compute six unique elements of the gravitational tensor. This may improve the computational speed, as the evaluation of band-limited spherical radial basis functions (the sum over harmonic degrees) is the most time-consuming part of the synthesis.

As an example of a band-limited SRBF, the figure below depicts the smoothed Shannon SRBF.

The package is an enhanced version of the MATLAB-based routines used by Bucha et al. (2016). It is available as MATLAB-based functions and Fortran 95 routines accompanied by a Python wrapper. The MATLAB synthesis is parallelized using "parfor" and the Fortran synthesis and analysis are parallelized using OpenMP (the "!$OMP PARALLEL DO" directive). When using the same number of CPUs, the Fortran routines can be up to about 10 times faster than their MATLAB-counterparts. The Python wrapper allows an easy access to the Fortran routines for Python users.

Attached to the package are also also test data, manuals with a detailed description of the package and scripts with two test computations to verify the correctness of the implementation in MATLAB, Fortran and Python. While the test examples are based on global gravity field modelling for clarity, SRBF_package was developed primarily for regional applications with a residual gravity field. This is because SRBF_package is significantly slower when compared to FFT-based SHS and SHA, both of which are designed for global applications. Although the package was designed mainly for gravity field modelling, other applications are possible. Further details can be found in the attached PDF.

Please use the following reference when using SRBF_package:

Bucha, B., Janák, J., Papčo, J., Bezděk, A., 2016. High-resolution regional gravity field modelling in a mountainous area from terrestrial gravity data. Geophysical Journal International 207, 949-966,

Go to top


GMSQ2019 is a Gravimetric Model of the Slovak Quasigeoid at 2 arc-sec spatial resolution. It is developed from a high-resolution regional gravity field model (Bucha et al. 2016), which relies on gravity information from (i) EIGEN-6C4 up to degree 2190 (~9 km resolution), (ii) dense terrestrial gravity data modelled by spherical radial basis functions up to degree 21600 (~30 arc-sec), and (iii) residual terrain modelling (~2 arc-sec).

A validation of GMSQ2019 based on 34 high-quality GNSS/levelling data revealed a standard deviation of 1.9 cm. When including reference data of worse quality (but higher quantity, in total 8838), the standard deviation reaches 4.3 cm and a systematic bias of -0.5 m is observed.

The development of GMSQ2019 is discussed in detail in a technical report (in Slovak), which is freely available at (~17 MB). Further technicalities are provided in the Bucha et al. (2016) paper. GMSQ2019 in 6 arc-sec spatial resolution is shown in the figure below (units are metres, figure size is ~2.8 MB). A figure at the original 2 arc-sec resolution is provided in the technical report (Obr. 1).

Go to top


Name: grav-sr-2arcsec
Description: 2-arcsec surface gravitational maps of the Slovak Republic
Version: 1.0

grav-sr-2arcsec is a suite of 2-arcsec surface gravitational maps of the Slovak Republic in terms of

Below is shown the $V_{zz}$ element of the gravitational tensor over the High and Low Tatra Mountains, representing the roughest part of Slovakia.

The next figure shows the $V_y$ element of the gravitational vector over the more or less same region.

High-resolution figures of the entire grids can be found for all ten quantities in the "figs" directory at link above (more than 30 MB per figure).

The suite was computed from a high-resolution gravity field model of Slovakia developed by Bucha et al. (2016). The maps rely on three sources of gravitational information:

An independent validation of grav-sr-2arcsec revealed the standard deviation of 0.789 mGal in terms of the gravity.

Since the residual terrain modelling technique is based on the constant-mass density assumption, grav-sr-2arcsec must *not* be geophysically interpreted at resolutions higher than ~30 arcsec (~950 m). This is because this information comes solely from the shape of the topography and a constant mass density, so the effect of anomalous densities is *not* captured beyond the ~30 arcsec resolution. In the case of the combination "EIGEN-6C4" + "residual terrain modelling" (that is, no spherical radial basis functions), the geophysically relevant information can be found only up to the resolution of EIGEN-6C4, that is, ~9 km. Nevertheless, the residual terrain modelling grids are still very useful for improving the accuracy with respect to independent ground-truth (e.g., observed terrestrial gravity).

Unfortunately, due to the license restrictions, the gravitational information captured by spherical radial basis functions is not allowed to be published as user-friendly grids. The information is, however, provided in terms of expansion coefficients of spherical radial basis functions, so that one can *exactly* reproduce their contribution, avoiding any approximations (attached is an efficient ready-to-use Fortran program that will do the work for you).

To name a few possible applications:

If you use the grav-sr-2arcsec suite, I will be happy if you let me know about your application!

The core data (~6.5 GB) are provided as HDF5 data files, which can be read in many commonly used programming languages. Further informations on grav-sr-2arcsec can be found in the attached README.txt file.

Go to top

GFM_DE_rule package

Name: GFM_DE_rule
Description: Spatial-domain gravity forward modelling
Programming Language: Fortran, MATLAB

The GFM_DE_rule package computes the gravitational potential of a body, the shape of which is given by a surface spherical harmonic expansion and the density of which is constant. The evaluation points can be located inside the body, on its surface or outside the body.

The package builds on the spatial-domain gravity forward modelling technique developed by Fukushima (2017). It evaluates the Newton integral in the spatial domain via the double exponential rule. The Fortran routines from that paper are utilized, while several modifications were necessary to fit the goals of the package.

The package is written in Fortran and MATLAB. The Fortran version supports computations in double as well as quadruple precision and is generally faster than the MATLAB version.

The accuracy of the output potential can reach 14 correct digits in double precision and ~30 digits in quadruple precision. The accuracy level is controlled via a relative tolerance error parameter. The higher the required accuracy, the longer the computation times can be expected. The package was designed mainly for modelling gravitational fields of bodies with shapes expanded to low-degree surface spherical harmonic series.

As an example, below are shown respectively degree-15 surface spherical harmonic expansion of shape of the asteroid Bennu (units are metres) and the implied surface gravitational potential ($\mathrm{m}^2 \, \mathrm{s}^{-2}$). At 99.4 % of the grid points, the gravitational potential is computed correctly up to the 14th digit. The figures are due to Bucha and Sanso (2021). The data used to prepare the figures can be found in the Asteroid (101955) Bennu and its gravitational field section.

Go to top

Asteroid Bennu and its gravitational field

Description: Dataset from the Bucha and Sanso (2021) study
Dataset Download:

In Bucha and Sanso (2021), we investigated three spherical-harmonic-based techniques to deliver external gravitational field models of the asteroid (101955) Bennu within its circumscribing sphere. This region is known to be peculiar for external spherical harmonic expansions, because it may lead to a divergent series. The studied approaches are

While the first method diverges beyond any reasonable doubts, we show that the other two methods may ensure relative accuracy from ${\sim}10^{-6}$ to $10^{-8}$ in the vicinity of Bennu. This is possible even with the second method, despite the fact that it relies on a single series of external spherical harmonics. Our main motivation was to study conceptual differences between spherical harmonic coefficients from satellite data (analogy to the first method) and from surface gravitational data (the second method). Such coefficients are available through the popular spherical-harmonic-based models of the Earth's gravitational field and often are combined together. We show that the coefficients from terrestrial data may lead to a divergence effect of partial sums, though excellent accuracy can be achieved when such model is used in full. Under (presently) extreme but realistic conditions, the divergence effect of partial sums may affect many near-surface geoscientific applications, such as the geoid/quasigeoid computation or residual terrain modelling.

The dataset is available at the link above. Next follows a brief description.

Surface spherical harmonic coefficients of Bennu's shape

The shape surface spherical harmonic coefficients of Bennu are provided up to maximum harmonic degree 15 and were computed from a polyhedral shape model. The physical model of Bennu was obtained after assigning a single-constant mass density of $\rho = 1260\ \mathrm{kg} \, \mathrm{m}^{-3}$ to this shape. A figure of the synthesized shape of Bennu is provided in the GFM_DE_rule section above. The original polyhedral model is shown below (colour bar in meters, Bucha and Sanso, 2021).

Spectral gravity forward modelling using external spherical harmonics

Using the degree-15 shape expansion of Bennu and assuming a constant-mass density, we computed spherical harmonic coefficients of the external gravitational potential using spectral gravity forward modelling. The maximum degree of the potential coefficients is 225. On a Brillouin sphere, they reproduce the gravitational potential with a 14-digit accuracy or better. However, on/below the sphere of convergence, including the surface of Bennu, the spherical harmonic series diverges and produces grossly invalid results. The figure below shows the standard deviations of the errors of the model as a function of maximum harmonic degree (Bucha and Sanso, 2021).

Least-squares estimation from surface gravitational data using external spherical harmonics

To allow for an accurate spherical-harmonic-based gravitational field modelling everywhere on and above the surface of Bennu while still relying on external spherical harmonics, we developed models from surface gravitational data using the least-squares estimation. These models rely on the Runge—Krarup theorem and provide an $\varepsilon$-accuracy, $\varepsilon > 0$, everywhere on and above Bennu. In our case, they outperform spectral gravity forward modelling on the surface of Bennu by more than 3 orders of magnitude (relative accuracy is ${\sim}10^{-6}$; Bucha and Sanso, 2021). The models were estimated up to maximum degrees 0, 5, 10, ..., 65. The figure below validates them on the surface of Bennu; the solid colour lines represent standard deviations of the errors of the least-squares models and the dashed black line represents spectral gravity forward modelling.

The computation of the models took 90 CPU years.

Combination of internal and external series expansions

Finally, spherical harmonic coefficients of the internal and external series expansions were computed. They do not suffer from the divergence, but their usage on the surface of Bennu is less straightforward as with external spherical harmonics. The figure below (Bucha and Sanso, 2021) shows the accuracy levels achieved with spectral gravity forward modelling (the orange curve), least-squares method (green) and the combination of internal and external series expansions (blue) as a function of harmonic degree.

Reference gravitational data

Reference gravitational data (accurate up to the 14th digit) on a Brillouin sphere, a mean sphere and the surface of Bennu are also available. The data were computed using the GFM_DE_rule package. The surface potential is shown in the GFM_DE_rule section above.

Go to top


Preprints are available at my ResearchGate profile. For some of the papers above, free full-text view-only versions are available.

Go to top


Go to top


Go to top


Go to top

Favourite websites

Go to top

Last modified: .