More generic x-point check (#4172)

* More general x-point lcfs check in finder

* Rename and remove repetition

* add array use

Author: geograham

bluemira/equilibria/find.py

@@ -22,7 +22,12 @@
 from bluemira.base.look_and_feel import bluemira_warn
 from bluemira.equilibria.constants import B_TOLERANCE, X_TOLERANCE
 from bluemira.equilibria.error import EquilibriaError
-from bluemira.geometry.coordinates import Coordinates, get_area_2d, in_polygon
+from bluemira.geometry.coordinates import (
+    Coordinates,
+    get_area_2d,
+    in_polygon,
+    join_intersect,
+)
 from bluemira.utilities.tools import floatify
 
 if TYPE_CHECKING:
@@ -819,13 +824,15 @@ def get_flux_loop(psi_norm):
     if o_points is None or x_points is None:
         o_points, x_points = find_OX_points(x, z, psi)
 
+    primary_xp, primary_op = x_points[0], o_points[0]
     delta = high - low
 
     while delta > psi_n_tol:
         middle = low + delta / 2
         flux_surface = get_flux_loop(middle)
 
-        if flux_surface.closed:
+        ixp = x_point_check(flux_surface, primary_op, primary_xp)
+        if flux_surface.closed and ixp is None:
             # Middle flux surface is still closed, shift search bounds
             low = middle
 
@@ -869,6 +876,86 @@ def get_flux_loop(psi_norm):
     return lcfs, separatrix
 
 
+def x_point_check(flux_surface: Coordinates, op: Opoint, xp: Xpoint):
+    """
+    Check if there are intersections of a flux surface with the
+    line tangent to the o-point-x-point vector at the x-point.
+
+    Parameters
+    ----------
+    flux_surface:
+        Flux surface of interest, e.g., candidate closed
+        flux surfaces when finding LCFS
+    op:
+        Primary o-point for an equilibrium
+    xp:
+        X-point of interest, usually primary
+
+    Returns
+    -------
+    arg_inters:
+        Intersection indices
+    """
+    length = np.hypot(np.max(flux_surface.x), np.max(np.abs(flux_surface.z)))
+    tanget_line = two_point_angled_line(op, xp, length)  # default theta is tangent
+    _, arg_inters = join_intersect(flux_surface, tanget_line, get_arg=True)
+    return arg_inters.sort()
+
+
+def two_point_angled_line(
+    centre_point: PsiPoint | Coordinates | np.ndarray,
+    edge_point: PsiPoint | Coordinates | np.ndarray,
+    length: float,
+    theta: float = np.pi / 2,
+):
+    """
+    Make a Coordinate object for a line of a given length that is at a given
+    angle (default is tangent) to a surface with a reference radial vector
+    specified by two points.
+
+    Parameters
+    ----------
+    centre_point:
+        Start point of reference vector
+    edge_point:
+        End point of reference vector, on the surface
+        where we will take the tangent.
+    length:
+        Length to make the tangent line Coordinate
+    theta:
+        CCW angle in radians
+
+    Returns
+    -------
+    :
+        Coordinates of the tangent line with length=length
+    """
+    cp = (
+        np.array([[centre_point.x], [0.0], [centre_point.z]])
+        if isinstance(centre_point, PsiPoint | Coordinates)
+        else centre_point
+    )
+    tp = (
+        np.array([[edge_point.x], [0.0], [edge_point.z]])
+        if isinstance(edge_point, PsiPoint | Coordinates)
+        else edge_point
+    )
+    a = cp - tp
+    a_hat = a / np.linalg.norm(a)
+
+    rot_matrix = np.array([
+        [np.cos(theta), 0, -np.sin(theta)],
+        [0, 0, 0],
+        [np.sin(theta), 0, np.cos(theta)],
+    ])  # ccw rotation about y-axis
+
+    b_hat = rot_matrix @ a_hat
+
+    p1 = tp - b_hat * length / 2
+    p2 = tp + b_hat * length / 2
+    return Coordinates(np.array([p1, p2]).T)
+
+
 def grid_2d_contour(x: np.ndarray, z: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
     """
     Grid a smooth contour and get the outline of the cells it encompasses.

eudemo/eudemo/ivc/divertor_silhouette.py

@@ -20,7 +20,10 @@
 from bluemira.base.parameter_frame import Parameter, ParameterFrame
 from bluemira.base.reactor import ComponentManager
 from bluemira.builders.divertor import DivertorBuilder
-from bluemira.equilibria.find import find_flux_surface_through_point
+from bluemira.equilibria.find import (
+    find_flux_surface_through_point,
+    two_point_angled_line,
+)
 from bluemira.equilibria.find_legs import LegFlux
 from bluemira.geometry.tools import (
     interpolate_bspline,
@@ -215,33 +218,22 @@ def _make_angled_target(self, leg: LegPosition):
         # i.e. towards the increasing leg length
         post_target_point = sol.value_at(distance=self.leg_length[leg].value + 0.1)
 
-        a = post_target_point - target_point
-        a_hat = a / np.linalg.norm(a)
-
         # ccw angle
         theta = (
             np.deg2rad(self.params.div_targ_angle_ib.value)
             if leg is LegPosition.INNER
             else np.deg2rad(self.params.div_targ_angle_ob.value)
         )
 
-        rot_matrix = np.array([
-            [np.cos(theta), 0, -np.sin(theta)],
-            [0, 0, 0],
-            [np.sin(theta), 0, np.cos(theta)],
-        ])  # ccw rotation about y-axis
-
-        b_hat = rot_matrix @ a_hat
-
-        target_half_length = (
-            0.5 * self.params.div_Ltarg_ib.value
+        target_length = (
+            self.params.div_Ltarg_ib.value
             if leg is LegPosition.INNER
-            else 0.5 * self.params.div_Ltarg_ob.value
+            else self.params.div_Ltarg_ob.value
         )
 
-        p1 = target_point - b_hat * target_half_length
-        p2 = target_point + b_hat * target_half_length
-        return np.array([p1, p2]).T
+        return two_point_angled_line(
+            post_target_point, target_point, target_length, theta
+        )
 
     def make_dome(self, start: np.ndarray, end: np.ndarray, label: str) -> BluemiraWire:
         """