Merge pull request #485 from hongyuchen1030/GCA_ConstantLat_Intersection

Gca constant lat intersection

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
 =========================

docs/user_api/index.rst

@@ -273,6 +273,7 @@ Intersections
    :toctree: _autosummary
 
    grid.intersections.gca_gca_intersection
+   grid.intersections.gca_constLat_intersection
 
 Utils
 -----

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,23 @@ 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)])))

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,79 @@ 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.
+
+    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
+
+    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)
+    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
+    if point_within_gca(p1, gca_cart):
+        converged_pt = _newton_raphson_solver_for_gca_constLat(p1,
+                                                               gca_cart,
+                                                               verbose=verbose)
+    elif point_within_gca(p2, gca_cart):
+        converged_pt = _newton_raphson_solver_for_gca_constLat(p2,
+                                                               gca_cart,
+                                                               verbose=verbose)
+    else:
+        converged_pt = np.array([])
+    return converged_pt

uxarray/grid/utils.py

@@ -1,4 +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):
@@ -121,3 +123,109 @@ 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 _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)
+    #
+    # # 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:
+        raise warnings("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.
+
+    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)