'''
This code presents an approach for implementing the gravitational field
produced by a sphere by using the analytical formulas.

References

'''

import numpy as np
from numba import jit

#: The gravitational constant in m^3 kg^{-1} s^{-1}
GRAVITATIONAL_CONST = 0.00000000006673

def gravity_gradient(coordinates,spheres,density,field):
    """
    Computing the gravitational field generated by a sphere.
    """

    kernels = {"g_xx":kernel_g_xx,
               "g_yy":kernel_g_yy,
               "g_zz":kernel_g_zz,
               "g_xy":kernel_g_xy,
               "g_yz":kernel_g_yz,
               "g_xz":kernel_g_xx
               }

    # Checking input field
    if field not in kernels:
        raise ValueError("Gravitational field {} not recognized".format(field))

    # Checking input parameters
    if coordinates.ndim != 2:
        raise ValueError(
            "coordinates ndim ({}) ".format(coordinates.ndim)
            + "not equal to 2"
        )
    if coordinates.shape[0] != 3:
        raise ValueError(
            "Number of lines in coordinates ({}) ".format(coordinates.shape[0])
            + "not equal to 3"
        )
    if spheres.ndim != 2:
        raise ValueError(
            "prisms ndim ({}) ".format(spheres.ndim)
            + "not equal to 2"
        )
    if spheres.shape[1] != 4:
        raise ValueError(
            "Number of columns in spheres ({}) ".format(spheres.shape[1])
            + "not equal to 6"
        )
    if density.ndim != 1:
        raise ValueError(
            "density ndim ({}) ".format(density.ndim)
            + "not equal to 1"
        )
    if density.size != spheres.shape[0]:
        raise ValueError(
            "Number of elements in density ({}) ".format(density.size)
            + "mismatch the number of spheres ({})".format(spheres.shape[0])
        )

    # Create the array to store the result
    result = np.zeros(coordinates[0].size, dtype="float64")

    # Compute gravitational field
    jit_gravitational(coordinates, spheres, density, kernels[field], result)
    result *= GRAVITATIONAL_CONST
    # Convert from SI to Eotvos
    if field in ["g_xx", "g_yy","g_zz","g_xy","g_xz","g_yz"]:
        result *= 1e9
    return result

@jit(nopython=True)
def jit_gravitational(coordinates,spheres,density,kernel,out):
    """
    Compute gravitational field at the observation points
    """

    # Iterate over coordinates
    for l in range(coordinates[0].size):
        # Iterate over spheres:
        for m in range(spheres.shape[0]):
            Y = spheres[m,0] - coordinates[0,l]
            X = spheres[m,1] - coordinates[1,l]
            Z = spheres[m,2] - coordinates[2,l]
            r3 = spheres[m,3]*spheres[m,3]*spheres[m,3]
            Volume = (4/3)*np.pi*r3
            out[l] += (density[m]*Volume*kernel(Y,X,Z))

@jit(nopython=True)
def kernel_g_xx(Y, X, Z):
    """
    Kernel for xx component of gravitational acceleration of a sphere
    """
    radius = np.sqrt(Y ** 2 + X ** 2 + Z ** 2)
    r2 = radius*radius
    r5 = r2*r2*radius
    kernel = (3*X**2 - r2)/r5
    return kernel

@jit(nopython=True)
def kernel_g_yy(Y, X, Z):
    """
    Kernel for yy component of gravitational acceleration of a sphere
    """
    radius = np.sqrt(Y ** 2 + X ** 2 + Z ** 2)
    r2 = radius*radius
    r5 = r2*r2*radius
    kernel = (3*Y**2 - r2)/r5
    return kernel

@jit(nopython=True)
def kernel_g_zz(Y, X, Z):
    """
    Kernel for xx component of gravitational acceleration of a sphere
    """
    radius = np.sqrt(Y ** 2 + X ** 2 + Z ** 2)
    r2 = radius*radius
    r5 = r2*r2*radius
    kernel = (3*Z**2 - r2)/r5
    return kernel

@jit(nopython=True)
def kernel_g_xy(Y, X, Z):
    """
    Kernel for xy component of gravitational acceleration of a sphere
    """
    radius = np.sqrt(Y ** 2 + X ** 2 + Z ** 2)
    r2 = radius*radius
    r5 = r2*r2*radius
    kernel = (3*X*Y)/r5
    return kernel

@jit(nopython=True)
def kernel_g_yz(Y, X, Z):
    """
    Kernel for xz component of gravitational acceleration of a sphere
    """
    radius = np.sqrt(Y ** 2 + X ** 2 + Z ** 2)
    r2 = radius*radius
    r5 = r2*r2*radius
    kernel = (3*Y*Z)/r5
    return kernel

'''
This code presents an approach for implementing the gravitational field
produced by a sphere by using the analytical formulas.

References

'''

import numpy as np
from numba import jit

#: The gravitational constant in m^3 kg^{-1} s^{-1}
GRAVITATIONAL_CONST = 0.00000000006673

def gravity_gradient(coordinates,spheres,density,field):
    """
    Computing the gravitational field generated by a sphere.
    """

    kernels = {"g_xx":kernel_g_xx,
               "g_yy":kernel_g_yy,
               "g_zz":kernel_g_zz,
               "g_xy":kernel_g_xy,
               "g_yz":kernel_g_yz,
               "g_xz":kernel_g_xx
               }

    # Checking input field
    if field not in kernels:
        raise ValueError("Gravitational field {} not recognized".format(field))

    # Checking input parameters
    if coordinates.ndim != 2:
        raise ValueError(
            "coordinates ndim ({}) ".format(coordinates.ndim)
            + "not equal to 2"
        )
    if coordinates.shape[0] != 3:
        raise ValueError(
            "Number of lines in coordinates ({}) ".format(coordinates.shape[0])
            + "not equal to 3"
        )
    if spheres.ndim != 2:
        raise ValueError(
            "prisms ndim ({}) ".format(spheres.ndim)
            + "not equal to 2"
        )
    if spheres.shape[1] != 4:
        raise ValueError(
            "Number of columns in spheres ({}) ".format(spheres.shape[1])
            + "not equal to 6"
        )
    if density.ndim != 1:
        raise ValueError(
            "density ndim ({}) ".format(density.ndim)
            + "not equal to 1"
        )
    if density.size != spheres.shape[0]:
        raise ValueError(
            "Number of elements in density ({}) ".format(density.size)
            + "mismatch the number of spheres ({})".format(spheres.shape[0])
        )

    # Create the array to store the result
    result = np.zeros(coordinates[0].size, dtype="float64")

    # Compute gravitational field
    jit_gravitational(coordinates, spheres, density, kernels[field], result)
    result *= GRAVITATIONAL_CONST
    # Convert from SI to Eotvos
    if field in ["g_xx", "g_yy","g_zz","g_xy","g_xz","g_yz"]:
        result *= 1e9
    return result

@jit(nopython=True)
def jit_gravitational(coordinates,spheres,density,kernel,out):
    """
    Compute gravitational field at the observation points
    """

    # Iterate over coordinates
    for l in range(coordinates[0].size):
        # Iterate over spheres:
        for m in range(spheres.shape[0]):
            Y = spheres[m,0] - coordinates[0,l]
            X = spheres[m,1] - coordinates[1,l]
            Z = spheres[m,2] - coordinates[2,l]
            r3 = spheres[m,3]*spheres[m,3]*spheres[m,3]
            Volume = (4/3)*np.pi*r3
            out[l] += (density[m]*Volume*kernel(Y,X,Z))

@jit(nopython=True)
def kernel_g_xx(Y, X, Z):
    """
    Kernel for xx component of gravitational acceleration of a sphere
    """
    radius = np.sqrt(Y ** 2 + X ** 2 + Z ** 2)
    r2 = radius*radius
    r5 = r2*r2*radius
    kernel = (3*X**2 - r2)/r5
    return kernel

@jit(nopython=True)
def kernel_g_yy(Y, X, Z):
    """
    Kernel for yy component of gravitational acceleration of a sphere
    """
    radius = np.sqrt(Y ** 2 + X ** 2 + Z ** 2)
    r2 = radius*radius
    r5 = r2*r2*radius
    kernel = (3*Y**2 - r2)/r5
    return kernel

@jit(nopython=True)
def kernel_g_zz(Y, X, Z):
    """
    Kernel for xx component of gravitational acceleration of a sphere
    """
    radius = np.sqrt(Y ** 2 + X ** 2 + Z ** 2)
    r2 = radius*radius
    r5 = r2*r2*radius
    kernel = (3*Z**2 - r2)/r5
    return kernel

@jit(nopython=True)
def kernel_g_xy(Y, X, Z):
    """
    Kernel for xy component of gravitational acceleration of a sphere
    """
    radius = np.sqrt(Y ** 2 + X ** 2 + Z ** 2)
    r2 = radius*radius
    r5 = r2*r2*radius
    kernel = (3*X*Y)/r5
    return kernel

@jit(nopython=True)
def kernel_g_yz(Y, X, Z):
    """
    Kernel for xz component of gravitational acceleration of a sphere
    """
    radius = np.sqrt(Y ** 2 + X ** 2 + Z ** 2)
    r2 = radius*radius
    r5 = r2*r2*radius
    kernel = (3*Y*Z)/r5
    return kernel

This code presents an approach for implementing the gravitational field produced by a sphere by using the analytical formulas.  References #: The gravitational constant in m^3 kg^{-1} s^{-1} Computing the gravitational field generated by a sphere. # Checking input field # Checking input parameters # Create the array to store the result # Compute gravitational field # Convert from SI to Eotvos Compute gravitational field at the observation points # Iterate over coordinates # Iterate over spheres: Kernel for xx component of gravitational acceleration of a sphere Kernel for yy component of gravitational acceleration of a sphere Kernel for xx component of gravitational acceleration of a sphere Kernel for xy component of gravitational acceleration of a sphere Kernel for xz component of gravitational acceleration of a sphere