From 4729b69d7bb0d229ce538969c71726965af0b175 Mon Sep 17 00:00:00 2001
From: Souza-Junior <gelsonfersouza1@gmail.com>
Date: Fri, 6 Sep 2024 16:25:50 -0300
Subject: [PATCH 1/3] Add functions to convert angles to vectors and vice versa

These functions differ slightly from Harmonica's due to the convention in
paleomagnetism, where a positive vector is defined as downward.
---
 magali/_utils.py | 74 ++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 74 insertions(+)
 create mode 100644 magali/_utils.py

diff --git a/magali/_utils.py b/magali/_utils.py
new file mode 100644
index 00000000..0b0f622f
--- /dev/null
+++ b/magali/_utils.py
@@ -0,0 +1,74 @@
+# Copyright (c) 2024 The Magali Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+#
+# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
+#
+
+
+import numpy as np
+
+def vector_to_angles(vector):
+    """
+    Generate inclination, declination, and amplitude from a 3-component vector
+
+    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
+        The inclination values in degrees.
+    declination : float or array
+        The declination values in degrees.
+    amplitude : float or array
+        The vector amplitude values.
+    """
+    vector = np.asarray(vector)
+    x, y, z = vector.T
+    amplitude = 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, amplitude
+
+def angles_to_vector(inclination, declination, amplitude):
+    """
+    Generate a 3-component vector from inclination, declination, and amplitude
+
+    Inclination is positive downwards and declination is the angle with the y
+    component. The vector has x, y, and z (upward) Cartesian components.
+
+    Parameters
+    ----------
+    inclination : float or array
+        The inclination values in degrees.
+    declination : float or array
+        The declination values in degrees.
+    amplitude : float or array
+        The vector amplitude values.
+
+    Returns
+    -------
+    vector : 1D or 2D array
+        The calculated x, y, z vector components. 1D if it's a single vector.
+        If N vectors are calculated, the "vector" will have shape (N, 3) with
+        each vector in a row of the array.
+    """
+    inclination = np.radians(inclination)
+    declination = np.radians(declination)
+    amplitude = np.asarray(amplitude)
+    sin_inc = np.sin(-inclination)
+    cos_inc = np.cos(-inclination)
+    sin_dec = np.sin(declination)
+    cos_dec = np.cos(declination)
+    x = cos_inc * sin_dec * amplitude
+    y = cos_inc * cos_dec * amplitude
+    z = sin_inc * amplitude
+    return np.transpose([x, y, z])
\ No newline at end of file

From b5b2283f872fef9337405a31931f3d252e417e71 Mon Sep 17 00:00:00 2001
From: Souza-Junior <gelsonfersouza1@gmail.com>
Date: Fri, 6 Sep 2024 18:00:03 -0300
Subject: [PATCH 2/3] Add function for generating stereographic projection
 plots

For now the function only plots lines
---
 magali/_stereograms.py | 176 +++++++++++++++++++++++++++++++++++++++++
 1 file changed, 176 insertions(+)
 create mode 100644 magali/_stereograms.py

diff --git a/magali/_stereograms.py b/magali/_stereograms.py
new file mode 100644
index 00000000..3723ada4
--- /dev/null
+++ b/magali/_stereograms.py
@@ -0,0 +1,176 @@
+# Copyright (c) 2024 The Magali Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+#
+# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
+#
+
+
+import numpy as np
+import matplotlib.pyplot as plt
+from ._utils import vector_to_angles, angles_to_vector
+
+def equal_area_projection(vectors):
+
+    """
+    Perform an equal-area projection of 3D vectors onto a 2D plane.
+
+    Parameters
+    ----------
+    vectors : ndarray of shape (n, 3)
+        Array of 3D vectors to be projected. Each row represents a vector 
+        with components [x, y, z].
+
+    Returns
+    -------
+    XY_projected : ndarray of shape (n, 3)
+        Array containing the 2D projected coordinates each vector.
+        The first two columns represent the x and y coordinates on the 
+        2D plane. The third column contains the amplitude, where the 
+        sign indicates whether the inclination is positive (downward) 
+        or negative (upward).
+    """
+
+    norm = np.linalg.norm(vectors, axis=1)
+    vectors_unitary = vectors / norm[:, np.newaxis]
+    inclinations, declinations, amplitudes = vector_to_angles(vectors)
+
+    XY_projected = np.zeros((len(vectors), 3))
+    for i, projected_vector in enumerate(vectors_unitary):
+        r = np.sqrt(1 - np.abs(vectors_unitary[i, 2])) / np.sqrt(vectors_unitary[i, 0]**2 + vectors_unitary[i, 1]**2)
+        XY_projected[i, 0] = r * vectors_unitary[i, 1]
+        XY_projected[i, 1] = r * vectors_unitary[i, 0]
+        XY_projected[i, 2] = amplitudes[i] if inclinations[i] >= 0 else -amplitudes[i]
+
+    return XY_projected
+
+
+
+class StereographicProjection:
+    def __init__(self, vectors):
+        self.vectors = vectors
+
+    def plot(self, ax=None, cmap="inferno", cmap_norm=plt.Normalize, vmin=None, vmax=None, 
+             label="", s=50, add_ticks=True, draw_cross=True, add_radial_grid=True, 
+             facecolor="#ffffff00", add_legend=False,
+             **kwargs):
+    """
+    Plot a stereographic projection with a radius of 1.
+
+    Parameters
+    ----------
+    ax : matplotlib.axes.Axes, optional
+        The axis on which to plot. If None, a new figure and axis will be created.
+    cmap : str, optional
+        Colormap to be used for the plotted data. Default is 'inferno'.
+    cmap_norm : matplotlib.colors.Normalize, optional
+        Normalization function for the colormap. Default is plt.Normalize.
+    vmin : float, optional
+        Minimum data value for colormap normalization. If None, the minimum value of the data is used.
+    vmax : float, optional
+        Maximum data value for colormap normalization. If None, the maximum value of the data is used.
+    label : str, optional
+        Label for the plotted data. Default is an empty string.
+    s : float, optional
+        Marker size for the scatter plot. Default is 50.
+    draw_cross : bool, optional
+        Whether to draw the central cross on the plot. Default is True.
+    add_radial_grid : bool, optional
+        Whether to add a radial grid to the plot. Default is True.
+    facecolor : str, optional
+        The facecolor of the plot background. Default is transparent ("#ffffff00").
+    add_legend : bool, optional
+        The legend of the plot. Default is False.
+    **kwargs : dict, optional
+        Additional keyword arguments passed to the scatter function.
+
+    Returns
+    -------
+    mappable : matplotlib.cm.ScalarMappable
+        The ScalarMappable object, which can be used to create a colorbar.
+    """
+    
+        if ax is None:
+            fig = plt.figure()
+            ax = fig.add_subplot(1, 1, 1)
+
+        # Add a face color
+        background_circle = plt.Circle((0, 0), 1, color=facecolor, zorder=-2)
+        ax.add_artist(background_circle)
+
+        # Draw the great circle
+        circle = plt.Circle((0, 0), 1, color='black', fill=False, zorder=3)
+        ax.add_artist(circle)
+
+        # Define the clipping area
+        clip_path = plt.Circle((0, 0), 1, transform=ax.transData)
+
+        if add_ticks:
+            # Add ticks 
+            ax.text(-0.025, 1.05, '0°')
+            ax.text(-0.05, -1.075, '180°')
+            ax.text(-1.15, -0.025, '270°')
+            ax.text(1.025, -0.025, '90°')
+        
+        if label and not label.endswith(" "):
+            label = label + " "
+            # Draw the central cross if requested
+            if draw_cross:
+                hline = ax.axhline(y=0, color='black', zorder=-1)
+                hline.set_clip_path(clip_path)
+                vline = ax.axvline(x=0, color='black', zorder=-1)
+                vline.set_clip_path(clip_path)
+
+
+        # Draw the radial grid
+        if add_radial_grid:
+            for rad_dec in range(0, 360, 10):  # Grid lines every 10 degrees
+                rad_inc = np.linspace(0, 90, 1000)  # Inclinations from 0 to 90 degrees
+                # Generate the radial unitary vectors
+                radial_vector = angles_to_vector(rad_inc, rad_dec, 1)
+                # Project the radial vectors
+                radial_projected = equal_area_projection(radial_vector)
+
+                # Plot each radial grid line
+                ax.plot(radial_projected[:, 1], radial_projected[:, 0], color='gray', zorder=-2, lw=0.5)
+            for circ_inc in range(0, 90, 10): # Isoinclination lines every 10 degrees
+                circ_dec = np.linspace(0, 360, 1000)  # Declinations from 0 to 360 degrees
+                # Generate the circular unitary vectors
+                circle_vector = angles_to_vector(circ_inc, circ_dec, 1)
+                # Project the radial vectors
+                circle_projected = equal_area_projection(circle_vector)
+
+                # Plot each radial grid line
+                ax.plot(circle_projected[:, 1], circle_projected[:, 0], color='gray', zorder=-2, lw=0.5)
+        
+        # Calculate the equal area projection
+        XY_projected = equal_area_projection(self.vectors)
+
+        # Generate colors based on the amplitude values
+        norm = cmap_norm(vmin=vmin, vmax=vmax)
+        colors = plt.colormaps[cmap](norm(abs(XY_projected[:, 2])))
+
+        # Plotting the data
+        positive_inc = XY_projected[:, 2] > 0
+        scatter_pos = ax.scatter(XY_projected[:, 1][positive_inc], XY_projected[:, 0][positive_inc], 
+                                 c=colors[positive_inc], edgecolors="#333333", label=f"{label}$I > 0$",
+                                **kwargs)
+        scatter_pos.set_clip_path(clip_path)
+
+        scatter_neg = ax.scatter(XY_projected[:, 1][~positive_inc], XY_projected[:, 0][~positive_inc],
+                                 c="#ffffff00", edgecolors=colors[~positive_inc], label=f"{label}$I \leq 0$",
+                                **kwargs)
+        scatter_neg.set_clip_path(clip_path)
+
+        # Configure the axis
+        ax.set_aspect('equal')
+        ax.set_xticks([])
+        ax.set_yticks([])
+        if add_legend:
+            ax.legend()
+        for spine in ax.spines.values():
+            spine.set_visible(False)
+
+        # To generate colorbar if necessary
+        mappable = plt.cm.ScalarMappable(norm=norm, cmap=cmap)
+        return mappable
\ No newline at end of file

From 168220050b330314240b834f21452434f1585678 Mon Sep 17 00:00:00 2001
From: Souza-Junior <gelsonfersouza1@gmail.com>
Date: Fri, 6 Sep 2024 18:51:12 -0300
Subject: [PATCH 3/3] Fix formatting issues

---
 magali/_stereograms.py | 200 ++++++++++++++++++++++++-----------------
 magali/_utils.py       |   4 +-
 2 files changed, 122 insertions(+), 82 deletions(-)

diff --git a/magali/_stereograms.py b/magali/_stereograms.py
index 3723ada4..81d1acef 100644
--- a/magali/_stereograms.py
+++ b/magali/_stereograms.py
@@ -6,28 +6,29 @@
 #
 
 
-import numpy as np
 import matplotlib.pyplot as plt
-from ._utils import vector_to_angles, angles_to_vector
+import numpy as np
+
+from ._utils import angles_to_vector, vector_to_angles
 
-def equal_area_projection(vectors):
 
+def equal_area_projection(vectors):
     """
     Perform an equal-area projection of 3D vectors onto a 2D plane.
 
     Parameters
     ----------
     vectors : ndarray of shape (n, 3)
-        Array of 3D vectors to be projected. Each row represents a vector 
+        Array of 3D vectors to be projected. Each row represents a vector
         with components [x, y, z].
 
     Returns
     -------
-    XY_projected : ndarray of shape (n, 3)
+    xy_projected : ndarray of shape (n, 3)
         Array containing the 2D projected coordinates each vector.
-        The first two columns represent the x and y coordinates on the 
-        2D plane. The third column contains the amplitude, where the 
-        sign indicates whether the inclination is positive (downward) 
+        The first two columns represent the x and y coordinates on the
+        2D plane. The third column contains the amplitude, where the
+        sign indicates whether the inclination is positive (downward)
         or negative (upward).
     """
 
@@ -35,61 +36,75 @@ def equal_area_projection(vectors):
     vectors_unitary = vectors / norm[:, np.newaxis]
     inclinations, declinations, amplitudes = vector_to_angles(vectors)
 
-    XY_projected = np.zeros((len(vectors), 3))
-    for i, projected_vector in enumerate(vectors_unitary):
-        r = np.sqrt(1 - np.abs(vectors_unitary[i, 2])) / np.sqrt(vectors_unitary[i, 0]**2 + vectors_unitary[i, 1]**2)
-        XY_projected[i, 0] = r * vectors_unitary[i, 1]
-        XY_projected[i, 1] = r * vectors_unitary[i, 0]
-        XY_projected[i, 2] = amplitudes[i] if inclinations[i] >= 0 else -amplitudes[i]
-
-    return XY_projected
+    xy_projected = np.zeros((len(vectors), 3))
+    for i, _ in enumerate(vectors_unitary):
+        r = np.sqrt(1 - np.abs(vectors_unitary[i, 2])) / np.sqrt(
+            vectors_unitary[i, 0] ** 2 + vectors_unitary[i, 1] ** 2
+        )
+        xy_projected[i, 0] = r * vectors_unitary[i, 1]
+        xy_projected[i, 1] = r * vectors_unitary[i, 0]
+        xy_projected[i, 2] = amplitudes[i] if inclinations[i] >= 0 else -amplitudes[i]
 
+    return xy_projected
 
 
 class StereographicProjection:
     def __init__(self, vectors):
         self.vectors = vectors
 
-    def plot(self, ax=None, cmap="inferno", cmap_norm=plt.Normalize, vmin=None, vmax=None, 
-             label="", s=50, add_ticks=True, draw_cross=True, add_radial_grid=True, 
-             facecolor="#ffffff00", add_legend=False,
-             **kwargs):
-    """
-    Plot a stereographic projection with a radius of 1.
-
-    Parameters
-    ----------
-    ax : matplotlib.axes.Axes, optional
-        The axis on which to plot. If None, a new figure and axis will be created.
-    cmap : str, optional
-        Colormap to be used for the plotted data. Default is 'inferno'.
-    cmap_norm : matplotlib.colors.Normalize, optional
-        Normalization function for the colormap. Default is plt.Normalize.
-    vmin : float, optional
-        Minimum data value for colormap normalization. If None, the minimum value of the data is used.
-    vmax : float, optional
-        Maximum data value for colormap normalization. If None, the maximum value of the data is used.
-    label : str, optional
-        Label for the plotted data. Default is an empty string.
-    s : float, optional
-        Marker size for the scatter plot. Default is 50.
-    draw_cross : bool, optional
-        Whether to draw the central cross on the plot. Default is True.
-    add_radial_grid : bool, optional
-        Whether to add a radial grid to the plot. Default is True.
-    facecolor : str, optional
-        The facecolor of the plot background. Default is transparent ("#ffffff00").
-    add_legend : bool, optional
-        The legend of the plot. Default is False.
-    **kwargs : dict, optional
-        Additional keyword arguments passed to the scatter function.
+    def plot(
+        self,
+        ax=None,
+        cmap="inferno",
+        cmap_norm=plt.Normalize,
+        vmin=None,
+        vmax=None,
+        label="",
+        add_ticks=True,
+        draw_cross=True,
+        add_radial_grid=True,
+        facecolor="#ffffff00",
+        add_legend=False,
+        **kwargs,
+    ):
+        r"""
+        Plot a stereographic projection with a radius of 1.
+
+        Parameters
+        ----------
+        ax : matplotlib.axes.Axes, optional
+            The axis on which to plot. If None, a new figure and axis will
+            be created.
+        cmap : str, optional
+            Colormap to be used for the plotted data. Default is 'inferno'.
+        cmap_norm : matplotlib.colors.Normalize, optional
+            Normalization function for the colormap. Default is plt.Normalize.
+        vmin : float, optional
+            Minimum data value for colormap normalization. If None,
+            the minimum value of the data is used.
+        vmax : float, optional
+            Maximum data value for colormap normalization. If None,
+            the maximum value of the data is used.
+        label : str, optional
+            Label for the plotted data. Default is an empty string.
+        draw_cross : bool, optional
+            Whether to draw the central cross on the plot. Default is True.
+        add_radial_grid : bool, optional
+            Whether to add a radial grid to the plot. Default is True.
+        facecolor : str, optional
+            The facecolor of the plot background. Default is
+            transparent ("#ffffff00").
+        add_legend : bool, optional
+            The legend of the plot. Default is False.
+        kwargs : dict, optional
+            Additional keyword arguments passed to the scatter function.
+
+        Returns
+        -------
+        mappable : matplotlib.cm.ScalarMappable
+            The ScalarMappable object, which can be used to create a colorbar.
+        """
 
-    Returns
-    -------
-    mappable : matplotlib.cm.ScalarMappable
-        The ScalarMappable object, which can be used to create a colorbar.
-    """
-    
         if ax is None:
             fig = plt.figure()
             ax = fig.add_subplot(1, 1, 1)
@@ -99,29 +114,28 @@ def plot(self, ax=None, cmap="inferno", cmap_norm=plt.Normalize, vmin=None, vmax
         ax.add_artist(background_circle)
 
         # Draw the great circle
-        circle = plt.Circle((0, 0), 1, color='black', fill=False, zorder=3)
+        circle = plt.Circle((0, 0), 1, color="black", fill=False, zorder=3)
         ax.add_artist(circle)
 
         # Define the clipping area
         clip_path = plt.Circle((0, 0), 1, transform=ax.transData)
 
         if add_ticks:
-            # Add ticks 
-            ax.text(-0.025, 1.05, '0°')
-            ax.text(-0.05, -1.075, '180°')
-            ax.text(-1.15, -0.025, '270°')
-            ax.text(1.025, -0.025, '90°')
-        
+            # Add ticks
+            ax.text(-0.025, 1.05, "0°")
+            ax.text(-0.05, -1.075, "180°")
+            ax.text(-1.15, -0.025, "270°")
+            ax.text(1.025, -0.025, "90°")
+
         if label and not label.endswith(" "):
             label = label + " "
             # Draw the central cross if requested
             if draw_cross:
-                hline = ax.axhline(y=0, color='black', zorder=-1)
+                hline = ax.axhline(y=0, color="black", zorder=-1)
                 hline.set_clip_path(clip_path)
-                vline = ax.axvline(x=0, color='black', zorder=-1)
+                vline = ax.axvline(x=0, color="black", zorder=-1)
                 vline.set_clip_path(clip_path)
 
-
         # Draw the radial grid
         if add_radial_grid:
             for rad_dec in range(0, 360, 10):  # Grid lines every 10 degrees
@@ -132,38 +146,62 @@ def plot(self, ax=None, cmap="inferno", cmap_norm=plt.Normalize, vmin=None, vmax
                 radial_projected = equal_area_projection(radial_vector)
 
                 # Plot each radial grid line
-                ax.plot(radial_projected[:, 1], radial_projected[:, 0], color='gray', zorder=-2, lw=0.5)
-            for circ_inc in range(0, 90, 10): # Isoinclination lines every 10 degrees
-                circ_dec = np.linspace(0, 360, 1000)  # Declinations from 0 to 360 degrees
+                ax.plot(
+                    radial_projected[:, 1],
+                    radial_projected[:, 0],
+                    color="gray",
+                    zorder=-2,
+                    lw=0.5,
+                )
+            for circ_inc in range(0, 90, 10):  # Isoinclination lines every 10 degrees
+                circ_dec = np.linspace(
+                    0, 360, 1000
+                )  # Declinations from 0 to 360 degrees
                 # Generate the circular unitary vectors
                 circle_vector = angles_to_vector(circ_inc, circ_dec, 1)
                 # Project the radial vectors
                 circle_projected = equal_area_projection(circle_vector)
 
                 # Plot each radial grid line
-                ax.plot(circle_projected[:, 1], circle_projected[:, 0], color='gray', zorder=-2, lw=0.5)
-        
+                ax.plot(
+                    circle_projected[:, 1],
+                    circle_projected[:, 0],
+                    color="gray",
+                    zorder=-2,
+                    lw=0.5,
+                )
+
         # Calculate the equal area projection
-        XY_projected = equal_area_projection(self.vectors)
+        xy_projected = equal_area_projection(self.vectors)
 
         # Generate colors based on the amplitude values
         norm = cmap_norm(vmin=vmin, vmax=vmax)
-        colors = plt.colormaps[cmap](norm(abs(XY_projected[:, 2])))
+        colors = plt.colormaps[cmap](norm(abs(xy_projected[:, 2])))
 
         # Plotting the data
-        positive_inc = XY_projected[:, 2] > 0
-        scatter_pos = ax.scatter(XY_projected[:, 1][positive_inc], XY_projected[:, 0][positive_inc], 
-                                 c=colors[positive_inc], edgecolors="#333333", label=f"{label}$I > 0$",
-                                **kwargs)
+        positive_inc = xy_projected[:, 2] > 0
+        scatter_pos = ax.scatter(
+            xy_projected[:, 1][positive_inc],
+            xy_projected[:, 0][positive_inc],
+            c=colors[positive_inc],
+            edgecolors="#333333",
+            label=f"{label}$I > 0$",
+            **kwargs,
+        )
         scatter_pos.set_clip_path(clip_path)
 
-        scatter_neg = ax.scatter(XY_projected[:, 1][~positive_inc], XY_projected[:, 0][~positive_inc],
-                                 c="#ffffff00", edgecolors=colors[~positive_inc], label=f"{label}$I \leq 0$",
-                                **kwargs)
+        scatter_neg = ax.scatter(
+            xy_projected[:, 1][~positive_inc],
+            xy_projected[:, 0][~positive_inc],
+            c="#ffffff00",
+            edgecolors=colors[~positive_inc],
+            label=rf"{label}$I \leq 0$",
+            **kwargs,
+        )
         scatter_neg.set_clip_path(clip_path)
 
         # Configure the axis
-        ax.set_aspect('equal')
+        ax.set_aspect("equal")
         ax.set_xticks([])
         ax.set_yticks([])
         if add_legend:
@@ -173,4 +211,4 @@ def plot(self, ax=None, cmap="inferno", cmap_norm=plt.Normalize, vmin=None, vmax
 
         # To generate colorbar if necessary
         mappable = plt.cm.ScalarMappable(norm=norm, cmap=cmap)
-        return mappable
\ No newline at end of file
+        return mappable
diff --git a/magali/_utils.py b/magali/_utils.py
index 0b0f622f..74a7a0d1 100644
--- a/magali/_utils.py
+++ b/magali/_utils.py
@@ -8,6 +8,7 @@
 
 import numpy as np
 
+
 def vector_to_angles(vector):
     """
     Generate inclination, declination, and amplitude from a 3-component vector
@@ -38,6 +39,7 @@ def vector_to_angles(vector):
     declination = np.degrees(np.arctan2(x, y))
     return inclination, declination, amplitude
 
+
 def angles_to_vector(inclination, declination, amplitude):
     """
     Generate a 3-component vector from inclination, declination, and amplitude
@@ -71,4 +73,4 @@ def angles_to_vector(inclination, declination, amplitude):
     x = cos_inc * sin_dec * amplitude
     y = cos_inc * cos_dec * amplitude
     z = sin_inc * amplitude
-    return np.transpose([x, y, z])
\ No newline at end of file
+    return np.transpose([x, y, z])