pyfracval.geometry#

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#

pyfracval.geometry.norm3(v)[source]#

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).

pyfracval.geometry.cross3(a, b)[source]#

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.

pyfracval.geometry.rodrigues_rotation(vectors, axis, angle)[source]#

Rotate vector(s) around an axis using Rodrigues’ rotation formula.

Parameters:

  • vectors (np.ndarray) – A single 3D vector or an Nx3 array of vectors to rotate.

  • axis (np.ndarray) – The 3D rotation axis (does not need to be normalized).

  • angle (float) – The rotation angle in radians.

Returns:

The rotated vector or Nx3 array of rotated vectors. Returns the original vectors if the axis norm is near zero.

Return type:

np.ndarray

Raises:

ValueError – If input vectors is not 1D (3,) or 2D (N, 3).

pyfracval.geometry.two_sphere_intersection(sphere_1, sphere_2, rng=None)[source]#

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.

Parameters:

  • sphere_1 (np.ndarray) – Definition of the first sphere: [x1, y1, z1, r1].

  • sphere_2 (np.ndarray) – Definition of the second sphere: [x2, y2, z2, r2].

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).

Return type:

tuple[float, float, float, float, np.ndarray, np.ndarray, np.ndarray, bool]

Note

https://mathworld.wolfram.com/Sphere-SphereIntersection.html

pyfracval.geometry.spherical_cap_angle(sphere_1, sphere_2)[source]#

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).

Parameters:

  • sphere_1 (np.ndarray) – [x, y, z, r] for each sphere.

  • sphere_2 (np.ndarray) – [x, y, z, r] for each sphere.

Returns:

The cap half-angle in radians, measured from the axis connecting the two sphere centers.

Return type:

float

pyfracval.geometry.random_point_sc(case, spheres_1_ext, spheres_2_ext, rng=None)[source]#

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 two_sphere_intersection()) doesn’t capture - see the case 1/2/3 branch in cca/sticking.py::_cca_sticking_v1.

Parameters:

  • case (int) – Which shell-overlap case applies (1, 2, or 3 - see caller).

  • spheres_1_ext (np.ndarray) – [x, y, z, d_min, d_max] for each cluster’s shell.

  • spheres_2_ext (np.ndarray) – [x, y, z, d_min, d_max] for each cluster’s shell.

  • rng (np.random.Generator, optional) – Random generator to use; a fresh default one if not given.

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.

Return type:

tuple[float, float, float, bool]