Generate random directions

magali/__init__.py

@@ -4,5 +4,5 @@
 #
 # This code is part of the Fatiando a Terra project (https://www.fatiando.org)
 #
-from ._input_output import read_qdm_harvard
+from ._input_output import random_unitary_vector, read_qdm_harvard
 from ._version import __version__

magali/_input_output.py

@@ -114,3 +114,68 @@ def _create_qdm_harvard_grid(coordinates, data_names, bz, path):
     qdm_data.bz.attrs = {"long_name": "vertical magnetic field", "units": "nT"}
     qdm_data.attrs = {"file_name": str(path)}
     return qdm_data
+
+
+def random_unitary_vector(dispersion_angle, loc=0, size=100):
+    """
+    Creates a random unitary vector from a defined dispersion angle.
+
+    Parameters
+    ----------
+    dispersion_angle: float
+        Dispersion angle that defines a region in a sphere surface.
+    loc:float or array_like of floats
+        Mean (“centre”) of the distribution.
+    size:int or tuple of ints, optional
+        Output shape. If the given shape is, e.g., (m, n, k), then m * n * k
+        samples are drawn. If size is None (default), a single value is
+        returned if loc and scale are both scalars. Otherwise,
+        np.broadcast(loc, scale).size samples are drawn.
+
+    Returns
+    -------
+    r_vector : :class:'numpy.ndarray'
+        A random unitary vector.
+    """
+    alpha = np.deg2rad(np.random.uniform(0, 360, size))
+    r = np.random.normal(loc, dispersion_angle, size)
+    x = np.sin(r) * np.cos(alpha)
+    y = np.sin(r) * np.sin(alpha)
+    z = np.cos(r)
+
+    r_vector = np.array([x, y, z])
+    return r_vector
+
+
+def _vector_to_angles(vector):
+    r"""
+    Generate inclination, declination, and intensity from a 3-component vector
+
+    .. note::
+
+        Inclination is positive downwards and declination is the angle with
+        the y component. The vector has x, y, and z (upward) Cartesian
+        components.
+
+    Parameters
+    ----------
+    vector : 1D or 2D array
+        The x, y, z vector components. Can be a 1D array for a single vector
+        or 2D for multiple. If 2D, then each vector should be a row of the
+        array.
+
+    Returns
+    -------
+    inclination : float or array
+        Inclination of the magnetic vector in degrees.
+    declination : float or array
+        Declination of the magnetic vector in degrees
+    intensity : float or array
+        Intensity of the magnetic vector.
+    """
+    vector = np.asarray(vector)
+    x, y, z = vector.T
+    intensity = np.sqrt(x**2 + y**2 + z**2)
+    inclination = -np.degrees(np.arctan2(z, np.hypot(x, y)))
+    declination = np.degrees(np.arctan2(x, y))
+    return inclination, declination, intensity

magali/tests/test_io.py

@@ -8,9 +8,11 @@
 Test the IO functions
 """
 
+import numpy as np
 import pooch
+import scipy as sp
 
-from .._input_output import read_qdm_harvard
+from .._input_output import random_unitary_vector, read_qdm_harvard
 
 
 def test_read_qdm_harvard():
@@ -21,3 +23,12 @@ def test_read_qdm_harvard():
     )
     data = read_qdm_harvard(path)
     assert data.bz.size == 576_000
+
+
+def test_random_unitary_vector():
+    "Perform a shapiro test, in order to check normal distribution"
+    r_vector = random_unitary_vector(np.deg2rad(50))
+
+    assert (sp.stats.shapiro(r_vector[0]).pvalue >= 0.05) and (
+        sp.stats.shapiro(r_vector[1]).pvalue >= 0.05
+    )