pyfracval.geometry ================== .. py:module:: pyfracval.geometry .. autoapi-nested-parse:: Geometric primitives for PyFracVAL. Rodrigues rotation, sphere intersection, and related constants. Functions --------- rodrigues_rotation Rotate vectors around an axis using Rodrigues' formula. two_sphere_intersection Compute the intersection circle of two overlapping spheres. spherical_cap_angle Critical polar angle of a spherical cap (``ext_case=1`` support). random_point_sc Sample a random point on a spherical cap (``ext_case=1`` support). Constants --------- FLOATING_POINT_ERROR Numerical tolerance for floating-point comparisons. Module Contents --------------- .. py:function:: norm3(v) Euclidean norm of a single 3-vector. ``np.linalg.norm`` is a general N-dimensional routine: it validates axes, coerces inputs and dispatches, which costs about 1.3 us against 0.6 us here. That is irrelevant once per array and very relevant at the ~20k calls per aggregate the CCA sticking path makes on plain 3-vectors (measured with benchmarks/profile_pipeline.py). .. py:function:: cross3(a, b) Cross product of two 3-vectors. ``np.cross`` is dramatically worse than ``np.linalg.norm`` for this: it supports 2- and 3-component inputs over arbitrary axes with broadcasting, and pays roughly 21 us per call against 1.7 us for the three explicit components. It is called from the rotation setup on every sticking attempt. .. py:function:: rodrigues_rotation(vectors, axis, angle) Rotate vector(s) around an axis using Rodrigues' rotation formula. :param vectors: A single 3D vector or an Nx3 array of vectors to rotate. :type vectors: np.ndarray :param axis: The 3D rotation axis (does not need to be normalized). :type axis: np.ndarray :param angle: The rotation angle in radians. :type angle: float :returns: The rotated vector or Nx3 array of rotated vectors. Returns the original vectors if the axis norm is near zero. :rtype: np.ndarray :raises ValueError: If input `vectors` is not 1D (3,) or 2D (N, 3). .. py:function:: two_sphere_intersection(sphere_1, sphere_2, rng = None) Find the intersection circle of two spheres and pick a random point. Calculates the center (x0, y0, z0) and radius (r0) of the intersection circle, defines basis vectors (i_vec, j_vec) for the circle's plane, and returns a random point (x, y, z) on that circle based on a random angle (theta). Handles edge cases: spheres too far, one contained, coincidence, touching. :param sphere_1: Definition of the first sphere: [x1, y1, z1, r1]. :type sphere_1: np.ndarray :param sphere_2: Definition of the second sphere: [x2, y2, z2, r2]. :type sphere_2: np.ndarray :returns: A tuple containing: - x, y, z (float): Coordinates of a random point on the intersection. - theta (float): Random angle (radians) used to generate the point. - vec_0 (np.ndarray): [x0, y0, z0, r0] - center and radius of the intersection circle (r0=0 if spheres touch at a point). - i_vec (np.ndarray): First basis vector of the intersection plane. - j_vec (np.ndarray): Second basis vector of the intersection plane. - valid (bool): True if a valid intersection (circle or point) exists, False otherwise (e.g., separate, contained, coincident). :rtype: tuple[float, float, float, float, np.ndarray, np.ndarray, np.ndarray, bool] .. note:: https://mathworld.wolfram.com/Sphere-SphereIntersection.html .. py:function:: spherical_cap_angle(sphere_1, sphere_2) Critical polar angle of the spherical cap of ``sphere_1`` that lies within ``sphere_2`` (Fortran ``Spherical_cap_angle``, see ``docs/FracVAL/CCA_module.f90``). :param sphere_1: [x, y, z, r] for each sphere. :type sphere_1: np.ndarray :param sphere_2: [x, y, z, r] for each sphere. :type sphere_2: np.ndarray :returns: The cap half-angle in radians, measured from the axis connecting the two sphere centers. :rtype: float .. py:function:: random_point_sc(case, spheres_1_ext, spheres_2_ext, rng = None) Sample a random point on the appropriate spherical cap for CCA's ``ext_case=1`` contact-point search (Fortran ``Random_point_SC``, see ``docs/FracVAL/CCA_module.f90``). Used when the "shell" spheres (defined by [Dmin, Dmax] distance from each cluster's center of mass) overlap in a way that a single Dmax/Dmax intersection circle (the ``ext_case=0`` path used by :func:`two_sphere_intersection`) doesn't capture - see the ``case`` 1/2/3 branch in ``cca/sticking.py::_cca_sticking_v1``. :param case: Which shell-overlap case applies (1, 2, or 3 - see caller). :type case: int :param spheres_1_ext: [x, y, z, d_min, d_max] for each cluster's shell. :type spheres_1_ext: np.ndarray :param spheres_2_ext: [x, y, z, d_min, d_max] for each cluster's shell. :type spheres_2_ext: np.ndarray :param rng: Random generator to use; a fresh default one if not given. :type rng: np.random.Generator, optional :returns: (x, y, z, valid) - a random point on the selected sphere's surface, within the appropriate cap angle. ``valid`` is False for an unrecognized case or a degenerate (coincident) center pair. :rtype: tuple[float, float, float, bool]