From f69c195af8f0f4df342367898f5b15f8a2d9dee7 Mon Sep 17 00:00:00 2001
From: Hongyu Chen <hyvchen@ucdavis.edu>
Date: Wed, 27 Sep 2023 10:49:43 -0700
Subject: [PATCH 1/5] Initial implementation

---
 test/test_intersections.py    | 19 ++++++++++++-
 uxarray/grid/intersections.py | 53 ++++++++++++++++++++++++++++++++++-
 uxarray/grid/utils.py         | 45 +++++++++++++++++++++++++++++
 3 files changed, 115 insertions(+), 2 deletions(-)

diff --git a/test/test_intersections.py b/test/test_intersections.py
index 15c3da534..72a525e63 100644
--- a/test/test_intersections.py
+++ b/test/test_intersections.py
@@ -4,7 +4,7 @@
 from uxarray.constants import ERROR_TOLERANCE
 
 from uxarray.grid.coordinates import node_lonlat_rad_to_xyz, node_xyz_to_lonlat_rad
-from uxarray.grid.intersections import gca_gca_intersection
+from uxarray.grid.intersections import gca_gca_intersection, gca_constLat_intersection
 
 
 class TestGCAGCAIntersection(TestCase):
@@ -90,3 +90,20 @@ def test_get_GCA_GCA_intersections_perpendicular(self):
             np.allclose(res_lonlat_rad,
                         np.array([np.deg2rad(170.0),
                                   np.deg2rad(0.0)])))
+
+class TestGCAconstLatIntersection(TestCase):
+
+    def test_GCA_constLat_intersections_antimeridian(self):
+        GCR1_cart = np.array([
+            node_lonlat_rad_to_xyz([np.deg2rad(170.0),
+                                    np.deg2rad(89.99)]),
+            node_lonlat_rad_to_xyz([np.deg2rad(170.0),
+                                    np.deg2rad(10.0)])
+        ])
+
+        res = gca_constLat_intersection(GCR1_cart, np.deg2rad(60.0), verbose=True)
+        res_lonlat_rad = node_xyz_to_lonlat_rad(res.tolist())
+        self.assertTrue(
+            np.allclose(res_lonlat_rad,
+                        np.array([np.deg2rad(170.0),
+                                  np.deg2rad(60.0)])))
\ No newline at end of file
diff --git a/uxarray/grid/intersections.py b/uxarray/grid/intersections.py
index c3b8cece0..81c356fa6 100644
--- a/uxarray/grid/intersections.py
+++ b/uxarray/grid/intersections.py
@@ -1,6 +1,6 @@
 import numpy as np
 from uxarray.constants import ERROR_TOLERANCE
-from uxarray.grid.utils import cross_fma
+from uxarray.grid.utils import cross_fma, _newton_raphson_solver_for_gca_constLat
 from uxarray.grid.lines import point_within_gca
 import platform
 import warnings
@@ -115,3 +115,54 @@ def gca_gca_intersection(gca1_cart, gca2_cart, fma_disabled=False):
         res = np.append(res, x2)
 
     return res
+
+def gca_constLat_intersection(gca_cart, constLat, fma_disabled=False, verbose=False):
+    """Calculate the intersection point(s) of a Great Circle Arc (GCA) and a constant latitude line in a
+    Cartesian coordinate system.
+
+    To reduce relative errors, the Fused Multiply-Add (FMA) operation is utilized.
+    A warning is raised if the given coordinates are not in the cartesian coordinates, or
+    they cannot be accurately handled using floating-point arithmetic.
+
+    Parameters
+    ----------
+    gca_cart : [2, 3] np.ndarray Cartesian coordinates of the two end points GCA.
+    constLat : float
+        The constant latitude of the latitude line.
+    fma_disabled : bool, optional (default=False)
+        If True, the FMA operation is disabled. And a naive `np.cross` is used instead.
+
+    Returns
+    -------
+    np.ndarray
+        Cartesian coordinates of the intersection point(s).
+
+    Raises
+    ------
+    ValueError
+        If the input GCA is not in the cartesian [x, y, z] format.
+    """
+    constZ = np.sin(constLat)
+    x1, x2 = gca_cart
+    n = cross_fma(x1, x2)
+    nx, ny, nz = n
+
+    s_tilde = np.sqrt(nx ** 2 + ny ** 2 - np.linalg.norm(n) ** 2 * constZ ** 2)
+    p1_x = -(1.0 / (nx ** 2 + ny ** 2)) * (constZ * nx * nz + s_tilde * ny)
+    p2_x = -(1.0 / (nx ** 2 + ny ** 2)) * (constZ * nx * nz - s_tilde * ny)
+    p1_y = -(1.0 / (nx ** 2 + ny ** 2)) * (constZ * ny * nz - s_tilde * nx)
+    p2_y = -(1.0 / (nx ** 2 + ny ** 2)) * (constZ * ny * nz + s_tilde * nx)
+
+    p1 = np.array([p1_x, p1_y, constZ])
+    p2 = np.array([p2_x, p2_y, constZ])
+
+    # Now test which intersection point is within the GCA range
+    res = np.array([])
+    if point_within_gca(p1, gca_cart):
+        converged_pt = _newton_raphson_solver_for_gca_constLat(p1, gca_cart, verbose=verbose)
+        return converged_pt
+    elif point_within_gca(p2, gca_cart):
+        converged_pt = _newton_raphson_solver_for_gca_constLat(p2, gca_cart, verbose=verbose)
+        return converged_pt
+
+    return res
diff --git a/uxarray/grid/utils.py b/uxarray/grid/utils.py
index 6524399b9..c6bf53b43 100644
--- a/uxarray/grid/utils.py
+++ b/uxarray/grid/utils.py
@@ -121,3 +121,48 @@ def cross_fma(v1, v2):
     y = _fmms(v1[2], v2[0], v1[0], v2[2])
     z = _fmms(v1[0], v2[1], v1[1], v2[0])
     return np.array([x, y, z])
+
+def _newton_raphson_solver_for_gca_constLat(init_cart, gca_cart, max_iter=1000, verbose=False):
+    """
+    Solve for the intersection point between a great circle arc and a constant latitude.
+
+    Args:
+        init_cart (np.ndarray): Initial guess for the intersection point.
+        w0_cart (np.ndarray): First vector defining the great circle arc.
+        w1_cart (np.ndarray): Second vector defining the great circle arc.
+        max_iter (int, optional): Maximum number of iterations. Defaults to 1000.
+        verbose (bool, optional): Whether to print verbose output. Defaults to False.
+
+    Returns:
+        np.ndarray or None: The intersection point or None if the solver fails to converge.
+    """
+    tolerance = ERROR_TOLERANCE
+    w0_cart, w1_cart = gca_cart
+    error = float('inf')
+    constZ = init_cart[2]
+    y_guess = np.array(init_cart[0:2])
+    y_new = y_guess
+
+    _iter = 0
+
+    while error > tolerance and _iter < max_iter:
+        f_vector = np.array([
+            np.dot(np.cross(w0_cart, w1_cart), np.array([y_guess[0], y_guess[1], constZ])),
+            y_guess[0] * y_guess[0] + y_guess[1] * y_guess[1] + constZ * constZ - 1.0
+        ])
+
+        j_inv = __inv_jacobian(w0_cart[0], w1_cart[0], w0_cart[1], w1_cart[1], w0_cart[2], w1_cart[2], y_guess[0],
+                                 y_guess[1])
+
+        if j_inv is None:
+            return None
+
+        y_new = y_guess - np.matmul(j_inv, f_vector)
+        error = np.max(np.abs(y_guess - y_new))
+        y_guess = y_new
+
+        if verbose:
+            print(f"Newton method iter: {_iter}, error: {error}")
+        _iter += 1
+
+    return np.append(y_new, constZ)
\ No newline at end of file

From a0d7ad80864e0f57adec157f64d5885965c3f213 Mon Sep 17 00:00:00 2001
From: Hongyu Chen <hyvchen@ucdavis.edu>
Date: Wed, 27 Sep 2023 12:11:46 -0700
Subject: [PATCH 2/5] inverse jacobian

---
 uxarray/grid/utils.py | 24 +++++++++++++++++++++++-
 1 file changed, 23 insertions(+), 1 deletion(-)

diff --git a/uxarray/grid/utils.py b/uxarray/grid/utils.py
index c6bf53b43..02a99e42c 100644
--- a/uxarray/grid/utils.py
+++ b/uxarray/grid/utils.py
@@ -1,5 +1,5 @@
 import numpy as np
-
+from uxarray.constants import ERROR_TOLERANCE
 
 def _replace_fill_values(grid_var, original_fill, new_fill, new_dtype=None):
     """Replaces all instances of the current fill value (``original_fill``) in
@@ -122,6 +122,28 @@ def cross_fma(v1, v2):
     z = _fmms(v1[0], v2[1], v1[1], v2[0])
     return np.array([x, y, z])
 
+def __inv_jacobian(x0, x1, y0, y1, z0, z1, x_i_old, y_i_old):
+    # d_dx = (x0 * x_i_old - x1 * x_i_old * z0 + y0 * y_i_old * z1 - y1 * y_i_old * z0 - y1 * y_i_old * z0)
+    # d_dy = 2 * (x0 * x_i_old * z1 - x1 * x_i_old * z0 + y0 * y_i_old * z1 - y1 * y_i_old * z0)
+    #
+    # # row 1
+    # J[0, 0] = y_i_old / d_dx
+    # J[0, 1] = (x0 * z1 - z0 * x1) / d_dy
+    # # row 2
+    # J[1, 0] = x_i_old / d_dx
+    # J[1, 1] = (y0 * z1 - z0 * y1) / d_dy
+
+    # The Jacobian Matrix
+    jacobian = [[_fmms(y0, z1, z0, y1), _fmms(x0, z1, z0, x1)],
+                [2 * x_i_old, 2 * y_i_old]]
+    try:
+        inverse_jacobian = np.linalg.inv(jacobian)
+    except np.linalg.LinAlgError:
+        print("Warning: Singular Jacobian matrix encountered.")
+        return None
+
+    return inverse_jacobian
+
 def _newton_raphson_solver_for_gca_constLat(init_cart, gca_cart, max_iter=1000, verbose=False):
     """
     Solve for the intersection point between a great circle arc and a constant latitude.

From f0ab8df7cc491d8688e5d5b1977d648e9ac4d168 Mon Sep 17 00:00:00 2001
From: Hongyu Chen <hyvchen@ucdavis.edu>
Date: Tue, 3 Oct 2023 14:34:49 -0700
Subject: [PATCH 3/5] Add the required documentation

---
 docs/user_api/index.rst       |  2 ++
 test/test_intersections.py    |  7 +++++--
 uxarray/grid/intersections.py | 28 ++++++++++++++++++----------
 uxarray/grid/utils.py         | 28 ++++++++++++++++++----------
 4 files changed, 43 insertions(+), 22 deletions(-)

diff --git a/docs/user_api/index.rst b/docs/user_api/index.rst
index 2e2652454..13593038c 100644
--- a/docs/user_api/index.rst
+++ b/docs/user_api/index.rst
@@ -273,6 +273,7 @@ Intersections
    :toctree: _autosummary
 
    grid.intersections.gca_gca_intersection
+   grid.intersections.gca_constLat_intersection
 
 Utils
 -----
@@ -280,6 +281,7 @@ Utils
    :toctree: _autosummary
 
    grid.utils.cross_fma
+   grid.utils._newton_raphson_solver_for_gca_constLat
 
 Numba
 -----
diff --git a/test/test_intersections.py b/test/test_intersections.py
index 72a525e63..e6543cc97 100644
--- a/test/test_intersections.py
+++ b/test/test_intersections.py
@@ -91,6 +91,7 @@ def test_get_GCA_GCA_intersections_perpendicular(self):
                         np.array([np.deg2rad(170.0),
                                   np.deg2rad(0.0)])))
 
+
 class TestGCAconstLatIntersection(TestCase):
 
     def test_GCA_constLat_intersections_antimeridian(self):
@@ -101,9 +102,11 @@ def test_GCA_constLat_intersections_antimeridian(self):
                                     np.deg2rad(10.0)])
         ])
 
-        res = gca_constLat_intersection(GCR1_cart, np.deg2rad(60.0), verbose=True)
+        res = gca_constLat_intersection(GCR1_cart,
+                                        np.deg2rad(60.0),
+                                        verbose=True)
         res_lonlat_rad = node_xyz_to_lonlat_rad(res.tolist())
         self.assertTrue(
             np.allclose(res_lonlat_rad,
                         np.array([np.deg2rad(170.0),
-                                  np.deg2rad(60.0)])))
\ No newline at end of file
+                                  np.deg2rad(60.0)])))
diff --git a/uxarray/grid/intersections.py b/uxarray/grid/intersections.py
index 81c356fa6..961884576 100644
--- a/uxarray/grid/intersections.py
+++ b/uxarray/grid/intersections.py
@@ -116,9 +116,13 @@ def gca_gca_intersection(gca1_cart, gca2_cart, fma_disabled=False):
 
     return res
 
-def gca_constLat_intersection(gca_cart, constLat, fma_disabled=False, verbose=False):
-    """Calculate the intersection point(s) of a Great Circle Arc (GCA) and a constant latitude line in a
-    Cartesian coordinate system.
+
+def gca_constLat_intersection(gca_cart,
+                              constLat,
+                              fma_disabled=False,
+                              verbose=False):
+    """Calculate the intersection point(s) of a Great Circle Arc (GCA) and a
+    constant latitude line in a Cartesian coordinate system.
 
     To reduce relative errors, the Fused Multiply-Add (FMA) operation is utilized.
     A warning is raised if the given coordinates are not in the cartesian coordinates, or
@@ -147,11 +151,11 @@ def gca_constLat_intersection(gca_cart, constLat, fma_disabled=False, verbose=Fa
     n = cross_fma(x1, x2)
     nx, ny, nz = n
 
-    s_tilde = np.sqrt(nx ** 2 + ny ** 2 - np.linalg.norm(n) ** 2 * constZ ** 2)
-    p1_x = -(1.0 / (nx ** 2 + ny ** 2)) * (constZ * nx * nz + s_tilde * ny)
-    p2_x = -(1.0 / (nx ** 2 + ny ** 2)) * (constZ * nx * nz - s_tilde * ny)
-    p1_y = -(1.0 / (nx ** 2 + ny ** 2)) * (constZ * ny * nz - s_tilde * nx)
-    p2_y = -(1.0 / (nx ** 2 + ny ** 2)) * (constZ * ny * nz + s_tilde * nx)
+    s_tilde = np.sqrt(nx**2 + ny**2 - np.linalg.norm(n)**2 * constZ**2)
+    p1_x = -(1.0 / (nx**2 + ny**2)) * (constZ * nx * nz + s_tilde * ny)
+    p2_x = -(1.0 / (nx**2 + ny**2)) * (constZ * nx * nz - s_tilde * ny)
+    p1_y = -(1.0 / (nx**2 + ny**2)) * (constZ * ny * nz - s_tilde * nx)
+    p2_y = -(1.0 / (nx**2 + ny**2)) * (constZ * ny * nz + s_tilde * nx)
 
     p1 = np.array([p1_x, p1_y, constZ])
     p2 = np.array([p2_x, p2_y, constZ])
@@ -159,10 +163,14 @@ def gca_constLat_intersection(gca_cart, constLat, fma_disabled=False, verbose=Fa
     # Now test which intersection point is within the GCA range
     res = np.array([])
     if point_within_gca(p1, gca_cart):
-        converged_pt = _newton_raphson_solver_for_gca_constLat(p1, gca_cart, verbose=verbose)
+        converged_pt = _newton_raphson_solver_for_gca_constLat(p1,
+                                                               gca_cart,
+                                                               verbose=verbose)
         return converged_pt
     elif point_within_gca(p2, gca_cart):
-        converged_pt = _newton_raphson_solver_for_gca_constLat(p2, gca_cart, verbose=verbose)
+        converged_pt = _newton_raphson_solver_for_gca_constLat(p2,
+                                                               gca_cart,
+                                                               verbose=verbose)
         return converged_pt
 
     return res
diff --git a/uxarray/grid/utils.py b/uxarray/grid/utils.py
index 02a99e42c..6ae662cf1 100644
--- a/uxarray/grid/utils.py
+++ b/uxarray/grid/utils.py
@@ -1,6 +1,7 @@
 import numpy as np
 from uxarray.constants import ERROR_TOLERANCE
 
+
 def _replace_fill_values(grid_var, original_fill, new_fill, new_dtype=None):
     """Replaces all instances of the current fill value (``original_fill``) in
     (``grid_var``) with (``new_fill``) and converts to the dtype defined by
@@ -122,6 +123,7 @@ def cross_fma(v1, v2):
     z = _fmms(v1[0], v2[1], v1[1], v2[0])
     return np.array([x, y, z])
 
+
 def __inv_jacobian(x0, x1, y0, y1, z0, z1, x_i_old, y_i_old):
     # d_dx = (x0 * x_i_old - x1 * x_i_old * z0 + y0 * y_i_old * z1 - y1 * y_i_old * z0 - y1 * y_i_old * z0)
     # d_dy = 2 * (x0 * x_i_old * z1 - x1 * x_i_old * z0 + y0 * y_i_old * z1 - y1 * y_i_old * z0)
@@ -134,8 +136,8 @@ def __inv_jacobian(x0, x1, y0, y1, z0, z1, x_i_old, y_i_old):
     # J[1, 1] = (y0 * z1 - z0 * y1) / d_dy
 
     # The Jacobian Matrix
-    jacobian = [[_fmms(y0, z1, z0, y1), _fmms(x0, z1, z0, x1)],
-                [2 * x_i_old, 2 * y_i_old]]
+    jacobian = [[_fmms(y0, z1, z0, y1),
+                 _fmms(x0, z1, z0, x1)], [2 * x_i_old, 2 * y_i_old]]
     try:
         inverse_jacobian = np.linalg.inv(jacobian)
     except np.linalg.LinAlgError:
@@ -144,9 +146,13 @@ def __inv_jacobian(x0, x1, y0, y1, z0, z1, x_i_old, y_i_old):
 
     return inverse_jacobian
 
-def _newton_raphson_solver_for_gca_constLat(init_cart, gca_cart, max_iter=1000, verbose=False):
-    """
-    Solve for the intersection point between a great circle arc and a constant latitude.
+
+def _newton_raphson_solver_for_gca_constLat(init_cart,
+                                            gca_cart,
+                                            max_iter=1000,
+                                            verbose=False):
+    """Solve for the intersection point between a great circle arc and a
+    constant latitude.
 
     Args:
         init_cart (np.ndarray): Initial guess for the intersection point.
@@ -169,12 +175,14 @@ def _newton_raphson_solver_for_gca_constLat(init_cart, gca_cart, max_iter=1000,
 
     while error > tolerance and _iter < max_iter:
         f_vector = np.array([
-            np.dot(np.cross(w0_cart, w1_cart), np.array([y_guess[0], y_guess[1], constZ])),
-            y_guess[0] * y_guess[0] + y_guess[1] * y_guess[1] + constZ * constZ - 1.0
+            np.dot(np.cross(w0_cart, w1_cart),
+                   np.array([y_guess[0], y_guess[1],
+                             constZ])), y_guess[0] * y_guess[0] +
+            y_guess[1] * y_guess[1] + constZ * constZ - 1.0
         ])
 
-        j_inv = __inv_jacobian(w0_cart[0], w1_cart[0], w0_cart[1], w1_cart[1], w0_cart[2], w1_cart[2], y_guess[0],
-                                 y_guess[1])
+        j_inv = __inv_jacobian(w0_cart[0], w1_cart[0], w0_cart[1], w1_cart[1],
+                               w0_cart[2], w1_cart[2], y_guess[0], y_guess[1])
 
         if j_inv is None:
             return None
@@ -187,4 +195,4 @@ def _newton_raphson_solver_for_gca_constLat(init_cart, gca_cart, max_iter=1000,
             print(f"Newton method iter: {_iter}, error: {error}")
         _iter += 1
 
-    return np.append(y_new, constZ)
\ No newline at end of file
+    return np.append(y_new, constZ)

From 82ff86eaa2bde9eab008ff4d38d1b0e3fefc40db Mon Sep 17 00:00:00 2001
From: Hongyu Chen <hyvchen@ucdavis.edu>
Date: Sun, 8 Oct 2023 14:28:32 -0700
Subject: [PATCH 4/5] Address PR review

---
 docs/internal_api/index.rst   |  2 ++
 uxarray/grid/intersections.py | 29 ++++++++++++++++++++-----
 uxarray/grid/utils.py         | 41 +++++++++++++++++++++++++++++++----
 3 files changed, 62 insertions(+), 10 deletions(-)

diff --git a/docs/internal_api/index.rst b/docs/internal_api/index.rst
index 02d6f1515..40bb660e9 100644
--- a/docs/internal_api/index.rst
+++ b/docs/internal_api/index.rst
@@ -146,6 +146,8 @@ Utils
    :toctree: _autosummary
 
    grid.utils._fmms
+   grid.utils._newton_raphson_solver_for_gca_constLat
+   grid.utils._inv_jacobian
 
 Grid Parsing and Encoding
 =========================
diff --git a/uxarray/grid/intersections.py b/uxarray/grid/intersections.py
index 961884576..dfccbaf81 100644
--- a/uxarray/grid/intersections.py
+++ b/uxarray/grid/intersections.py
@@ -145,10 +145,29 @@ def gca_constLat_intersection(gca_cart,
     ------
     ValueError
         If the input GCA is not in the cartesian [x, y, z] format.
+
+    Warning
+    -------
+        If running on the Windows system with fma_disabled=False since the C/C++ implementation of FMA in MS Windows
+        is fundamentally broken. (bug report: https://bugs.python.org/msg312480)
     """
     constZ = np.sin(constLat)
     x1, x2 = gca_cart
-    n = cross_fma(x1, x2)
+
+    if fma_disabled:
+        n = np.cross(x1, x2)
+
+    else:
+        # Raise a warning for Windows users
+        if platform.system() == "Windows":
+
+            warnings.warn(
+                "The C/C++ implementation of FMA in MS Windows is reportedly broken. Use with care. (bug report: "
+                "https://bugs.python.org/msg312480)"
+                "The single rounding cannot be guaranteed, hence the relative error bound of 3u cannot be guaranteed."
+            )
+        n = cross_fma(x1, x2)
+
     nx, ny, nz = n
 
     s_tilde = np.sqrt(nx**2 + ny**2 - np.linalg.norm(n)**2 * constZ**2)
@@ -161,16 +180,14 @@ def gca_constLat_intersection(gca_cart,
     p2 = np.array([p2_x, p2_y, constZ])
 
     # Now test which intersection point is within the GCA range
-    res = np.array([])
     if point_within_gca(p1, gca_cart):
         converged_pt = _newton_raphson_solver_for_gca_constLat(p1,
                                                                gca_cart,
                                                                verbose=verbose)
-        return converged_pt
     elif point_within_gca(p2, gca_cart):
         converged_pt = _newton_raphson_solver_for_gca_constLat(p2,
                                                                gca_cart,
                                                                verbose=verbose)
-        return converged_pt
-
-    return res
+    else:
+        converged_pt = np.array([])
+    return converged_pt
diff --git a/uxarray/grid/utils.py b/uxarray/grid/utils.py
index 6ae662cf1..6bfeb38b0 100644
--- a/uxarray/grid/utils.py
+++ b/uxarray/grid/utils.py
@@ -1,5 +1,6 @@
 import numpy as np
 from uxarray.constants import ERROR_TOLERANCE
+import warnings
 
 
 def _replace_fill_values(grid_var, original_fill, new_fill, new_dtype=None):
@@ -124,7 +125,39 @@ def cross_fma(v1, v2):
     return np.array([x, y, z])
 
 
-def __inv_jacobian(x0, x1, y0, y1, z0, z1, x_i_old, y_i_old):
+def _inv_jacobian(x0, x1, y0, y1, z0, z1, x_i_old, y_i_old):
+    """Calculate the inverse Jacobian matrix for a given set of parameters.
+
+    Parameters
+    ----------
+    x0 : float
+        Description of x0.
+    x1 : float
+        Description of x1.
+    y0 : float
+        Description of y0.
+    y1 : float
+        Description of y1.
+    z0 : float
+        Description of z0.
+    z1 : float
+        Description of z1.
+    x_i_old : float
+        Description of x_i_old.
+    y_i_old : float
+        Description of y_i_old.
+
+    Returns
+    -------
+    numpy.ndarray or None
+        The inverse Jacobian matrix if it is non-singular, or None if a singular matrix is encountered.
+
+    Notes
+    -----
+    This function calculates the inverse Jacobian matrix based on the provided parameters. If the Jacobian matrix
+    is singular, a warning is printed, and None is returned.
+    """
+
     # d_dx = (x0 * x_i_old - x1 * x_i_old * z0 + y0 * y_i_old * z1 - y1 * y_i_old * z0 - y1 * y_i_old * z0)
     # d_dy = 2 * (x0 * x_i_old * z1 - x1 * x_i_old * z0 + y0 * y_i_old * z1 - y1 * y_i_old * z0)
     #
@@ -141,7 +174,7 @@ def __inv_jacobian(x0, x1, y0, y1, z0, z1, x_i_old, y_i_old):
     try:
         inverse_jacobian = np.linalg.inv(jacobian)
     except np.linalg.LinAlgError:
-        print("Warning: Singular Jacobian matrix encountered.")
+        raise warnings("Warning: Singular Jacobian matrix encountered.")
         return None
 
     return inverse_jacobian
@@ -181,8 +214,8 @@ def _newton_raphson_solver_for_gca_constLat(init_cart,
             y_guess[1] * y_guess[1] + constZ * constZ - 1.0
         ])
 
-        j_inv = __inv_jacobian(w0_cart[0], w1_cart[0], w0_cart[1], w1_cart[1],
-                               w0_cart[2], w1_cart[2], y_guess[0], y_guess[1])
+        j_inv = _inv_jacobian(w0_cart[0], w1_cart[0], w0_cart[1], w1_cart[1],
+                              w0_cart[2], w1_cart[2], y_guess[0], y_guess[1])
 
         if j_inv is None:
             return None

From 625152b3d49b874e1d66265340d6c5cabcb0c8e0 Mon Sep 17 00:00:00 2001
From: Hongyu Chen <hyvchen@ucdavis.edu>
Date: Mon, 9 Oct 2023 11:27:30 -0700
Subject: [PATCH 5/5] Remove newton method from users' API

---
 docs/user_api/index.rst | 1 -
 1 file changed, 1 deletion(-)

diff --git a/docs/user_api/index.rst b/docs/user_api/index.rst
index 13593038c..a564aa75b 100644
--- a/docs/user_api/index.rst
+++ b/docs/user_api/index.rst
@@ -281,7 +281,6 @@ Utils
    :toctree: _autosummary
 
    grid.utils.cross_fma
-   grid.utils._newton_raphson_solver_for_gca_constLat
 
 Numba
 -----