How to use the harmonica.constants.GRAVITATIONAL_CONST function in harmonica

result = np.zeros(cast.size, dtype=dtype)
    # Prepare arrays to be passed to the jitted functions
    coordinates = tuple(np.atleast_1d(i).ravel() for i in coordinates[:3])
    points = tuple(np.atleast_1d(i).ravel() for i in points[:3])
    masses = np.atleast_1d(masses).ravel()
    # Sanity checks
    if masses.size != points[0].size:
        raise ValueError(
            "Number of elements in masses ({}) ".format(masses.size)
            + "mismatch the number of points ({})".format(points[0].size)
        )
    # Compute gravitational field
    dispatchers[coordinate_system](
        *coordinates, *points, masses, result, kernels[coordinate_system][field]
    )
    result *= GRAVITATIONAL_CONST
    # Convert to more convenient units
    if field in ("g_easting", "g_northing", "g_z"):
        result *= 1e5  # SI to mGal
    return result.reshape(cast.shape)

coordinates,
        tesseroids,
        density,
        stack,
        small_tesseroids,
        point_masses,
        weights,
        result,
        distance_size_ratio,
        radial_adaptive_discretization,
        n_nodes,
        glq_nodes,
        glq_weights,
        kernels[field],
    )
    result *= GRAVITATIONAL_CONST
    # Convert to more convenient units
    if field == "g_z":
        result *= 1e5  # SI to mGal
    return result.reshape(cast.shape)

result = np.zeros(cast.size, dtype=dtype)
    # Convert coordinates, prisms and density to arrays with proper shape
    coordinates = tuple(np.atleast_1d(i).ravel() for i in coordinates[:3])
    prisms = np.atleast_2d(prisms)
    density = np.atleast_1d(density).ravel()
    # Sanity checks
    if not disable_checks:
        if density.size != prisms.shape[0]:
            raise ValueError(
                "Number of elements in density ({}) ".format(density.size)
                + "mismatch the number of prisms ({})".format(prisms.shape[0])
            )
        _check_prisms(prisms)
    # Compute gravitational field
    jit_prism_gravity(coordinates, prisms, density, kernels[field], result)
    result *= GRAVITATIONAL_CONST
    # Convert to more convenient units
    if field == "g_z":
        result *= 1e5  # SI to mGal
    return result.reshape(cast.shape)

-------
    grav_bouguer : array or :class:`xarray.DataArray`
        The gravitational effect of topography and residual bathymetry in mGal.

    """
    # Need to cast to array to make sure numpy indexing works as expected for
    # 1D DataArray topography
    oceans = np.array(topography < 0)
    continent = np.logical_not(oceans)
    density = np.full(topography.shape, np.nan, dtype="float")
    density[continent] = density_crust
    # The minus sign is used to negate the bathymetry (which is negative and
    # the equation calls for "thickness", not height). This is more practical
    # than taking the absolute value of the topography.
    density[oceans] = -1 * (density_water - density_crust)
    bouguer = 1e5 * 2 * np.pi * GRAVITATIONAL_CONST * density * topography
    return bouguer