<commit_before>r"""
Calculate the potential fields of a homogeneous sphere.

**Magnetic**

Calculates the magnetic effect produced by an sphere. The functions are
based on Blakely (1995).

* :func:`~fatiando.gravmag.sphere.tf`: calculates the total-field anomaly
* :func:`~fatiando.gravmag.sphere.bx`: calculates the x component of the
  induction
* :func:`~fatiando.gravmag.sphere.by`: calculates the y component of the
  induction
* :func:`~fatiando.gravmag.sphere.bz`: calculates the z component of the
  induction

Remember that:

The magnetization :math:`\mathbf{M}` and the dipole moment :math:`\mathbf{m}`
are related with the volume V:

.. math::

    \mathbf{M} = \dfrac{\mathbf{m}}{V}.

The total-field anomaly is:

.. math::

    \Delta T = |\mathbf{T}| - |\mathbf{F}|,

where :math:`\mathbf{T}` is the measured field and :math:`\mathbf{F}` is a
reference (regional) field. The forward modeling functions
:func:`~fatiando.gravmag.sphere.bx`, :func:`~fatiando.gravmag.sphere.by`,
and :func:`~fatiando.gravmag.sphere.bz` calculate the 3 components of the
field perturbation :math:`\Delta\mathbf{F}`

.. math::

    \Delta\mathbf{F} = \mathbf{T} - \mathbf{F}.

Then the total-field anomaly caused by the sphere is

.. math::

    \Delta T \approx \hat{\mathbf{F}}\cdot\Delta\mathbf{F}.

**Gravity**

Calculates the gravitational acceleration and gravity gradient tensor
components.

* :func:`~fatiando.gravmag.sphere.gz`
* :func:`~fatiando.gravmag.sphere.gxx`
* :func:`~fatiando.gravmag.sphere.gxy`
* :func:`~fatiando.gravmag.sphere.gxz`
* :func:`~fatiando.gravmag.sphere.gyy`
* :func:`~fatiando.gravmag.sphere.gyz`
* :func:`~fatiando.gravmag.sphere.gzz`


**Auxiliary Functions**

Calculates the second derivatives of the function

.. math::

    \phi(x,y,z) = \frac{4}{3} \pi R^3 \frac{1}{r}

with respect to the variables :math:`x`, :math:`y`, and :math:`z`. In
this equation,

.. math::

    r = \sqrt{(x - \nu)^2 + (y - \eta)^2 + (z - \zeta)^2},

and :math:`R` is the radius of a sphere with centre at the Cartesian
coordinates :math:`\nu`, :math:`\eta` and :math:`\zeta`.

These second derivatives are used to calculate the total field magnetic anomaly
and the gravity gradient tensor components.

* :func:`~fatiando.gravmag.sphere.kernelxx`
* :func:`~fatiando.gravmag.sphere.kernelxy`
* :func:`~fatiando.gravmag.sphere.kernelxz`
* :func:`~fatiando.gravmag.sphere.kernelyy`
* :func:`~fatiando.gravmag.sphere.kernelyz`
* :func:`~fatiando.gravmag.sphere.kernelzz`

**References**

<NAME>. (1995), Potential Theory in Gravity and Magnetic Applications,
Cambridge University Press.

----

"""
from __future__ import division

import numpy

from ..constants import SI2MGAL, G, CM, T2NT, SI2EOTVOS
from .. import utils


def tf(xp, yp, zp, spheres, inc, dec, pmag=None):
    """
    Calculate the total-field anomaly of spheres.

    .. note:: Input units are SI. Output is in nT

    Parameters:

    * xp, yp, zp : arrays
        The x, y, and z coordinates where the anomaly will be calculated
    * spheres : list of :class:`fatiando.mesher.Sphere`
        The spheres. Spheres must have the physical property
        ``'magnetization'``. Spheres without ``'magnetization'`` will be
        ignored.
    * inc : float
        The inclination of the regional field (in degrees)
    * dec : float
        The declination of the regional field (in degrees)
    * pmag : [mx, my, mz] or None
        A magnetization vector. If not None, will use this value instead of the
        ``'magnetization'`` property of the spheres. Use this, e.g., for
        sensitivity matrix building.

    Returns:

    * tf : array
        The total-field anomaly

    Example:

    >>> from fatiando import mesher, gridder, utils
    >>> # Set the inclination and declination of the regional field
    >>> inc, dec = -30, 45
    >>> # Create a sphere model
    >>> model = [
    ...         # One with induced magnetization
    ...         mesher.Sphere(1000, 1000, 600, 500, {'magnetization':5}),
    ...         # and one with remanent
    ...         mesher.Sphere(-1000, -1000, 600, 500,
    ...             {'magnetization':utils.ang2vec(10, 70, -5)})]
    >>> # Create a regular grid at 100m height
    >>> shape = (4, 4)
    >>> area = (-3000, 3000, -3000, 3000)
    >>> xp, yp, zp = gridder.regular(area, shape, z=-100)
    >>> # Calculate the anomaly for a given regional field
    >>> for t in tf(xp, yp, zp, model, inc, dec):
    ...     print '%15.8e' % t
     2.72951375e+01
     3.63637351e+01
     5.35842876e+00
     6.65189557e-02
     6.01998831e+01
    -1.71920499e+03
     3.30025228e+00
    -5.13176612e+00
     1.47871812e+00
    -3.96103758e+01
    -1.85654021e+02
     2.18002960e+01
    -2.82713826e+00
    -1.10341542e+01
     1.93353982e+01
     2.03174254e+01

    """
    if xp.shape != yp.shape != zp.shape:
        raise ValueError("Input arrays xp, yp, and zp must have same shape!")
    size = len(xp)
    res = numpy.zeros(size, dtype=numpy.float)
    # Calculate the 3 components of the unit vector in the direction of the
    # regional field
    fx, fy, fz = utils.dircos(inc, dec)
    if pmag is not None:
        if isinstance(pmag, float) or isinstance(pmag, int):
            pmx, pmy, pmz = pmag*fx, pmag*fy, pmag*fz
        else:
            pmx, pmy, pmz = pmag
    for sphere in spheres:
        if sphere is None or ('magnetization' not in sphere.props
                              and pmag is None):
            continue
        radius = sphere.radius
        # Get the intensity and unit vector from the magnetization
        if pmag is None:
            mag = sphere.props['magnetization']
            if isinstance(mag, float) or isinstance(mag, int):
                mx, my, mz = mag*fx, mag*fy, mag*fz
            else:
                mx, my, mz = mag
        else:
            mx, my, mz = pmx, pmy, pmz
        v1 = kernelxx(xp, yp, zp, sphere)
        v2 = kernelxy(xp, yp, zp, sphere)
        v3 = kernelxz(xp, yp, zp, sphere)
        v4 = kernelyy(xp, yp, zp, sphere)
        v5 = kernelyz(xp, yp, zp, sphere)
        v6 = kernelzz(xp, yp, zp, sphere)
        bx = (v1*mx + v2*my + v3*mz)
        by = (v2*mx + v4*my + v5*mz)
        bz = (v3*mx + v5*my + v6*mz)
        res += fx*bx + fy*by + fz*bz
    res *= CM*T2NT
    return res

def bx(xp, yp, zp, spheres):
    """
    Calculates the x component of the magnetic induction produced by spheres.

    .. note:: Input units are SI. Output is in nT

    Parameters:

    * xp, yp, zp : arrays
        The x, y, and z coordinates where the anomaly will be calculated
    * spheres : list of :class:`fatiando.mesher.Sphere`
        The spheres. Spheres must have the physical property<commit_msg>Revert sphere.tf to old implementation<commit_after>r"""
Calculate the potential fields of a homogeneous sphere.

**Magnetic**

Calculates the magnetic effect produced by an sphere. The functions are
based on Blakely (1995).

* :func:`~fatiando.gravmag.sphere.tf`: calculates the total-field anomaly
* :func:`~fatiando.gravmag.sphere.bx`: calculates the x component of the
  induction
* :func:`~fatiando.gravmag.sphere.by`: calculates the y component of the
  induction
* :func:`~fatiando.gravmag.sphere.bz`: calculates the z component of the
  induction

Remember that:

The magnetization :math:`\mathbf{M}` and the dipole moment :math:`\mathbf{m}`
are related with the volume V:

.. math::

    \mathbf{M} = \dfrac{\mathbf{m}}{V}.

The total-field anomaly is:

.. math::

    \Delta T = |\mathbf{T}| - |\mathbf{F}|,

where :math:`\mathbf{T}` is the measured field and :math:`\mathbf{F}` is a
reference (regional) field. The forward modeling functions
:func:`~fatiando.gravmag.sphere.bx`, :func:`~fatiando.gravmag.sphere.by`,
and :func:`~fatiando.gravmag.sphere.bz` calculate the 3 components of the
field perturbation :math:`\Delta\mathbf{F}`

.. math::

    \Delta\mathbf{F} = \mathbf{T} - \mathbf{F}.

Then the total-field anomaly caused by the sphere is

.. math::

    \Delta T \approx \hat{\mathbf{F}}\cdot\Delta\mathbf{F}.

**Gravity**

Calculates the gravitational acceleration and gravity gradient tensor
components.

* :func:`~fatiando.gravmag.sphere.gz`
* :func:`~fatiando.gravmag.sphere.gxx`
* :func:`~fatiando.gravmag.sphere.gxy`
* :func:`~fatiando.gravmag.sphere.gxz`
* :func:`~fatiando.gravmag.sphere.gyy`
* :func:`~fatiando.gravmag.sphere.gyz`
* :func:`~fatiando.gravmag.sphere.gzz`


**Auxiliary Functions**

Calculates the second derivatives of the function

.. math::

    \phi(x,y,z) = \frac{4}{3} \pi R^3 \frac{1}{r}

with respect to the variables :math:`x`, :math:`y`, and :math:`z`. In
this equation,

.. math::

    r = \sqrt{(x - \nu)^2 + (y - \eta)^2 + (z - \zeta)^2},

and :math:`R` is the radius of a sphere with centre at the Cartesian
coordinates :math:`\nu`, :math:`\eta` and :math:`\zeta`.

These second derivatives are used to calculate the total field magnetic anomaly
and the gravity gradient tensor components.

* :func:`~fatiando.gravmag.sphere.kernelxx`
* :func:`~fatiando.gravmag.sphere.kernelxy`
* :func:`~fatiando.gravmag.sphere.kernelxz`
* :func:`~fatiando.gravmag.sphere.kernelyy`
* :func:`~fatiando.gravmag.sphere.kernelyz`
* :func:`~fatiando.gravmag.sphere.kernelzz`

**References**

<NAME>. (1995), Potential Theory in Gravity and Magnetic Applications,
Cambridge University Press.

----

"""
from __future__ import division

import numpy

from ..constants import SI2MGAL, G, CM, T2NT, SI2EOTVOS
from .. import utils


def tf(xp, yp, zp, spheres, inc, dec, pmag=None):
    """
    Calculate the total-field anomaly of spheres.

    .. note:: Input units are SI. Output is in nT

    Parameters:

    * xp, yp, zp : arrays
        The x, y, and z coordinates where the anomaly will be calculated
    * spheres : list of :class:`fatiando.mesher.Sphere`
        The spheres. Spheres must have the physical property
        ``'magnetization'``. Spheres without ``'magnetization'`` will be
        ignored.
    * inc : float
        The inclination of the regional field (in degrees)
    * dec : float
        The declination of the regional field (in degrees)
    * pmag : [mx, my, mz] or None
        A magnetization vector. If not None, will use this value instead of the
        ``'magnetization'`` property of the spheres. Use this, e.g., for
        sensitivity matrix building.

    Returns:

    * tf : array
        The total-field anomaly

    Example:

    >>> from fatiando import mesher, gridder, utils
    >>> # Set the inclination and declination of the regional field
    >>> inc, dec = -30, 45
    >>> # Create a sphere model
    >>> model = [
    ...         # One with induced magnetization
    ...         mesher.Sphere(1000, 1000, 600, 500, {'magnetization':5}),
    ...         # and one with remanent
    ...         mesher.Sphere(-1000, -1000, 600, 500,
    ...             {'magnetization':utils.ang2vec(10, 70, -5)})]
    >>> # Create a regular grid at 100m height
    >>> shape = (4, 4)
    >>> area = (-3000, 3000, -3000, 3000)
    >>> xp, yp, zp = gridder.regular(area, shape, z=-100)
    >>> # Calculate the anomaly for a given regional field
    >>> for t in tf(xp, yp, zp, model, inc, dec):
    ...     print '%15.8e' % t
     2.72951375e+01
     3.63637351e+01
     5.35842876e+00
     6.65189557e-02
     6.01998831e+01
    -1.71920499e+03
     3.30025228e+00
    -5.13176612e+00
     1.47871812e+00
    -3.96103758e+01
    -1.85654021e+02
     2.18002960e+01
    -2.82713826e+00
    -1.10341542e+01
     1.93353982e+01
     2.03174254e+01

    """
    if xp.shape != yp.shape != zp.shape:
        raise ValueError("Input arrays xp, yp, and zp must have same shape!")
    size = len(xp)
    res = numpy.zeros(size, dtype=numpy.float)
    # Calculate the 3 components of the unit vector in the direction of the
    # regional field
    fx, fy, fz = utils.dircos(inc, dec)
    if pmag is not None:
        if isinstance(pmag, float) or isinstance(pmag, int):
            pmx, pmy, pmz = pmag*fx, pmag*fy, pmag*fz
        else:
            pmx, pmy, pmz = pmag
    for sphere in spheres:
        if sphere is None or ('magnetization' not in sphere.props
                              and pmag is None):
            continue
        radius = sphere.radius
        # Get the intensity and unit vector from the magnetization
        if pmag is None:
            mag = sphere.props['magnetization']
            if isinstance(mag, float) or isinstance(mag, int):
                mx, my, mz = mag*fx, mag*fy, mag*fz
            else:
                mx, my, mz = mag
        else:
            mx, my, mz = pmx, pmy, pmz
        # First thing to do is make the computation point P the origin of the
        # coordinate system
        x = sphere.x - xp
        y = sphere.y - yp
        z = sphere.z - zp
        # Calculate the 3 components of B
        dotprod = mx*x + my*y + mz*z
        r_sqr = x**2 + y**2 + z**2
        r5 = r_sqr**(2.5)
        moment = 4.*numpy.pi*(radius**3)/3.
        bx = moment*(3*dotprod*x - r_sqr*mx)/r5
        by = moment*(3*dotprod*y - r_sqr*my)/r5
        bz = moment*(3*dotprod*z - r_sqr*mz)/r5
        #v1 = kernelxx(xp, yp, zp, sphere)
        #v2 = kernelxy(xp, yp, zp, sphere)
        #v3 = kernelxz(xp, yp, zp, sphere)
        #v4 = kernelyy(xp, yp, zp, sphere)
        #v5 = kernelyz(xp, yp, zp, sphere)
        #v6 = kernelzz(xp, yp, zp, sphere)
        #bx = (v1*mx + v2*my + v3*mz)
        #by = (v2*mx + v4*my + v5*mz)
        #bz = (v3*mx + v5*my + v6*mz)
        res += fx*bx + fy*by + fz*bz
    res *= CM*T2NT
    return res

def bx(xp, yp, zp, spheres):
    """
    Calculates the x component of the magnetic induction produced by spheres.

    .. note:: Input units are SI. Output is in nT

    Parameters:

    * xp, yp, zp : arrays
        The x, y, and z coordinates where the anomaly will be calculated
    * spheres : list of :class:`fatiando.mesher.Sphere`
        The spheres. Spheres must have the physical property
{'fraction_non_alphanumeric': 0.1335198737265031, 'fraction_numerical': 0.041038886497345384, 'mean_word_length': 2.6451359832635983, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 71, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '7289273', 'n_tokens_mistral': 5268, 'n_tokens_neox': 4774, 'n_words': 1449}