use {cgmath};

use {collision, component, constraint, math, object};

use geometry::shape;

/// Relation configuration of a pair of bounded objects.
///
/// - 'Disjoint' corresponds to a positive distance `> CONTACT_DISTANCE`
/// - 'Contact' corresponds to a positive distance `< CONTACT_DISTANCE`
/// - 'Inersect' corresponds to a strictly negative contact distance `< 0.0`
#[derive(Clone,Copy,Debug,Eq,PartialEq)]
pub enum Relation {
  Disjoint,
  Contact,
  Intersect
}

/// Result of a proximity query between a pair of bounded objects.
///
/// Note that the *axis vector* points from object A to object B, while the
/// *normal vector* points from object B to object A, but they should be
/// parallel when the axis vector is non-zero.
#[derive(Clone,Debug,PartialEq)]
pub struct Proximity {
  /// Signed distance.
  ///
  /// A positive distance indicates the separating distance for a disjoint
  /// result, and a negative distance indicates a penetration depth of an
  /// intersection result.
  pub distance  : f64,
  /// Non-normalized vector pointing along the axis of separation or penetration
  /// from object A to object B.
  ///
  /// For a disjoint configuration, adding this vector to the midpoint gives the
  /// nearest point on object B and subtracting this vector from the midpoint
  /// gives the nearest point on object A. Translating object A by this vector
  /// and translating objet B by its inverse will bring the objects into
  /// contact.
  ///
  /// For an intersecting configuration, translating object A by this vector and
  /// translating object B by its inverse will resolve the inter-penetration.
  ///
  /// Note this vector can be zero if the objects are in contact.
  pub half_axis : cgmath::Vector3 <f64>,
  /// Indicates the midpoint of the proximity query
  pub midpoint  : cgmath::Point3  <f64>,
  /// The unit normal to the separating plane directed from object B to object
  /// A
  pub normal    : math::Unit3 <f64>
}

impl Proximity {
  /// Proximity query.
  ///
  /// # Examples
  ///
  /// **Capsule v. capsule**
  ///
  /// By convention if the axes of capsules A and B intersect, the +X axis will
  /// be chosen for the axis vector and normal:
  ///
  /// ```
  /// # extern crate linear_sim;
  /// # use linear_sim::*;
  /// # use collision::*;
  /// # use geometry::shape;
  /// # use math;
  /// # fn main() {
  /// let a = object::Static {
  ///   position: component::Position ([0.0, 0.0, 1.0].into()),
  ///   bound:    component::Bound (shape::Bounded::from (
  ///     shape::Capsule::noisy (2.0, 3.0)).into()),
  ///   material: component::MATERIAL_STONE
  /// };
  /// let b = object::Static {
  ///   position: component::Position ([0.0, 0.0, -1.0].into()),
  ///   bound:    component::Bound (shape::Bounded::from (
  ///     shape::Capsule::noisy (1.0, 2.0)).into()),
  ///   material: component::MATERIAL_STONE
  /// };
  /// assert_eq!(
  ///   Proximity::query (&a, &b),
  ///   Proximity {
  ///     distance:  -3.0,
  ///     half_axis: [ 1.5, 0.0,  0.0].into(),
  ///     normal:    math::Unit3::axis_x(),
  ///     midpoint:  [-0.5, 0.0, -0.5].into()
  ///   }
  /// );
  /// # }
  /// ```
  pub fn query <A, B> (object_a : &A, object_b : &B) -> Self where
    A : object::Bounded, B : object::Bounded
  {
    use geometry::shape;
    let component::Position (ref position_a) = object_a.position();
    let component::Position (ref position_b) = object_b.position();
    let component::Bound    (ref shape_a)    = object_a.bound();
    let component::Bound    (ref shape_b)    = object_b.bound();

    match (shape_a, shape_b) {
      //
      // capsule vs. capsule
      //
      ( shape::Variant::Bounded (shape::Bounded::Capsule (ref capsule_a)),
        shape::Variant::Bounded (shape::Bounded::Capsule (ref capsule_b))
      ) => Self::query_capsule_capsule (
        position_a, capsule_a, position_b, capsule_b),
      //
      // cuboid vs. cuboid
      //
      ( shape::Variant::Bounded (shape::Bounded::Cuboid (ref cuboid_a)),
        shape::Variant::Bounded (shape::Bounded::Cuboid (ref cuboid_b))
      ) => Self::query_cuboid_cuboid (
        position_a, cuboid_a, position_b, cuboid_b),
      //
      // capsule vs. cuboid
      //
      ( shape::Variant::Bounded (shape::Bounded::Capsule (ref capsule_a)),
        shape::Variant::Bounded (shape::Bounded::Cuboid  (ref cuboid_b))
      ) => Self::query_capsule_cuboid (
        position_a, capsule_a, position_b, cuboid_b),
      ( shape::Variant::Bounded (shape::Bounded::Cuboid  (ref cuboid_a)),
        shape::Variant::Bounded (shape::Bounded::Capsule (ref capsule_b))
      ) => {
        let mut proximity = Self::query_capsule_cuboid (
          position_b, capsule_b, position_a, cuboid_a);
        proximity.half_axis *= -1.0;
        proximity.normal    = proximity.normal.invert();
        proximity
      }
      //
      // capsule vs. orthant
      //
      ( shape::Variant::Bounded   (shape::Bounded::Capsule   (ref capsule_a)),
        shape::Variant::Unbounded (shape::Unbounded::Orthant (ref orthant_b))
      ) => Self::query_capsule_orthant (
        position_a, capsule_a, position_b, orthant_b),
      ( shape::Variant::Unbounded (shape::Unbounded::Orthant (ref orthant_a)),
        shape::Variant::Bounded   (shape::Bounded::Capsule   (ref capsule_b))
      ) => {
        let mut proximity = Self::query_capsule_orthant (
          position_b, capsule_b, position_a, orthant_a);
        proximity.half_axis *= -1.0;
        proximity.normal    = proximity.normal.invert();
        proximity
      }
      _ => unimplemented!("proximity query not implemented for (A,B): {:?}",
        (shape_a, shape_b))
    }
  }

  pub fn relation (&self) -> Relation {
    if self.distance < 0.0 {
      Relation::Intersect
    } else if self.distance <= collision::CONTACT_DISTANCE {
      Relation::Contact
    } else {
      debug_assert!(self.distance > collision::CONTACT_DISTANCE);
      Relation::Disjoint
    }
  }

  #[inline]
  pub fn constraint_planar (&self) -> constraint::Planar {
    constraint::Planar::new (self.midpoint, self.normal)
  }

  pub fn query_capsule_capsule (
    position_a : &cgmath::Point3 <f64>, capsule_a  : &shape::Capsule <f64>,
    position_b : &cgmath::Point3 <f64>, capsule_b  : &shape::Capsule <f64>
  ) -> Self {
    use cgmath::{EuclideanSpace, InnerSpace, Zero};
    let (half_height_a, half_height_b) =
      (*capsule_a.half_height, *capsule_b.half_height);
    let (radius_a, radius_b) = (*capsule_a.radius, *capsule_b.radius);
    let shaft_max_a = *position_a + cgmath::vec3 (0.0, 0.0, half_height_a);
    let shaft_min_a = *position_a - cgmath::vec3 (0.0, 0.0, half_height_a);
    let shaft_max_b = *position_b + cgmath::vec3 (0.0, 0.0, half_height_b);
    let shaft_min_b = *position_b - cgmath::vec3 (0.0, 0.0, half_height_b);
    let distance_normal_half_axis = | axis : &cgmath::Vector3 <f64> | {
      let axis_distance = axis.magnitude();
      let distance      = axis_distance - radius_a - radius_b;
      let normal        = if axis.is_zero() {
        math::Unit3::axis_x()
      } else {
        math::Unit3::unchecked (-(axis / axis_distance))
      };
      let half_axis     = -0.5 * distance * *normal;
      (distance, normal, half_axis)
    };
    if shaft_max_a.z <= shaft_min_b.z {
      // cylinder shaft of B is above cylinder shaft of A
      let axis     = shaft_min_b - shaft_max_a;
      let (distance, normal, half_axis) = distance_normal_half_axis (&axis);
      let midpoint = shaft_max_a - radius_a * *normal + half_axis;
      Proximity { distance, half_axis, midpoint, normal }
    } else if shaft_max_b.z <= shaft_min_a.z {
      // cylinder shaft of B is below cylinder shaft of A
      let axis     = shaft_max_b.to_vec() - shaft_min_a.to_vec();
      let (distance, normal, half_axis) = distance_normal_half_axis (&axis);
      let midpoint = shaft_min_a - radius_a * *normal + half_axis;
      Proximity { distance, half_axis, midpoint, normal }
    } else {
      // cylinder shafts are adjacent
      let axis     =
        cgmath::Point3::from ([position_b.x, position_b.y, 0.0]) -
        cgmath::Point3::from ([position_a.x, position_a.y, 0.0]);
      let (distance, normal, half_axis) = distance_normal_half_axis (&axis);
      let midpoint = {
        let max_z  = f64::min (shaft_max_a.z, shaft_max_b.z);
        let min_z  = f64::max (shaft_min_a.z, shaft_min_b.z);
        let mid_z  = 0.5 * (max_z + min_z);
        let mid_a  = position_a - radius_a * *normal + half_axis;
        [mid_a.x, mid_a.y, mid_z].into()
      };
      Proximity { distance, half_axis, midpoint, normal }
    }
  }

  pub fn query_cuboid_cuboid (
    position_a : &cgmath::Point3 <f64>, cuboid_a : &shape::Cuboid <f64>,
    position_b : &cgmath::Point3 <f64>, cuboid_b : &shape::Cuboid <f64>
  ) -> Self {
    use cgmath::{EuclideanSpace, InnerSpace, Zero};

    fn nearest_ab_to_distance (
      nearest_a      : &cgmath::Point3  <f64>,
      nearest_b      : &cgmath::Point3  <f64>,
      default_normal : &math::Unit3 <f64>
    ) -> Proximity {
      let axis      = nearest_b - nearest_a;
      let distance  = axis.magnitude();
      let half_axis = 0.5 * (nearest_b - nearest_a);
      let midpoint  = nearest_a + half_axis;
      let normal    = if axis.is_zero() {
        *default_normal
      } else {
        math::Unit3::unchecked (-axis / distance)
      };
      Proximity { distance, half_axis, midpoint, normal }
    }

    let max_a = position_a + cuboid_a.max().to_vec();
    let min_a = position_a + cuboid_a.min().to_vec();
    let max_b = position_b + cuboid_b.max().to_vec();
    let min_b = position_b + cuboid_b.min().to_vec();
    let overlap_x = min_b.x < max_a.x && min_a.x < max_b.x;
    let overlap_y = min_b.y < max_a.y && min_a.y < max_b.y;
    let overlap_z = min_b.z < max_a.z && min_a.z < max_b.z;
    match (overlap_x, overlap_y, overlap_z) {
      (false, false, false) => {
        use cgmath::{ElementWise, EuclideanSpace};
        // no overlaps: nearest features are corners
        let a_to_b = position_b - position_a;
        debug_assert!(!a_to_b.x.is_zero());
        debug_assert!(!a_to_b.y.is_zero());
        debug_assert!(!a_to_b.z.is_zero());
        let a_to_b_sigvec = a_to_b.map (f64::signum);
        let nearest_a = position_a +
          cuboid_a.max().to_vec().mul_element_wise (a_to_b_sigvec);
        let nearest_b = position_b +
          cuboid_b.max().to_vec().mul_element_wise (-a_to_b_sigvec);
        nearest_ab_to_distance (&nearest_a, &nearest_b,
          &math::Unit3::normalize (-a_to_b_sigvec))
      }
      (true, false, false)  => {
        // overlap x only: nearest features are edges
        let max_x  = f64::min (max_a.x, max_b.x);
        let min_x  = f64::max (min_a.x, min_b.x);
        debug_assert_ne!(min_x, max_x);
        let mid_x  = 0.5 * (min_x + max_x);
        let a_to_b = position_b - position_a;
        let a_to_b_yz = cgmath::vec3 (0.0, a_to_b.y, a_to_b.z);
        let a_to_b_yz_sigvec = math::sigvec (a_to_b_yz);
        debug_assert!(a_to_b_yz_sigvec.x.is_zero());
        debug_assert!(!a_to_b_yz_sigvec.y.is_zero());
        debug_assert!(!a_to_b_yz_sigvec.z.is_zero());
        let nearest_a = {
          let y = *cuboid_a.half_extent_y * a_to_b_yz_sigvec.y;
          let z = *cuboid_a.half_extent_z * a_to_b_yz_sigvec.z;
          [mid_x, position_a.y + y, position_a.z + z].into()
        };
        let nearest_b = {
          let y = -*cuboid_b.half_extent_y * a_to_b_yz_sigvec.y;
          let z = -*cuboid_b.half_extent_z * a_to_b_yz_sigvec.z;
          [mid_x, position_b.y + y, position_b.z + z].into()
        };
        nearest_ab_to_distance (&nearest_a, &nearest_b,
          &math::Unit3::normalize (-a_to_b_yz_sigvec))
      }
      (false, true, false)  => {
        // overlap y only: nearest features are edges
        let max_y  = f64::min (max_a.y, max_b.y);
        let min_y  = f64::max (min_a.y, min_b.y);
        debug_assert_ne!(min_y, max_y);
        let mid_y  = 0.5 * (min_y + max_y);
        let a_to_b = position_b - position_a;
        let a_to_b_xz = cgmath::vec3 (a_to_b.x, 0.0, a_to_b.z);
        let a_to_b_xz_sigvec = math::sigvec (a_to_b_xz);
        debug_assert!(!a_to_b_xz_sigvec.x.is_zero());
        debug_assert!(a_to_b_xz_sigvec.y.is_zero());
        debug_assert!(!a_to_b_xz_sigvec.z.is_zero());
        let nearest_a = {
          let x = *cuboid_a.half_extent_x * a_to_b_xz_sigvec.x;
          let z = *cuboid_a.half_extent_z * a_to_b_xz_sigvec.z;
          [position_a.x + x, mid_y, position_a.z + z].into()
        };
        let nearest_b = {
          let x = -*cuboid_b.half_extent_x * a_to_b_xz_sigvec.x;
          let z = -*cuboid_b.half_extent_z * a_to_b_xz_sigvec.z;
          [position_b.x + x, mid_y, position_b.z + z].into()
        };
        nearest_ab_to_distance (&nearest_a, &nearest_b,
          &math::Unit3::normalize (-a_to_b_xz_sigvec))
      }
      (false, false, true)  => {
        // overlap z only: nearest features are edges
        let max_z  = f64::min (max_a.z, max_b.z);
        let min_z  = f64::max (min_a.z, min_b.z);
        debug_assert_ne!(min_z, max_z);
        let mid_z  = 0.5 * (min_z + max_z);
        let a_to_b = position_b - position_a;
        let a_to_b_xy = cgmath::vec3 (a_to_b.x, a_to_b.y, 0.0);
        let a_to_b_xy_sigvec = math::sigvec (a_to_b_xy);
        debug_assert!(!a_to_b_xy_sigvec.x.is_zero());
        debug_assert!(!a_to_b_xy_sigvec.y.is_zero());
        debug_assert!(a_to_b_xy_sigvec.z.is_zero());
        let nearest_a = {
          let x = *cuboid_a.half_extent_x * a_to_b_xy_sigvec.x;
          let y = *cuboid_a.half_extent_y * a_to_b_xy_sigvec.y;
          [position_a.x + x, position_a.y + y, mid_z].into()
        };
        let nearest_b = {
          let x = -*cuboid_b.half_extent_x * a_to_b_xy_sigvec.x;
          let y = -*cuboid_b.half_extent_y * a_to_b_xy_sigvec.y;
          [position_b.x + x, position_b.y + y, mid_z].into()
        };
        nearest_ab_to_distance (&nearest_a, &nearest_b,
          &math::Unit3::normalize (-a_to_b_xy_sigvec))
      }
      (true, true, false)   => {
        // overlap x and y: nearest features are faces
        let max_x    = f64::min (max_a.x, max_b.x);
        let min_x    = f64::max (min_a.x, min_b.x);
        let mid_x    = 0.5 * (min_x + max_x);
        let max_y    = f64::min (max_a.y, max_b.y);
        let min_y    = f64::max (min_a.y, min_b.y);
        let mid_y    = 0.5 * (min_y + max_y);
        let a_to_b_z = position_b.z - position_a.z;
        let a_to_b_signum_z = a_to_b_z.signum();
        debug_assert_ne!(a_to_b_signum_z, 0.0);
        let nearest_a = cgmath::Point3::new (
          mid_x, mid_y,
          position_a.z + a_to_b_signum_z * *cuboid_a.half_extent_z
        );
        let nearest_b = cgmath::Point3::new (
          mid_x, mid_y,
          position_b.z - a_to_b_signum_z * *cuboid_b.half_extent_z
        );
        let nearest_a_to_b_z = nearest_b.z - nearest_a.z;
        let distance  = nearest_a_to_b_z.abs();
        let half_axis = [0.0, 0.0, 0.5 * nearest_a_to_b_z].into();
        let normal    = math::Unit3::unchecked (
          [0.0, 0.0, -a_to_b_signum_z].into());
        let midpoint  = nearest_a + half_axis;
        Proximity { distance, half_axis, midpoint, normal }
      }
      (true, false, true)   => {
        // overlap x and z: nearest features are faces
        let max_x    = f64::min (max_a.x, max_b.x);
        let min_x    = f64::max (min_a.x, min_b.x);
        let mid_x    = 0.5 * (min_x + max_x);
        let max_z    = f64::min (max_a.z, max_b.z);
        let min_z    = f64::max (min_a.z, min_b.z);
        let mid_z    = 0.5 * (min_z + max_z);
        let a_to_b_y = position_b.y - position_a.y;
        let a_to_b_signum_y = a_to_b_y.signum();
        debug_assert_ne!(a_to_b_signum_y, 0.0);
        let nearest_a = cgmath::Point3::new (
          mid_x,
          position_a.y + a_to_b_signum_y * *cuboid_a.half_extent_y,
          mid_z
        );
        let nearest_b = cgmath::Point3::new (
          mid_x,
          position_b.y - a_to_b_signum_y * *cuboid_b.half_extent_y,
          mid_z
        );
        let nearest_a_to_b_y = nearest_b.y - nearest_a.y;
        let distance  = nearest_a_to_b_y.abs();
        let half_axis = [0.0, 0.5 * nearest_a_to_b_y, 0.0].into();
        let normal    = math::Unit3::unchecked (
          [0.0, -a_to_b_signum_y, 0.0].into());
        let midpoint  = nearest_a + half_axis;
        Proximity { distance, half_axis, midpoint, normal }
      }
      (false, true, true)   => {
        // overlap y and z: nearest features are faces
        let max_y    = f64::min (max_a.y, max_b.y);
        let min_y    = f64::max (min_a.y, min_b.y);
        let mid_y    = 0.5 * (min_y + max_y);
        let max_z    = f64::min (max_a.z, max_b.z);
        let min_z    = f64::max (min_a.z, min_b.z);
        let mid_z    = 0.5 * (min_z + max_z);
        let a_to_b_x = position_b.x - position_a.x;
        let a_to_b_signum_x = a_to_b_x.signum();
        debug_assert_ne!(a_to_b_signum_x, 0.0);
        let nearest_a = cgmath::Point3::new (
          position_a.x + a_to_b_signum_x * *cuboid_a.half_extent_x,
          mid_y,
          mid_z
        );
        let nearest_b = cgmath::Point3::new (
          position_b.x - a_to_b_signum_x * *cuboid_b.half_extent_x,
          mid_y,
          mid_z
        );
        let nearest_a_to_b_x = nearest_b.x - nearest_a.x;
        let distance  = nearest_a_to_b_x.abs();
        let half_axis = [0.5 * nearest_a_to_b_x, 0.0, 0.0].into();
        let normal    = math::Unit3::unchecked (
          [-a_to_b_signum_x, 0.0, 0.0].into());
        let midpoint  = nearest_a + half_axis;
        Proximity { distance, half_axis, midpoint, normal }
      }
      (true, true, true)    => {
        // overlap on all axes: intersection
        let max_x     = f64::min (max_a.x, max_b.x);
        let min_x     = f64::max (min_a.x, min_b.x);
        let mid_x     = 0.5 * (min_x + max_x);
        let max_y     = f64::min (max_a.y, max_b.y);
        let min_y     = f64::max (min_a.y, min_b.y);
        let mid_y     = 0.5 * (min_y + max_y);
        let max_z     = f64::min (max_a.z, max_b.z);
        let min_z     = f64::max (min_a.z, min_b.z);
        let mid_z     = 0.5 * (min_z + max_z);
        let midpoint  = [mid_x, mid_y, mid_z].into();
        let resolve_x = if position_a.x < position_b.x {
          min_b.x - max_a.x
        } else {
          debug_assert!(position_a.x >= position_b.x);
          max_b.x - min_a.x
        };
        let resolve_y = if position_a.y < position_b.y {
          min_b.y - max_a.y
        } else {
          debug_assert!(position_a.y >= position_b.y);
          max_b.y - min_a.y
        };
        let resolve_z = if position_a.z < position_b.z {
          min_b.z - max_a.z
        } else {
          debug_assert!(position_a.z >= position_b.z);
          max_b.z - min_a.z
        };
        let resolve_x_abs = resolve_x.abs();
        let resolve_y_abs = resolve_y.abs();
        let resolve_z_abs = resolve_z.abs();
        let (distance, half_axis, normal) : (
          f64, cgmath::Vector3 <f64>, math::Unit3 <f64>
        ) = if resolve_x_abs <= resolve_y_abs && resolve_x_abs <= resolve_z_abs {
          ( -resolve_x_abs,
            0.5 * cgmath::vec3 (resolve_x, 0.0, 0.0),
            math::Unit3::unchecked ([resolve_x.signum(), 0.0, 0.0].into())
          )
        } else if
          resolve_y_abs <= resolve_x_abs && resolve_y_abs <= resolve_z_abs
        {
          ( -resolve_y_abs,
            0.5 * cgmath::vec3 (0.0, resolve_y, 0.0),
            math::Unit3::unchecked ([0.0, resolve_y.signum(), 0.0].into())
          )
        } else {
          debug_assert!(resolve_z_abs <= resolve_x_abs);
          debug_assert!(resolve_z_abs <= resolve_y_abs);
          ( -resolve_z_abs,
            0.5 * cgmath::vec3 (0.0, 0.0, resolve_z),
            math::Unit3::unchecked ([0.0, 0.0, resolve_z.signum()].into())
          )
        };
        debug_assert!(!half_axis.is_zero());
        Proximity { distance, half_axis, midpoint, normal }
      }
    }
  }

  /*
  /// Uses pre-computed sigvec normals with array lookup.
  ///
  /// This benchmarked 2-3ns slower than the version using `normalize`.
  pub fn query_cuboid_cuboid_array (
    position_a : &cgmath::Point3 <f64>, cuboid_a : &shape::Cuboid <f64>,
    position_b : &cgmath::Point3 <f64>, cuboid_b : &shape::Cuboid <f64>
  ) -> Self {
    use cgmath::{EuclideanSpace, InnerSpace, Zero};

    fn nearest_ab_to_distance (
      nearest_a      : &cgmath::Point3  <f64>,
      nearest_b      : &cgmath::Point3  <f64>,
      default_normal : &cgmath::Vector3 <f64>
    ) -> Proximity {
      let axis      = nearest_b - nearest_a;
      let distance  = axis.magnitude();
      let half_axis = 0.5 * (nearest_b - nearest_a);
      let midpoint  = nearest_a + half_axis;
      let normal    = if axis.is_zero() {
        *default_normal
      } else {
        -axis / distance
      };
      Proximity { distance, half_axis, midpoint, normal }
    }

    let max_a = position_a + cuboid_a.max().to_vec();
    let min_a = position_a + cuboid_a.min().to_vec();
    let max_b = position_b + cuboid_b.max().to_vec();
    let min_b = position_b + cuboid_b.min().to_vec();
    let overlap_x = min_b.x < max_a.x && min_a.x < max_b.x;
    let overlap_y = min_b.y < max_a.y && min_a.y < max_b.y;
    let overlap_z = min_b.z < max_a.z && min_a.z < max_b.z;
    match (overlap_x, overlap_y, overlap_z) {
      (false, false, false) => {
        use cgmath::ElementWise;
        // no overlaps: nearest features are corners
        let a_to_b = position_b - position_a;
        debug_assert!(!a_to_b.x.is_zero());
        debug_assert!(!a_to_b.y.is_zero());
        debug_assert!(!a_to_b.z.is_zero());
        let a_to_b_sigvec = a_to_b.map (f64::signum);
        let nearest_a = position_a +
          cuboid_a.half_extents().mul_element_wise (a_to_b_sigvec);
        let nearest_b = position_b +
          cuboid_b.half_extents().mul_element_wise (-a_to_b_sigvec);
        nearest_ab_to_distance (&nearest_a, &nearest_b,
          &-math::sigvec_unit_f64_unchecked (a_to_b_sigvec))
      }
      (true, false, false)  => {
        // overlap x only: nearest features are edges
        let max_x  = f64::min (max_a.x, max_b.x);
        let min_x  = f64::max (min_a.x, min_b.x);
        debug_assert_ne!(min_x, max_x);
        let mid_x  = 0.5 * (min_x + max_x);
        let a_to_b = position_b - position_a;
        let a_to_b_yz = cgmath::vec3 (0.0, a_to_b.y, a_to_b.z);
        let a_to_b_yz_sigvec = math::sigvec (a_to_b_yz);
        debug_assert!(a_to_b_yz_sigvec.x.is_zero());
        debug_assert!(!a_to_b_yz_sigvec.y.is_zero());
        debug_assert!(!a_to_b_yz_sigvec.z.is_zero());
        let nearest_a = {
          let y = cuboid_a.half_extent_y * a_to_b_yz_sigvec.y;
          let z = cuboid_a.half_extent_z * a_to_b_yz_sigvec.z;
          [mid_x, position_a.y + y, position_a.z + z].into()
        };
        let nearest_b = {
          let y = -cuboid_b.half_extent_y * a_to_b_yz_sigvec.y;
          let z = -cuboid_b.half_extent_z * a_to_b_yz_sigvec.z;
          [mid_x, position_b.y + y, position_b.z + z].into()
        };
        nearest_ab_to_distance (&nearest_a, &nearest_b,
          &-math::sigvec_unit_f64_unchecked (a_to_b_yz_sigvec))
      }
      (false, true, false)  => {
        // overlap y only: nearest features are edges
        let max_y  = f64::min (max_a.y, max_b.y);
        let min_y  = f64::max (min_a.y, min_b.y);
        debug_assert_ne!(min_y, max_y);
        let mid_y  = 0.5 * (min_y + max_y);
        let a_to_b = position_b - position_a;
        let a_to_b_xz = cgmath::vec3 (a_to_b.x, 0.0, a_to_b.z);
        let a_to_b_xz_sigvec = math::sigvec (a_to_b_xz);
        debug_assert!(!a_to_b_xz_sigvec.x.is_zero());
        debug_assert!(a_to_b_xz_sigvec.y.is_zero());
        debug_assert!(!a_to_b_xz_sigvec.z.is_zero());
        let nearest_a = {
          let x = cuboid_a.half_extent_x * a_to_b_xz_sigvec.x;
          let z = cuboid_a.half_extent_z * a_to_b_xz_sigvec.z;
          [position_a.x + x, mid_y, position_a.z + z].into()
        };
        let nearest_b = {
          let x = -cuboid_b.half_extent_x * a_to_b_xz_sigvec.x;
          let z = -cuboid_b.half_extent_z * a_to_b_xz_sigvec.z;
          [position_b.x + x, mid_y, position_b.z + z].into()
        };
        nearest_ab_to_distance (&nearest_a, &nearest_b,
          &-math::sigvec_unit_f64_unchecked (a_to_b_xz_sigvec))
      }
      (false, false, true)  => {
        // overlap z only: nearest features are edges
        let max_z  = f64::min (max_a.z, max_b.z);
        let min_z  = f64::max (min_a.z, min_b.z);
        debug_assert_ne!(min_z, max_z);
        let mid_z  = 0.5 * (min_z + max_z);
        let a_to_b = position_b - position_a;
        let a_to_b_xy = cgmath::vec3 (a_to_b.x, a_to_b.y, 0.0);
        let a_to_b_xy_sigvec = math::sigvec (a_to_b_xy);
        debug_assert!(!a_to_b_xy_sigvec.x.is_zero());
        debug_assert!(!a_to_b_xy_sigvec.y.is_zero());
        debug_assert!(a_to_b_xy_sigvec.z.is_zero());
        let nearest_a = {
          let x = cuboid_a.half_extent_x * a_to_b_xy_sigvec.x;
          let y = cuboid_a.half_extent_y * a_to_b_xy_sigvec.y;
          [position_a.x + x, position_a.y + y, mid_z].into()
        };
        let nearest_b = {
          let x = -cuboid_b.half_extent_x * a_to_b_xy_sigvec.x;
          let y = -cuboid_b.half_extent_y * a_to_b_xy_sigvec.y;
          [position_b.x + x, position_b.y + y, mid_z].into()
        };
        nearest_ab_to_distance (&nearest_a, &nearest_b,
          &-math::sigvec_unit_f64_unchecked (a_to_b_xy_sigvec))
      }
      (true, true, false)   => {
        // overlap x and y: nearest features are faces
        let max_x    = f64::min (max_a.x, max_b.x);
        let min_x    = f64::max (min_a.x, min_b.x);
        let mid_x    = 0.5 * (min_x + max_x);
        let max_y    = f64::min (max_a.y, max_b.y);
        let min_y    = f64::max (min_a.y, min_b.y);
        let mid_y    = 0.5 * (min_y + max_y);
        let a_to_b_z = position_b.z - position_a.z;
        let a_to_b_signum_z = a_to_b_z.signum();
        debug_assert_ne!(a_to_b_signum_z, 0.0);
        let nearest_a = cgmath::Point3::new (
          mid_x, mid_y,
          position_a.z + a_to_b_signum_z * cuboid_a.half_extent_z
        );
        let nearest_b = cgmath::Point3::new (
          mid_x, mid_y,
          position_b.z - a_to_b_signum_z * cuboid_b.half_extent_z
        );
        let nearest_a_to_b_z = nearest_b.z - nearest_a.z;
        let distance  = nearest_a_to_b_z.abs();
        let half_axis = [0.0, 0.0, 0.5 * nearest_a_to_b_z].into();
        let normal    = [0.0, 0.0, -a_to_b_signum_z].into();
        let midpoint  = nearest_a + half_axis;
        Proximity { distance, half_axis, midpoint, normal }
      }
      (true, false, true)   => {
        // overlap x and z: nearest features are faces
        let max_x    = f64::min (max_a.x, max_b.x);
        let min_x    = f64::max (min_a.x, min_b.x);
        let mid_x    = 0.5 * (min_x + max_x);
        let max_z    = f64::min (max_a.z, max_b.z);
        let min_z    = f64::max (min_a.z, min_b.z);
        let mid_z    = 0.5 * (min_z + max_z);
        let a_to_b_y = position_b.y - position_a.y;
        let a_to_b_signum_y = a_to_b_y.signum();
        debug_assert_ne!(a_to_b_signum_y, 0.0);
        let nearest_a = cgmath::Point3::new (
          mid_x,
          position_a.y + a_to_b_signum_y * cuboid_a.half_extent_y,
          mid_z
        );
        let nearest_b = cgmath::Point3::new (
          mid_x,
          position_b.y - a_to_b_signum_y * cuboid_b.half_extent_y,
          mid_z
        );
        let nearest_a_to_b_y = nearest_b.y - nearest_a.y;
        let distance  = nearest_a_to_b_y.abs();
        let half_axis = [0.0, 0.5 * nearest_a_to_b_y, 0.0].into();
        let normal    = [0.0, -a_to_b_signum_y, 0.0].into();
        let midpoint  = nearest_a + half_axis;
        Proximity { distance, half_axis, midpoint, normal }
      }
      (false, true, true)   => {
        // overlap y and z: nearest features are faces
        let max_y    = f64::min (max_a.y, max_b.y);
        let min_y    = f64::max (min_a.y, min_b.y);
        let mid_y    = 0.5 * (min_y + max_y);
        let max_z    = f64::min (max_a.z, max_b.z);
        let min_z    = f64::max (min_a.z, min_b.z);
        let mid_z    = 0.5 * (min_z + max_z);
        let a_to_b_x = position_b.x - position_a.x;
        let a_to_b_signum_x = a_to_b_x.signum();
        debug_assert_ne!(a_to_b_signum_x, 0.0);
        let nearest_a = cgmath::Point3::new (
          position_a.x + a_to_b_signum_x * cuboid_a.half_extent_x,
          mid_y,
          mid_z
        );
        let nearest_b = cgmath::Point3::new (
          position_b.x - a_to_b_signum_x * cuboid_b.half_extent_x,
          mid_y,
          mid_z
        );
        let nearest_a_to_b_x = nearest_b.x - nearest_a.x;
        let distance  = nearest_a_to_b_x.abs();
        let half_axis = [0.5 * nearest_a_to_b_x, 0.0, 0.0].into();
        let normal    = [-a_to_b_signum_x, 0.0, 0.0].into();
        let midpoint  = nearest_a + half_axis;
        Proximity { distance, half_axis, midpoint, normal }
      }
      (true, true, true)    => {
        // overlap on all axes: intersection
        let max_x     = f64::min (max_a.x, max_b.x);
        let min_x     = f64::max (min_a.x, min_b.x);
        let mid_x     = 0.5 * (min_x + max_x);
        let max_y     = f64::min (max_a.y, max_b.y);
        let min_y     = f64::max (min_a.y, min_b.y);
        let mid_y     = 0.5 * (min_y + max_y);
        let max_z     = f64::min (max_a.z, max_b.z);
        let min_z     = f64::max (min_a.z, min_b.z);
        let mid_z     = 0.5 * (min_z + max_z);
        let midpoint  = [mid_x, mid_y, mid_z].into();
        let resolve_x = if position_a.x < position_b.x {
          min_b.x - max_a.x
        } else {
          debug_assert!(position_a.x >= position_b.x);
          max_b.x - min_a.x
        };
        let resolve_y = if position_a.y < position_b.y {
          min_b.y - max_a.y
        } else {
          debug_assert!(position_a.y >= position_b.y);
          max_b.y - min_a.y
        };
        let resolve_z = if position_a.z < position_b.z {
          min_b.z - max_a.z
        } else {
          debug_assert!(position_a.z >= position_b.z);
          max_b.z - min_a.z
        };
        let resolve_x_abs = resolve_x.abs();
        let resolve_y_abs = resolve_y.abs();
        let resolve_z_abs = resolve_z.abs();
        let (distance, half_axis, normal) : (
          f64, cgmath::Vector3 <f64>, cgmath::Vector3 <f64>
        ) = if resolve_x_abs <= resolve_y_abs && resolve_x_abs <= resolve_z_abs {
          ( -resolve_x_abs,
            0.5 * cgmath::vec3 (resolve_x, 0.0, 0.0),
            [resolve_x.signum(), 0.0, 0.0].into()
          )
        } else if
          resolve_y_abs <= resolve_x_abs && resolve_y_abs <= resolve_z_abs
        {
          ( -resolve_y_abs,
            0.5 * cgmath::vec3 (0.0, resolve_y, 0.0),
            [0.0, resolve_y.signum(), 0.0].into()
          )
        } else {
          debug_assert!(resolve_z_abs <= resolve_x_abs);
          debug_assert!(resolve_z_abs <= resolve_y_abs);
          ( -resolve_z_abs,
            0.5 * cgmath::vec3 (0.0, 0.0, resolve_z),
            [0.0, 0.0, resolve_z.signum()].into()
          )
        };
        debug_assert!(!half_axis.is_zero());
        Proximity { distance, half_axis, midpoint, normal }
      }
    }
  } // end query_cuboid_cuboid_array

  /// This is a version of the cuboid distance query with the three edge
  /// cases factored out into a function taking vector indices. This turns out
  /// to be slower in benchmarks (190ns vs. the 180ns inline versions above).
  pub fn query_cuboid_cuboid_refactor (
    position_a : &cgmath::Point3 <f64>, cuboid_a : &shape::Cuboid <f64>,
    position_b : &cgmath::Point3 <f64>, cuboid_b : &shape::Cuboid <f64>
  ) -> Self {
    use cgmath::{EuclideanSpace, InnerSpace, Zero};

    #[inline]
    fn distance_edge (
      edge : usize, other1 : usize, other2 : usize,
      position_a : &cgmath::Point3 <f64>, cuboid_a : &shape::Cuboid <f64>,
      position_b : &cgmath::Point3 <f64>, cuboid_b : &shape::Cuboid <f64>,
      max_a      : &cgmath::Point3 <f64>, min_a    : &cgmath::Point3 <f64>,
      max_b      : &cgmath::Point3 <f64>, min_b    : &cgmath::Point3 <f64>
    ) -> Proximity {
      let max_edge  = f64::min (max_a[edge], max_b[edge]);
      let min_edge  = f64::max (min_a[edge], min_b[edge]);
      debug_assert_ne!(min_edge, max_edge);
      let mid_edge  = 0.5 * (min_edge + max_edge);
      let a_to_b = position_b - position_a;
      let a_to_b_others = {
        let mut a_to_b_others = a_to_b;
        a_to_b_others[edge] = 0.0;
        a_to_b_others
      };
      let a_to_b_others_signum = a_to_b_others.map (math::signum_or_zero);
      debug_assert!(a_to_b_others_signum[edge].is_zero());
      debug_assert!(!a_to_b_others_signum[other1].is_zero());
      debug_assert!(!a_to_b_others_signum[other2].is_zero());
      let nearest_a = {
        let other_component1 =
          cuboid_a.half_extents()[other1] * a_to_b_others_signum[other1];
        let other_component2 =
          cuboid_a.half_extents()[other2] * a_to_b_others_signum[other2];
        let mut nearest_a = [0.0; 3];
        nearest_a[edge]   = mid_edge;
        nearest_a[other1] = position_a[other1] + other_component1;
        nearest_a[other2] = position_a[other2] + other_component2;
        nearest_a.into()
      };
      let nearest_b = {
        let other_component1 =
          -cuboid_b.half_extents()[other1] * a_to_b_others_signum[other1];
        let other_component2 =
          -cuboid_b.half_extents()[other2] * a_to_b_others_signum[other2];
        let mut nearest_b = [0.0; 3];
        nearest_b[edge]   = mid_edge;
        nearest_b[other1] = position_b[other1] + other_component1;
        nearest_b[other2] = position_b[other2] + other_component2;
        nearest_b.into()
      };
      nearest_ab_to_distance (&nearest_a, &nearest_b, &a_to_b_others_signum)
    }

    fn nearest_ab_to_distance (
      nearest_a     : &cgmath::Point3 <f64>,
      nearest_b     : &cgmath::Point3 <f64>,
      a_to_b_signum : &cgmath::Vector3 <f64>
    ) -> Proximity {
      let axis      = nearest_b - nearest_a;
      let distance  = axis.magnitude();
      let half_axis = 0.5 * (nearest_b - nearest_a);
      let midpoint  = nearest_a + half_axis;
      let normal    = if axis.is_zero() {
        -a_to_b_signum.normalize()
      } else {
        -axis / distance
      };
      Proximity { distance, half_axis, midpoint, normal }
    }

    let max_a = position_a + cuboid_a.max().to_vec();
    let min_a = position_a + cuboid_a.min().to_vec();
    let max_b = position_b + cuboid_b.max().to_vec();
    let min_b = position_b + cuboid_b.min().to_vec();

    let overlap_x = min_b.x < max_a.x && min_a.x < max_b.x;
    let overlap_y = min_b.y < max_a.y && min_a.y < max_b.y;
    let overlap_z = min_b.z < max_a.z && min_a.z < max_b.z;
    match (overlap_x, overlap_y, overlap_z) {
      (false, false, false) => {
        use cgmath::ElementWise;
        // no overlaps: nearest features are corners
        let a_to_b = position_b - position_a;
        debug_assert!(!a_to_b.x.is_zero());
        debug_assert!(!a_to_b.y.is_zero());
        debug_assert!(!a_to_b.z.is_zero());
        let a_to_b_signum = a_to_b.map (f64::signum);
        let nearest_a = position_a +
          cuboid_a.half_extents().mul_element_wise (a_to_b_signum);
        let nearest_b = position_b +
          cuboid_b.half_extents().mul_element_wise (-a_to_b_signum);
        nearest_ab_to_distance (&nearest_a, &nearest_b, &a_to_b_signum)
      }
      (true, false, false)  => {
        // overlap x only: nearest features are edges
        distance_edge (0, 1, 2,
          position_a, cuboid_a, position_b, cuboid_b,
          &max_a, &min_a, &max_b, &min_b
        )
      }
      (false, true, false)  => {
        // overlap y only: nearest features are edges
        distance_edge (1, 0, 2,
          position_a, cuboid_a, position_b, cuboid_b,
          &max_a, &min_a, &max_b, &min_b
        )
      }
      (false, false, true)  => {
        // overlap z only: nearest features are edges
        distance_edge (2, 1, 2,
          position_a, cuboid_a, position_b, cuboid_b,
          &max_a, &min_a, &max_b, &min_b
        )
      }
      (true, true, false)   => {
        // overlap x and y: nearest features are faces
        unimplemented!()
      }
      (true, false, true)   => {
        // overlap x and z: nearest features are faces
        unimplemented!()
      }
      (false, true, true)   => {
        // overlap y and z: nearest features are faces
        unimplemented!()
      }
      (true, true, true)    => {
        // overlap on all axes: intersection
        unimplemented!()
      }
    }
  } // end query_cuboid_cuboid_refactor
  */

  pub fn query_capsule_cuboid (
    position_a : &cgmath::Point3 <f64>, capsule_a : &shape::Capsule <f64>,
    position_b : &cgmath::Point3 <f64>, cuboid_b  : &shape::Cuboid  <f64>
  ) -> Self {
    use cgmath::{ElementWise, EuclideanSpace, InnerSpace, Zero};

    let shaft_max_a  = *position_a +
      cgmath::vec3 (0.0, 0.0, *capsule_a.half_height);
    let shaft_min_a  = *position_a -
      cgmath::vec3 (0.0, 0.0, *capsule_a.half_height);
    let cuboid_max_b = position_b + cuboid_b.max().to_vec();
    let cuboid_min_b = position_b - cuboid_b.max().to_vec();

    let overlap_x = cuboid_min_b.x < position_a.x &&
      position_a.x < cuboid_max_b.x;
    let overlap_y = cuboid_min_b.y < position_a.y &&
      position_a.y < cuboid_max_b.y;

    let capsule_a_radius       = *capsule_a.radius;
    let capsule_a_half_height  = *capsule_a.half_height;
    let cuboid_b_half_extent_x = *cuboid_b.half_extent_x;
    let cuboid_b_half_extent_y = *cuboid_b.half_extent_y;
    let cuboid_b_half_extent_z = *cuboid_b.half_extent_z;
    if shaft_min_a.z >= cuboid_max_b.z {
      // capsule above cuboid
      match (overlap_x, overlap_y) {
        (false, false) => {
          // nearest/deepest point is a corner of the cuboid
          let b_to_a        = shaft_min_a - position_b;
          debug_assert!(!b_to_a.x.is_zero());
          debug_assert!(!b_to_a.y.is_zero());
          debug_assert!(!b_to_a.z.is_zero());
          let b_to_a_sigvec = b_to_a.map (f64::signum);
          let nearest_b     = position_b +
            cuboid_b.max().to_vec().mul_element_wise (b_to_a_sigvec);
          let a_to_b        = nearest_b - shaft_min_a;
          let nearest_a     = shaft_min_a + if !a_to_b.is_zero() {
            capsule_a_radius * a_to_b.normalize()
          } else {
            [-b_to_a_sigvec.x * capsule_a_radius, 0.0, 0.0].into()
          };
          let axis          = nearest_b - nearest_a;
          let distance      = axis.magnitude() * if !a_to_b.is_zero() {
            a_to_b.dot (axis).signum()
          } else {
            -1.0
          };
          let half_axis     = 0.5 * axis;
          let normal        = if axis.is_zero() {
            math::Unit3::normalize (b_to_a_sigvec)
          } else {
            math::Unit3::unchecked (-axis / distance)
          };
          let midpoint      = nearest_a + half_axis;
          Proximity { distance, half_axis, midpoint, normal }
        }
        (true, false) => {
          // nearest/deepest point is on an edge parallel to the X axis
          let signum_y  = (shaft_min_a.y - position_b.y).signum();
          let nearest_b = cgmath::Point3::new (
            shaft_min_a.x,
            position_b.y + signum_y * cuboid_b_half_extent_y,
            cuboid_max_b.z
          );
          let a_to_b    = nearest_b - shaft_min_a;
          let nearest_a = shaft_min_a + if a_to_b.is_zero() {
            -signum_y * cgmath::vec3 (0.0, capsule_a_radius, 0.0)
          } else {
            capsule_a_radius * a_to_b.normalize()
          };
          let axis      = nearest_b - nearest_a;
          let distance  = axis.magnitude() * if !a_to_b.is_zero() {
            a_to_b.dot (axis).signum() * axis.magnitude()
          } else {
            -1.0
          };
          let half_axis = 0.5 * axis;
          let normal    = if axis.is_zero() {
            math::Unit3::normalize (cgmath::vec3 (0.0, signum_y, 1.0))
          } else {
            math::Unit3::unchecked (-axis / distance)
          };
          let midpoint  = nearest_a + half_axis;
          Proximity { distance, half_axis, midpoint, normal }
        }
        (false, true) => {
          // nearest/deepest point is on an edge parallel to the Y axis
          let signum_x  = (shaft_min_a.x - position_b.x).signum();
          let nearest_b = cgmath::Point3::new (
            position_b.x + signum_x * cuboid_b_half_extent_x,
            shaft_min_a.y,
            cuboid_max_b.z
          );
          let a_to_b    = nearest_b - shaft_min_a;
          let nearest_a = shaft_min_a + if a_to_b.is_zero() {
            -signum_x * cgmath::vec3 (capsule_a_radius, 0.0, 0.0)
          } else {
            capsule_a_radius * a_to_b.normalize()
          };
          let axis      = nearest_b - nearest_a;
          let distance  = axis.magnitude() * if !a_to_b.is_zero() {
            a_to_b.dot (axis).signum()
          } else {
            -1.0
          };
          let half_axis = 0.5 * axis;
          let normal    = if axis.is_zero() {
            math::Unit3::normalize (cgmath::vec3 (signum_x, 0.0, 1.0))
          } else {
            math::Unit3::unchecked (-axis / distance)
          };
          let midpoint  = nearest_a + half_axis;
          Proximity { distance, half_axis, midpoint, normal }
        }
        (true, true) => {
          // nearest/deepest point is in the -Z direction
          let nearest_a : cgmath::Point3 <f64> =
            shaft_min_a - cgmath::vec3 (0.0, 0.0, capsule_a_radius);
          let nearest_b : cgmath::Point3 <f64> =
            [position_a.x, position_a.y, cuboid_max_b.z].into();
          let axis      = nearest_b - nearest_a;
          let normal    = math::Unit3::axis_z();
          let half_axis = 0.5 * axis;
          let distance  = nearest_a.z - nearest_b.z;
          let midpoint  = nearest_a + half_axis;
          Proximity { distance, half_axis, midpoint, normal }
        }
      }
    } else if shaft_max_a.z <= cuboid_min_b.z {
      // capsule below cuboid
      match (overlap_x, overlap_y) {
        (false, false) => {
          // nearest/deepest point is a corner of the cuboid
          let b_to_a        = shaft_max_a - position_b;
          debug_assert!(!b_to_a.x.is_zero());
          debug_assert!(!b_to_a.y.is_zero());
          debug_assert!(!b_to_a.z.is_zero());
          let b_to_a_sigvec = b_to_a.map (f64::signum);
          let nearest_b     = position_b +
            cuboid_b.max().to_vec().mul_element_wise (b_to_a_sigvec);
          let a_to_b        = nearest_b - shaft_max_a;
          let nearest_a     = shaft_max_a + if !a_to_b.is_zero() {
            capsule_a_radius * a_to_b.normalize()
          } else {
            [-b_to_a_sigvec.x * capsule_a_radius, 0.0, 0.0].into()
          };
          let axis          = nearest_b - nearest_a;
          let distance      = axis.magnitude() * if !a_to_b.is_zero() {
            a_to_b.dot (axis).signum()
          } else {
            -1.0
          };
          let half_axis     = 0.5 * axis;
          let normal        = if axis.is_zero() {
            math::Unit3::normalize (b_to_a_sigvec)
          } else {
            math::Unit3::unchecked (-axis / distance)
          };
          let midpoint      = nearest_a + half_axis;
          Proximity { distance, half_axis, midpoint, normal }
        }
        (true, false) => {
          // nearest/deepest point is on an edge parallel to the X axis
          let signum_y  = (shaft_max_a.y - position_b.y).signum();
          let nearest_b = cgmath::Point3::new (
            shaft_max_a.x,
            position_b.y + signum_y * cuboid_b_half_extent_y,
            cuboid_min_b.z
          );
          let a_to_b    = nearest_b - shaft_max_a;
          let nearest_a = shaft_max_a + if a_to_b.is_zero() {
            -signum_y * cgmath::vec3 (0.0, capsule_a_radius, 0.0)
          } else {
            capsule_a_radius * a_to_b.normalize()
          };
          let axis      = nearest_b - nearest_a;
          let distance  = axis.magnitude() * if !a_to_b.is_zero() {
            a_to_b.dot (axis).signum()
          } else {
            -1.0
          };
          let half_axis = 0.5 * axis;
          let normal    = if axis.is_zero() {
            math::Unit3::normalize (cgmath::vec3 (0.0, signum_y, -1.0))
          } else {
            math::Unit3::unchecked (-axis / distance)
          };
          let midpoint  = nearest_a + half_axis;
          Proximity { distance, half_axis, midpoint, normal }
        }
        (false, true) => {
          // nearest/deepest point is on an edge parallel to the Y axis
          let signum_x  = (shaft_max_a.x - position_b.x).signum();
          let nearest_b = cgmath::Point3::new (
            position_b.x + signum_x * cuboid_b_half_extent_x,
            shaft_max_a.y,
            cuboid_min_b.z
          );
          let a_to_b    = nearest_b - shaft_max_a;
          let nearest_a = shaft_max_a + if a_to_b.is_zero() {
            -signum_x * cgmath::vec3 (capsule_a_radius, 0.0, 0.0)
          } else {
            capsule_a_radius * a_to_b.normalize()
          };
          let axis      = nearest_b - nearest_a;
          let distance  = axis.magnitude() * if !a_to_b.is_zero() {
            a_to_b.dot (axis).signum()
          } else {
            -1.0
          };
          let half_axis = 0.5 * axis;
          let normal    = if axis.is_zero() {
            math::Unit3::normalize (cgmath::vec3 (signum_x, 0.0, -1.0))
          } else {
            math::Unit3::unchecked (-axis / distance)
          };
          let midpoint  = nearest_a + half_axis;
          Proximity { distance, half_axis, midpoint, normal }
        }
        (true, true) => {
          // nearest/deepest point is in the -Z direction
          let nearest_a : cgmath::Point3 <f64> =
            shaft_max_a + cgmath::vec3 (0.0, 0.0, capsule_a_radius);
          let nearest_b : cgmath::Point3 <f64> =
            [position_a.x, position_a.y, cuboid_min_b.z].into();
          let axis      = nearest_b - nearest_a;
          let normal    = math::Unit3::axis_z().invert();
          let half_axis = 0.5 * axis;
          let distance  = nearest_b.z - nearest_a.z;
          let midpoint  = nearest_a + half_axis;
          Proximity { distance, half_axis, midpoint, normal }
        }
      }
    } else {
      // capsule cylinder shaft overlaps on z axis
      let max_z = f64::min (shaft_max_a.z, cuboid_max_b.z);
      let min_z = f64::max (shaft_min_a.z, cuboid_min_b.z);
      debug_assert_ne!(min_z, max_z);
      let mid_z = 0.5 * (max_z + min_z);
      match (overlap_x, overlap_y) {
        (false, false) => {
          // nearest point is on a Z parallel edge of the cuboid
          let b_to_a_signum_x = (position_a.x - position_b.x).signum();
          let b_to_a_signum_y = (position_a.y - position_b.y).signum();
          let b_to_a_unit_sigvec =
            cgmath::vec3 (b_to_a_signum_x, b_to_a_signum_y, 0.0).normalize();
          let nearest_b = cgmath::Point3::from ([
            position_b.x + b_to_a_signum_x * cuboid_b_half_extent_x,
            position_b.y + b_to_a_signum_y * cuboid_b_half_extent_y,
            mid_z
          ]);
          let a_to_b    = cgmath::vec3 (
            nearest_b.x - position_a.x,
            nearest_b.y - position_a.y,
            0.0
          );
          let nearest_a = cgmath::Point3::new (
            position_a.x, position_a.y, mid_z
          ) + capsule_a_radius * if !a_to_b.is_zero() {
            a_to_b.normalize()
          } else {
            -b_to_a_unit_sigvec
          };
          let axis      = nearest_b - nearest_a;
          let distance  = axis.magnitude() * if !a_to_b.is_zero() {
            a_to_b.dot (axis).signum()
          } else {
            -1.0
          };
          let half_axis = 0.5 * axis;
          let normal    = math::Unit3::unchecked (if axis.is_zero() {
            b_to_a_unit_sigvec
          } else {
            -axis / distance
          });
          let midpoint  = nearest_a + half_axis;
          Proximity { distance, half_axis, midpoint, normal }
        }
        (true, false)  => {
          // nearest point is on a Z/X parallel face of the cuboid
          let b_to_a_signum_y    = (position_a.y - position_b.y).signum();
          let b_to_a_unit_sigvec = cgmath::vec3 (0.0, b_to_a_signum_y, 0.0);
          let nearest_b = cgmath::Point3::from ([
            position_a.x,
            position_b.y + b_to_a_signum_y * cuboid_b_half_extent_y,
            mid_z
          ]);
          let a_to_b    = cgmath::vec3 (0.0, nearest_b.y - position_a.y, 0.0);
          let nearest_a = cgmath::Point3::new (
            position_a.x, position_a.y, mid_z
          ) + cgmath::vec3 (0.0, -b_to_a_signum_y * capsule_a_radius, 0.0);
          let axis      = nearest_b - nearest_a;
          let distance  = axis.magnitude() * if !a_to_b.is_zero() {
            a_to_b.dot (axis).signum()
          } else {
            -1.0
          };
          let half_axis = 0.5 * axis;
          let normal    = math::Unit3::unchecked (if axis.is_zero() {
            b_to_a_unit_sigvec
          } else {
            -axis / distance
          });
          let midpoint  = nearest_a + half_axis;
          Proximity { distance, half_axis, midpoint, normal }
        }
        (false, true)  => {
          // nearest point is on a Z/Y parallel face of the cuboid
          let b_to_a_signum_x    = (position_a.x - position_b.x).signum();
          let b_to_a_unit_sigvec = cgmath::vec3 (b_to_a_signum_x, 0.0, 0.0);
          let nearest_b = cgmath::Point3::from ([
            position_b.x + b_to_a_signum_x * cuboid_b_half_extent_x,
            position_a.y,
            mid_z
          ]);
          let a_to_b    = cgmath::vec3 (nearest_b.x - position_a.x, 0.0, 0.0);
          let nearest_a = cgmath::Point3::new (
            position_a.x, position_a.y, mid_z
          ) + cgmath::vec3 (-b_to_a_signum_x * capsule_a_radius, 0.0, 0.0);
          let axis      = nearest_b - nearest_a;
          let distance  = axis.magnitude() * if !a_to_b.is_zero() {
            a_to_b.dot (axis).signum()
          } else {
            -1.0
          };
          let half_axis = 0.5 * axis;
          let normal    = math::Unit3::unchecked (if axis.is_zero() {
            b_to_a_unit_sigvec
          } else {
            -axis / distance
          });
          let midpoint  = nearest_a + half_axis;
          Proximity { distance, half_axis, midpoint, normal }
        }
        (true, true) => {
          // cylinder shaft inside X/Y bounds
          if position_a == position_b {
            let depth_x = cuboid_b_half_extent_x + capsule_a_radius;
            let depth_y = cuboid_b_half_extent_y + capsule_a_radius;
            let depth_z = cuboid_b_half_extent_z + capsule_a_half_height +
              capsule_a_radius;
            if depth_x <= depth_y && depth_x <= depth_z {
              let distance  = -depth_x;
              let half_axis = [0.5 * depth_x, 0.0, 0.0].into();
              let normal    = math::Unit3::axis_x();
              let midpoint  = position_a + half_axis -
                cgmath::vec3 (capsule_a_radius, 0.0, 0.0);
              Proximity { distance, half_axis, midpoint, normal }
            } else if depth_y <= depth_x && depth_y <= depth_z {
              let distance = -depth_y;
              let half_axis = [0.0, 0.5 * depth_y, 0.0].into();
              let normal    = math::Unit3::axis_y();
              let midpoint  = position_a + half_axis -
                cgmath::vec3 (0.0, capsule_a_radius, 0.0);
              Proximity { distance, half_axis, midpoint, normal }
            } else {
              debug_assert!(depth_z < depth_x);
              debug_assert!(depth_z < depth_y);
              let distance = -depth_z;
              let half_axis = [0.0, 0.0, 0.5 * depth_z].into();
              let normal    = math::Unit3::axis_z();
              let midpoint  = position_a + half_axis - cgmath::vec3 (
                0.0, 0.0, capsule_a_radius + capsule_a_half_height);
              Proximity { distance, half_axis, midpoint, normal }
            }
          } else {
            // NB: we want the following to map zero values to positive 1.0
            // so we use f64::signum instead of math::sigvec
            let b_to_a_sigvec = (position_a - position_b).map (f64::signum);
            let corner_b      = position_b +
              cuboid_b.max().to_vec().mul_element_wise (b_to_a_sigvec);
            let depth_x       = (corner_b.x - position_a.x).abs() +
              capsule_a_radius;
            let depth_y       = (corner_b.y - position_a.y).abs() +
              capsule_a_radius;
            let depth_z       = (corner_b.z - position_a.z).abs() +
              capsule_a_radius + capsule_a_half_height -
              if position_a.z.abs() > corner_b.z.abs() {
                2.0 * (position_a.z.abs() - corner_b.z.abs())
              } else {
                0.0
              };

            if depth_x <= depth_y && depth_x <= depth_z {
              let distance  = -depth_x;
              let half_axis = [0.5 * b_to_a_sigvec.x * depth_x, 0.0, 0.0].into();
              let normal    = math::Unit3::unchecked (
                [b_to_a_sigvec.x, 0.0, 0.0].into());
              let midpoint  = position_a + half_axis +
                cgmath::vec3 (-b_to_a_sigvec.x * capsule_a_radius, 0.0, 0.0);
              Proximity { distance, half_axis, midpoint, normal }
            } else if depth_y <= depth_x && depth_y <= depth_z {
              let distance  = -depth_y;
              let half_axis = [0.0, 0.5 * b_to_a_sigvec.y * depth_y, 0.0].into();
              let normal    = math::Unit3::unchecked (
                [0.0, b_to_a_sigvec.y, 0.0].into());
              let midpoint  = position_a + half_axis +
                cgmath::vec3 (0.0, -b_to_a_sigvec.y * capsule_a_radius, 0.0);
              Proximity { distance, half_axis, midpoint, normal }
            } else {
              debug_assert!(depth_z < depth_x);
              debug_assert!(depth_z < depth_y);
              let nearest_a = cgmath::Point3::new (
                position_a.x,
                position_a.y,
                position_a.z - b_to_a_sigvec.z *
                  (capsule_a_radius + capsule_a_half_height)
              );
              let nearest_b = cgmath::Point3::new (
                position_a.x,
                position_a.y,
                position_b.z + b_to_a_sigvec.z * cuboid_b_half_extent_z
              );
              let axis      = nearest_b - nearest_a;
              debug_assert_eq!(axis.x, 0.0);
              debug_assert_eq!(axis.y, 0.0);
              let distance  = -axis.z.abs();
              let half_axis = 0.5 * axis;
              let normal    = math::Unit3::unchecked (
                [0.0, 0.0, b_to_a_sigvec.z].into());
              let midpoint  = nearest_a + half_axis;
              Proximity { distance, half_axis, midpoint, normal }
            }
          }
        }
      }
    }
  } // end query_capsule_cuboid

  pub fn query_capsule_orthant (
    position_a : &cgmath::Point3 <f64>, capsule_a : &shape::Capsule <f64>,
    position_b : &cgmath::Point3 <f64>, orthant_b : &shape::Orthant
  ) -> Self {
    use cgmath::{Array, ElementWise, EuclideanSpace, InnerSpace};
    let radius_a    = *capsule_a.radius;
    let half_height = *capsule_a.half_height;
    let normal      = orthant_b.normal_axis.to_vec::<f64>();
    let normal_abs  = normal.map (f64::abs);
    let surface_sigvec =
      cgmath::vec3 (1.0, 1.0, 1.0).sub_element_wise (normal_abs);
    let nearest_a   = position_a -
      cgmath::vec3 (radius_a, radius_a, radius_a + half_height)
        .mul_element_wise (normal);
    let nearest_b   = cgmath::Point3::from_vec (
      position_a.to_vec().mul_element_wise (surface_sigvec) +
      position_b.to_vec().mul_element_wise (normal_abs));
    let axis        = nearest_b - nearest_a;
    debug_assert_eq!(axis.sum().abs(), axis.magnitude());
    let distance    = -normal.dot (axis).signum() * axis.sum().abs();
    let half_axis   = 0.5 * axis;
    let midpoint    = nearest_a + half_axis;
    let normal      = math::Unit3::unchecked (normal);
    Proximity { distance, half_axis, midpoint, normal }
  } // end query_capsule_orthant
} // end impl Distance

#[cfg(test)]
mod tests {
  use {std, test, rand};
  use {rs_utils};
  use std::f64::consts::*;
  use super::*;
  //
  //  test_distance_query
  //
  #[test]
  fn test_distance_query() {
    use geometry::shape;
    use cgmath::InnerSpace;

    //
    //  capsule v. capsule
    //

    // intersection: positions of capsules A and B are identical-- half_axis and
    // normal parallel to +X
    let a = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound (shape::Bounded::from (
        shape::Capsule::noisy (2.0, 3.0)).into()),
      material: component::MATERIAL_STONE
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound (shape::Bounded::from (
        shape::Capsule::noisy (1.0, 2.0)).into()),
      material: component::MATERIAL_STONE
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -3.0,
        half_axis: [ 1.5, 0.0, 0.0].into(),
        normal:    math::Unit3::axis_x(),
        midpoint:  [-0.5, 0.0, 0.0].into()
      }
    );

    // intersection: cylinder shafts of A and B overlap-- half_axis and normal
    // parallel to +X
    let a = object::Static {
      position: component::Position ([0.0, 0.0, 1.0].into()),  .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, -1.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -3.0,
        half_axis: [ 1.5, 0.0,  0.0].into(),
        normal:    math::Unit3::axis_x(),
        midpoint:  [-0.5, 0.0, -0.5].into()
      }
    );

    // intersection: cylinder shaft max of A coincides with cylinder shaft min
    // of B
    let a = object::Static {
      position: component::Position ([0.0, 0.0, -3.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 2.0].into()),  .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -3.0,
        half_axis: [ 1.5, 0.0, 0.0].into(),
        normal:    math::Unit3::axis_x(),
        midpoint:  [-0.5, 0.0, 0.0].into()
      }
    );

    // intersection: cylinder shaft min of A coincides with cylinder shaft max of B
    let a = object::Static {
      position: component::Position ([0.0, 0.0, 3.0].into()),  .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, -2.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -3.0,
        half_axis: [ 1.5, 0.0, 0.0].into(),
        normal:    math::Unit3::axis_x(),
        midpoint:  [-0.5, 0.0, 0.0].into()
      }
    );

    // intersection: cylinder shaft of A adjacent to cylinder shaft of B
    let a = object::Static {
      position: component::Position ([-1.0, 0.0, 1.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, -1.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -2.0,
        half_axis: [-1.0, 0.0,  0.0].into(),
        normal:    math::Unit3::axis_x().invert(),
        midpoint:  [ 0.0, 0.0, -0.5].into()
      }
    );

    // intersection: cylinder shaft of B above cylinder shaft of A
    let a = object::Static {
      position: component::Position ([0.0, 0.0, -3.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0,  3.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -2.0,
        half_axis: [ 0.0, 0.0, -1.0].into(),
        normal:    math::Unit3::axis_z().invert(),
        midpoint:  [ 0.0, 0.0,  1.0].into()
      }
    );

    // intersection: cylinder shaft of A above cylinder shaft of B
    let a = object::Static {
      position: component::Position ([0.0, 0.0,  3.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, -3.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -2.0,
        half_axis: [ 0.0, 0.0,  1.0].into(),
        normal:    math::Unit3::axis_z(),
        midpoint:  [ 0.0, 0.0, -1.0].into()
      }
    );

    // intersection: cylinder shafts of A and B are adjacent
    let a = object::Static {
      position: component::Position ([ 1.0, 0.0, 0.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([-1.0, 0.0, 0.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -1.0,
        half_axis: [ 0.5, 0.0, 0.0].into(),
        normal:    math::Unit3::axis_x(),
        midpoint:  [-0.5, 0.0, 0.0].into()
      }
    );

    // disjoint: capsule A above and to the right (+X) of capsule B
    let a = object::Static {
      position: component::Position ([ 2.0, 0.0,  5.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([-1.0, 0.0, -3.0].into()), .. b
    };

    let proximity = Proximity::query (&a, &b);
    let distance_manual = {
      let distance  = cgmath::vec3 (3.0, 0.0, 3.0).magnitude() - 3.0;
      let half_axis = 0.5 * distance * cgmath::vec3 (-1.0, 0.0, -1.0).normalize();
      let normal    = math::Unit3::normalize (cgmath::vec3 (1.0, 0.0, 1.0));
      let midpoint  = cgmath::Point3::from ([-1.0, 0.0, -1.0]) +
        *normal - half_axis;
      Proximity { distance, half_axis, midpoint, normal }
    };
    assert_eq!(proximity.distance, distance_manual.distance);
    assert_relative_eq!(proximity.half_axis, distance_manual.half_axis);
    assert_relative_eq!(*proximity.normal,   *distance_manual.normal);
    assert_relative_eq!(proximity.midpoint,  distance_manual.midpoint,
      epsilon = 2.0 * std::f64::EPSILON);

    // disjoint: cylinder shafts of A and B are adjacent
    let a = object::Static {
      position: component::Position ([ 3.0, 0.0, 3.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([-1.0, 0.0, 0.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [-0.5, 0.0, 0.0].into(),
        normal:    math::Unit3::axis_x(),
        midpoint:  [ 0.5, 0.0, 1.0].into()
      }
    );

    //
    //  cuboid v. cuboid
    //

    // disjoint: nearest points are corners
    let a = object::Static {
      position: component::Position ([-5.0, -5.0, -5.0].into()),
      bound:    component::Bound (shape::Bounded::from (
        shape::Cuboid::noisy ([1.0, 2.0, 3.0])).into()),
      .. a
    };
    let b = object::Static {
      position: component::Position ([5.0, 5.0, 5.0].into()),
      bound:    component::Bound (shape::Bounded::from (
        shape::Cuboid::noisy ([3.0, 1.0, 2.0])).into()),
      .. b
    };
    let proximity = Proximity::query (&a, &b);
    assert_eq!(proximity.distance, f64::sqrt (110.0));
    assert_eq!(proximity.half_axis, [3.0, 3.5, 2.5].into());
    assert_ulps_eq!(*proximity.normal, -cgmath::vec3 (3.0, 3.5, 2.5).normalize());
    assert_eq!(proximity.midpoint,  [-1.0, 0.5, 0.5].into());

    // disjoint: overlap on X axis only
    let a = object::Static {
      position: component::Position ([-5.0, 0.0, -5.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([5.0, 0.0, 5.0].into()), .. b
    };
    let proximity = Proximity::query (&a, &b);
    assert_eq!(proximity.distance, f64::sqrt (61.0));
    assert_eq!(proximity.half_axis, [3.0, 0.0, 2.5].into());
    assert_ulps_eq!(*proximity.normal, -cgmath::vec3 (3.0, 0.0, 2.5).normalize());
    assert_eq!(proximity.midpoint,  [-1.0, 0.0, 0.5].into());

    //
    //  capsule v. cuboid
    //

    // disjoint: capsule nearest cuboid max corner
    let a = object::Static {
      position: component::Position ([5.0, 5.0, 3.0].into()),
      bound:    component::Bound (shape::Bounded::from (
        shape::Capsule::noisy (SQRT_2, 2.0)).into()),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound (shape::Bounded::from (
        shape::Cuboid::noisy ([1.0, 1.0, 1.0])).into()),
      .. b
    };
    let proximity = Proximity::query (&a, &b);
    assert_ulps_eq!(proximity.distance, 3.0 * SQRT_2);
    assert_eq!(proximity.half_axis, [-1.5, -1.5, 0.0].into());
    assert_ulps_eq!(*proximity.normal, cgmath::vec3 (1.0, 1.0, 0.0).normalize());
    assert_eq!(proximity.midpoint,  [ 2.5, 2.5, 1.0].into());

    // disjoint: capsule nearest cuboid min corner
    let a = object::Static {
      position: component::Position ([-5.0, -5.0, -3.0].into()),
      bound:    component::Bound (shape::Bounded::from (
        shape::Capsule::noisy (SQRT_2, 2.0)).into()),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound (shape::Bounded::from (
        shape::Cuboid::noisy ([1.0, 1.0, 1.0])).into()),
      .. b
    };
    let proximity = Proximity::query (&a, &b);
    assert_ulps_eq!(proximity.distance, 3.0 * SQRT_2);
    assert_eq!(proximity.half_axis, [ 1.5,  1.5, 0.0].into());
    assert_ulps_eq!(*proximity.normal, cgmath::vec3 (-1.0, -1.0, 0.0).normalize());
    assert_eq!(proximity.midpoint,  [-2.5,-2.5,-1.0].into());

    // disjoint: capsule nearest cuboid max Y edge
    let a = object::Static {
      position: component::Position ([4.0, 0.5, 4.0].into()),
      bound:    component::Bound (shape::Bounded::from (
        shape::Capsule::noisy (1.0, 2.0)).into()),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound (shape::Bounded::from (
        shape::Cuboid::noisy ([2.0, 2.0, 2.0])).into()),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [-0.5, 0.0, 0.0].into(),
        normal:    math::Unit3::axis_x(),
        midpoint:  [ 2.5, 0.5, 2.0].into()
      }
    );

    // disjoint: capsule nearest cuboid min Y edge
    let a = object::Static {
      position: component::Position ([4.0, 0.5, -4.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [-0.5, 0.0,  0.0].into(),
        normal:    math::Unit3::axis_x(),
        midpoint:  [ 2.5, 0.5, -2.0].into()
      }
    );

    // disjoint: capsule nearest cuboid max X edge
    let a = object::Static {
      position: component::Position ([0.5, 4.0, 4.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [ 0.0, -0.5, 0.0].into(),
        normal:    math::Unit3::axis_y(),
        midpoint:  [ 0.5,  2.5, 2.0].into()
      }
    );

    // disjoint: capsule nearest cuboid min X edge
    let a = object::Static {
      position: component::Position ([0.5, 4.0, -4.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [ 0.0, -0.5,  0.0].into(),
        normal:    math::Unit3::axis_y(),
        midpoint:  [ 0.5,  2.5, -2.0].into()
      }
    );

    // disjoint: capsule A above cuboid B
    let a = object::Static {
      position: component::Position ([0.5, 0.5, 6.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [ 0.0, 0.0, -0.5].into(),
        normal:    math::Unit3::axis_z(),
        midpoint:  [ 0.5, 0.5,  2.5].into()
      }
    );

    // disjoint: capsule A below cuboid B
    let a = object::Static {
      position: component::Position ([0.5, 0.5, -6.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [ 0.0, 0.0,  0.5].into(),
        normal:    math::Unit3::axis_z().invert(),
        midpoint:  [ 0.5, 0.5, -2.5].into()
      }
    );

    // disjoint: capsule A cylinder nearest max Z edge of cuboid B
    let a = object::Static {
      position: component::Position ([3.0, 3.0, 2.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()), .. b
    };
    let proximity = Proximity::query (&a, &b);
    assert_eq!(proximity.distance, SQRT_2 - 1.0);
    assert_eq!(proximity.half_axis, [
      -0.5 * f64::sqrt (0.5 * (SQRT_2 - 1.0).powf (2.0)),
      -0.5 * f64::sqrt (0.5 * (SQRT_2 - 1.0).powf (2.0)),
       0.0
    ].into());
    assert_ulps_eq!(*proximity.normal,
      [FRAC_1_SQRT_2, FRAC_1_SQRT_2, 0.0].into());
    assert_eq!(proximity.midpoint, [
      2.0 + 0.5 * f64::sqrt (0.5 * (SQRT_2 - 1.0).powf (2.0)),
      2.0 + 0.5 * f64::sqrt (0.5 * (SQRT_2 - 1.0).powf (2.0)),
      1.0
    ].into());

    // disjoint: capsule A cylinder nearest min Z edge of cuboid B
    let a = object::Static {
      position: component::Position ([-3.0, -3.0, -2.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()), .. b
    };
    let proximity = Proximity::query (&a, &b);
    assert_eq!(proximity.distance, SQRT_2 - 1.0);
    assert_eq!(proximity.half_axis, [
       0.5 * f64::sqrt (0.5 * (SQRT_2 - 1.0).powf (2.0)),
       0.5 * f64::sqrt (0.5 * (SQRT_2 - 1.0).powf (2.0)),
       0.0
    ].into());
    assert_ulps_eq!(*proximity.normal,
      [-FRAC_1_SQRT_2, -FRAC_1_SQRT_2, 0.0].into());
    assert_eq!(proximity.midpoint, [
      -2.0 - 0.5 * f64::sqrt (0.5 * (SQRT_2 - 1.0).powf (2.0)),
      -2.0 - 0.5 * f64::sqrt (0.5 * (SQRT_2 - 1.0).powf (2.0)),
      -1.0
    ].into());

    // disjoint: capsule A cylinder nearest max X/Z face of cuboid B
    let a = object::Static {
      position: component::Position ([1.0, 4.0, 2.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [ 0.0, -0.5,  0.0].into(),
        normal:    math::Unit3::axis_y(),
        midpoint:  [ 1.0,  2.5,  1.0].into()
      }
    );

    // disjoint: capsule A cylinder nearest min X/Z face of cuboid B
    let a = object::Static {
      position: component::Position ([-1.0, -4.0, -2.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),    .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [ 0.0,  0.5,  0.0].into(),
        normal:    math::Unit3::axis_y().invert(),
        midpoint:  [-1.0, -2.5, -1.0].into()
      }
    );

    // disjoint: capsule A cylinder nearest max Y/Z face of cuboid B
    let a = object::Static {
      position: component::Position ([4.0, 1.0, 2.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [-0.5, 0.0, 0.0].into(),
        normal:    math::Unit3::axis_x(),
        midpoint:  [ 2.5, 1.0, 1.0].into()
      }
    );

    // disjoint: capsule A cylinder nearest min Y/Z face of cuboid B
    let a = object::Static {
      position: component::Position ([-4.0, -1.0, -2.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),    .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [ 0.5,  0.0,  0.0].into(),
        normal:    math::Unit3::axis_x().invert(),
        midpoint:  [-2.5, -1.0, -1.0].into()
      }
    );

    // intersect: capsule cylinder shaft min equal to cuboid max corner
    let a = object::Static {
      position: component::Position ([2.0, 2.0, 4.0].into()),
      bound:    component::Bound (shape::Bounded::from (
        shape::Capsule::noisy (1.0, 2.0)).into()),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound (shape::Bounded::from (
        shape::Cuboid::noisy ([2.0, 2.0, 2.0])).into()),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -1.0,
        half_axis: [ 0.5, 0.0, 0.0].into(),
        normal:    math::Unit3::axis_x(),
        midpoint:  [ 1.5, 2.0, 2.0].into()
      }
    );

    // intersect: capsule cylinder shaft max equal to cuboid min corner
    let a = object::Static {
      position: component::Position ([-2.0, -2.0, -4.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),  .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -1.0,
        half_axis: [-0.5,  0.0, 0.0].into(),
        normal:    math::Unit3::axis_x().invert(),
        midpoint:  [-1.5, -2.0,-2.0].into()
      }
    );

    // intersect: capsule cylinder shaft min equal to cuboid max Y edge
    let a = object::Static {
      position: component::Position ([2.0, 1.0, 4.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -1.0,
        half_axis: [ 0.5, 0.0, 0.0].into(),
        normal:    math::Unit3::axis_x(),
        midpoint:  [ 1.5, 1.0, 2.0].into()
      }
    );

    // intersect: capsule cylinder shaft max equal to cuboid min Y edge
    let a = object::Static {
      position: component::Position ([-2.0, 1.0, -4.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),  .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -1.0,
        half_axis: [-0.5, 0.0,  0.0].into(),
        normal:    math::Unit3::axis_x().invert(),
        midpoint:  [-1.5, 1.0, -2.0].into()
      }
    );

    // intersect: capsule cylinder shaft min equal to cuboid max X edge
    let a = object::Static {
      position: component::Position ([1.0, 2.0, 4.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -1.0,
        half_axis: [ 0.0, 0.5, 0.0].into(),
        normal:    math::Unit3::axis_y(),
        midpoint:  [ 1.0, 1.5, 2.0].into()
      }
    );

    // intersect: capsule cylinder shaft max equal to cuboid min X edge
    let a = object::Static {
      position: component::Position ([ 1.0, -2.0, -4.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),  .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -1.0,
        half_axis: [0.0, -0.5,  0.0].into(),
        normal:    math::Unit3::axis_y().invert(),
        midpoint:  [1.0, -1.5, -2.0].into()
      }
    );

    // intersect: capsule cylinder shaft min equal to cuboid max Z face
    let a = object::Static {
      position: component::Position ([1.0, 1.0, 4.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -1.0,
        half_axis: [ 0.0, 0.0, 0.5].into(),
        normal:    math::Unit3::axis_z(),
        midpoint:  [ 1.0, 1.0, 1.5].into()
      }
    );

    // intersect: capsule cylinder shaft max equal to cuboid min Z face
    let a = object::Static {
      position: component::Position ([-1.0, -1.0, -4.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -1.0,
        half_axis: [ 0.0,  0.0, -0.5].into(),
        normal:    math::Unit3::axis_z().invert(),
        midpoint:  [-1.0, -1.0, -1.5].into()
      }
    );

    // intersect: capsule and cuboid position coincide
    let a = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -3.0,
        half_axis: [ 1.5,  0.0,  0.0].into(),
        normal:    math::Unit3::axis_x(),
        midpoint:  [ 0.5,  0.0,  0.0].into()
      }
    );

    // intersect: capsule nearest positive X face of cuboid
    let a = object::Static {
      position: component::Position ([1.0, 0.0, 0.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -2.0,
        half_axis: [ 1.0,  0.0,  0.0].into(),
        normal:    math::Unit3::axis_x(),
        midpoint:  [ 1.0,  0.0,  0.0].into()
      }
    );

    // intersect: capsule nearest negative X face of cuboid
    let a = object::Static {
      position: component::Position ([-1.0, 0.0, 0.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -2.0,
        half_axis: [-1.0,  0.0,  0.0].into(),
        normal:    math::Unit3::axis_x().invert(),
        midpoint:  [-1.0,  0.0,  0.0].into()
      }
    );

    // intersect: capsule nearest positive Y face of cuboid
    let a = object::Static {
      position: component::Position ([0.0, 1.0, 0.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -2.0,
        half_axis: [0.0, 1.0, 0.0].into(),
        normal:    math::Unit3::axis_y(),
        midpoint:  [0.0, 1.0, 0.0].into()
      }
    );

    // intersect: capsule nearest negative Y face of cuboid
    let a = object::Static {
      position: component::Position ([0.0, -1.0, 0.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -2.0,
        half_axis: [0.0, -1.0, 0.0].into(),
        normal:    math::Unit3::axis_y().invert(),
        midpoint:  [0.0, -1.0, 0.0].into()
      }
    );

    // intersect: capsule nearest positive Z face of cuboid
    let a = object::Static {
      position: component::Position ([0.0, 0.0, 3.0].into()),
      bound:    component::Bound (shape::Bounded::from (
        shape::Capsule::noisy (1.0, 2.0)).into()),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound (shape::Bounded::from (
        shape::Cuboid::noisy ([4.0, 4.0, 2.0])).into()),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -2.0,
        half_axis: [0.0, 0.0, 1.0].into(),
        normal:    math::Unit3::axis_z(),
        midpoint:  [0.0, 0.0, 1.0].into()
      }
    );

    // intersect: capsule nearest positive Z face of cuboid
    let a = object::Static {
      position: component::Position ([0.0, 0.0, 1.0].into()),
      bound:    component::Bound (shape::Bounded::from (
        shape::Capsule::noisy (1.0, 2.0)).into()),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound (shape::Bounded::from (
        shape::Cuboid::noisy ([4.0, 4.0, 2.0])).into()),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -4.0,
        half_axis: [0.0, 0.0, 2.0].into(),
        normal:    math::Unit3::axis_z(),
        midpoint:  [0.0, 0.0, 0.0].into()
      }
    );

    // intersect: capsule nearest negative Z face of cuboid
    let a = object::Static {
      position: component::Position ([0.0, 0.0, -3.0].into()),
      bound:    component::Bound (shape::Bounded::from (
        shape::Capsule::noisy (1.0, 2.0)).into()),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound (shape::Bounded::from (
        shape::Cuboid::noisy ([4.0, 4.0, 2.0])).into()),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -2.0,
        half_axis: [0.0, 0.0, -1.0].into(),
        normal:    math::Unit3::axis_z().invert(),
        midpoint:  [0.0, 0.0, -1.0].into()
      }
    );

    // intersect: capsule nearest negative Z face of cuboid
    let a = object::Static {
      position: component::Position ([0.0, 0.0, -1.0].into()),
      bound:    component::Bound (shape::Bounded::from (
        shape::Capsule::noisy (1.0, 2.0)).into()),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound (shape::Bounded::from (
        shape::Cuboid::noisy ([4.0, 4.0, 2.0])).into()),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -4.0,
        half_axis: [0.0, 0.0, -2.0].into(),
        normal:    math::Unit3::axis_z().invert(),
        midpoint:  [0.0, 0.0,  0.0].into()
      }
    );

    //
    //  capsule v. orthant
    //
    use math::SignedAxis3;

    // disjoint: capsule outside of orthant +Z
    let a = object::Static {
      position: component::Position ([0.0, 0.0, 4.0].into()),
      bound:    component::Bound (shape::Bounded::from (
        shape::Capsule::noisy (1.0, 2.0)).into()),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound (shape::Unbounded::from (
        shape::Orthant { normal_axis: SignedAxis3::PosZ }).into()),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [0.0, 0.0, -0.5].into(),
        normal:    math::Unit3::axis_z(),
        midpoint:  [0.0, 0.0,  0.5].into()
      }
    );

    // disjoint: capsule outside of orthant -Z
    let a = object::Static {
      position: component::Position ([0.0, 0.0, -4.0].into()),
      bound:    component::Bound (shape::Bounded::from (
        shape::Capsule::noisy (1.0, 2.0)).into()),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound (shape::Unbounded::from (
        shape::Orthant { normal_axis: SignedAxis3::NegZ }).into()),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [0.0, 0.0,  0.5].into(),
        normal:    math::Unit3::axis_z().invert(),
        midpoint:  [0.0, 0.0, -0.5].into()
      }
    );

    // disjoint: capsule outside orthant +Y
    let a = object::Static {
      position: component::Position ([0.0, 2.0, 0.0].into()),
      bound:    component::Bound (shape::Bounded::from (
        shape::Capsule::noisy (1.0, 2.0)).into()),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound (shape::Unbounded::from (
        shape::Orthant { normal_axis: SignedAxis3::PosY }).into()),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [0.0, -0.5, 0.0].into(),
        normal:    math::Unit3::axis_y(),
        midpoint:  [0.0,  0.5, 0.0].into()
      }
    );

    // disjoint: capsule outside orthant -Y
    let a = object::Static {
      position: component::Position ([0.0, -2.0, 0.0].into()),
      bound:    component::Bound (shape::Bounded::from (
        shape::Capsule::noisy (1.0, 2.0)).into()),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound (shape::Unbounded::from (
        shape::Orthant { normal_axis: SignedAxis3::NegY }).into()),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [0.0,  0.5, 0.0].into(),
        normal:    math::Unit3::axis_y().invert(),
        midpoint:  [0.0, -0.5, 0.0].into()
      }
    );

    // disjoint: capsule outside orthant +X
    let a = object::Static {
      position: component::Position ([2.0, 0.0, 0.0].into()),
      bound:    component::Bound (shape::Bounded::from (
        shape::Capsule::noisy (1.0, 2.0)).into()),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound (shape::Unbounded::from (
        shape::Orthant { normal_axis: SignedAxis3::PosX }).into()),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [-0.5, 0.0, 0.0].into(),
        normal:    math::Unit3::axis_x(),
        midpoint:  [ 0.5, 0.0, 0.0].into()
      }
    );

    // disjoint: outside orthant -X
    let a = object::Static {
      position: component::Position ([-2.0, 0.0, 0.0].into()),
      bound:    component::Bound (shape::Bounded::from (
        shape::Capsule::noisy (1.0, 2.0)).into()),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound (shape::Unbounded::from (
        shape::Orthant { normal_axis: SignedAxis3::NegX }).into()),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [ 0.5, 0.0, 0.0].into(),
        normal:    math::Unit3::axis_x().invert(),
        midpoint:  [-0.5, 0.0, 0.0].into()
      }
    );

    // intersect: capsule and orthant position coincide; normal +Z
    let a = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound (shape::Bounded::from (
        shape::Capsule::noisy (1.0, 2.0)).into()),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound (shape::Unbounded::from (
        shape::Orthant { normal_axis: SignedAxis3::PosZ }).into()),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -3.0,
        half_axis: [0.0, 0.0,  1.5].into(),
        normal:    math::Unit3::axis_z(),
        midpoint:  [0.0, 0.0, -1.5].into()
      }
    );

    // intersect: capsule and orthant position coincide; normal -Z
    let a = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound (shape::Bounded::from (
        shape::Capsule::noisy (1.0, 2.0)).into()),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound (shape::Unbounded::from (
        shape::Orthant { normal_axis: SignedAxis3::NegZ }).into()),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -3.0,
        half_axis: [0.0, 0.0, -1.5].into(),
        normal:    math::Unit3::axis_z().invert(),
        midpoint:  [0.0, 0.0,  1.5].into()
      }
    );

    // intersect: capsule and orthant position coincide; normal +Y
    let a = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound (shape::Bounded::from (
        shape::Capsule::noisy (1.0, 2.0)).into()),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound (shape::Unbounded::from (
        shape::Orthant { normal_axis: SignedAxis3::PosY }).into()),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -1.0,
        half_axis: [0.0,  0.5, 0.0].into(),
        normal:    math::Unit3::axis_y(),
        midpoint:  [0.0, -0.5, 0.0].into()
      }
    );

    // intersect: capsule and orthant position coincide; normal -Y
    let a = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound (shape::Bounded::from (
        shape::Capsule::noisy (1.0, 2.0)).into()),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound (shape::Unbounded::from (
        shape::Orthant { normal_axis: SignedAxis3::NegY }).into()),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -1.0,
        half_axis: [0.0, -0.5, 0.0].into(),
        normal:    math::Unit3::axis_y().invert(),
        midpoint:  [0.0,  0.5, 0.0].into()
      }
    );

    // intersect: capsule and orthant position coincide; normal +X
    let a = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound (shape::Bounded::from (
        shape::Capsule::noisy (1.0, 2.0)).into()),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound (shape::Unbounded::from (
        shape::Orthant { normal_axis: SignedAxis3::PosX }).into()),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -1.0,
        half_axis: [ 0.5, 0.0, 0.0].into(),
        normal:    math::Unit3::axis_x(),
        midpoint:  [-0.5, 0.0, 0.0].into()
      }
    );

    // intersect: capsule and orthant position coincide; normal -X
    let a = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound (shape::Bounded::from (
        shape::Capsule::noisy (1.0, 2.0)).into()),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound (shape::Unbounded::from (
        shape::Orthant { normal_axis: SignedAxis3::NegX }).into()),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -1.0,
        half_axis: [-0.5, 0.0, 0.0].into(),
        normal:    math::Unit3::axis_x().invert(),
        midpoint:  [ 0.5, 0.0, 0.0].into()
      }
    );

  } // end test_distance_query

  #[bench]
  /// 180ns
  fn bench_distance_query_cuboid_edges_normalize (b : &mut test::Bencher) {
    let mut rng = rs_utils::numeric::xorshift_rng_unseeded();
    b.iter (||{
      let (position_a, cuboid_a, position_b, cuboid_b) =
        distance_query_cuboid_edge_case (&mut rng);
      let proximity = Proximity::query_cuboid_cuboid (
        &position_a, &cuboid_a, &position_b, &cuboid_b);
      assert!(proximity.distance > 0.0);
    });
  }

  /*
  #[bench]
  /// 180ns
  fn bench_distance_query_cuboid_edges_array (b : &mut test::Bencher) {
    let mut rng = rs_utils::numeric::xorshift_rng_unseeded();
    b.iter (||{
      let (position_a, cuboid_a, position_b, cuboid_b) =
        distance_query_cuboid_edge_case (&mut rng);
      let proximity = Proximity::query_cuboid_cuboid_array (
        &position_a, &cuboid_a, &position_b, &cuboid_b);
      assert!(proximity.distance > 0.0);
    });
  }

  /// 190ns
  #[bench]
  fn bench_distance_query_cuboid_edges_refactor (b : &mut test::Bencher) {
    let mut rng = rs_utils::numeric::xorshift_rng_unseeded();
    b.iter (||{
      let (position_a, cuboid_a, position_b, cuboid_b) =
        distance_query_cuboid_edge_case (&mut rng);
      let proximity = Proximity::query_cuboid_cuboid_refactor (
        &position_a, &cuboid_a, &position_b, &cuboid_b);
      assert!(proximity.distance > 0.0);
    });
  }
  */

  /// Generate an edge case for cuboid/cuboid distance benchmarks
  fn distance_query_cuboid_edge_case <RNG : rand::Rng> (rng : &mut RNG) -> (
    cgmath::Point3 <f64>, shape::Cuboid <f64>,
    cgmath::Point3 <f64>, shape::Cuboid <f64>
  ) {
    // overlapping range
    let base_min_a : f64 = rng.gen_range (-10.0, 10.0);
    let base_max_a : f64 = rng.gen_range (base_min_a + 0.01, 10.01);
    let base_min_b : f64 = rng.gen_range (-10.01, base_max_a);
    let base_max_b : f64 = if base_min_b < base_min_a {
      rng.gen_range (base_min_a + 0.01, 10.02)
    } else {
      rng.gen_range (base_min_b + 0.01, 10.02)
    };
    debug_assert!(base_min_a < base_max_b);
    debug_assert!(base_min_b < base_max_a);
    let base_pos_a = 0.5 * (base_min_a + base_max_a);
    let base_half_extent_a = base_max_a - base_pos_a;
    let base_pos_b = 0.5 * (base_min_b + base_max_b);
    let base_half_extent_b = base_max_b - base_pos_b;
    debug_assert!(base_half_extent_a > 0.0);
    debug_assert!(base_half_extent_b > 0.0);
    // non-overlapping ranges
    let other1_min_a = rng.gen_range (-10.0, 10.0);
    let other1_max_a = rng.gen_range (other1_min_a + 0.01, 10.01);
    let (other1_min_b, other1_max_b) = if rng.gen_bool (0.5) {
      // greater than
      let other1_min_b = rng.gen_range (other1_max_a + 0.01, 10.02);
      let other1_max_b = rng.gen_range (other1_min_b + 0.01, 10.03);
      debug_assert!(other1_max_a < other1_min_b);
      (other1_min_b, other1_max_b)
    } else {
      // less than
      let other1_min_b = rng.gen_range (-10.03, other1_min_a - 0.03);
      let other1_max_b = rng.gen_range (other1_min_b + 0.01, other1_min_a - 0.01);
      debug_assert!(other1_max_b < other1_min_a);
      (other1_min_b, other1_max_b)
    };
    let other1_pos_a = 0.5 * (other1_min_a + other1_max_a);
    let other1_half_extent_a = other1_max_a - other1_pos_a;
    let other1_pos_b = 0.5 * (other1_min_b + other1_max_b);
    let other1_half_extent_b = other1_max_b - other1_pos_b;
    debug_assert!(other1_half_extent_a > 0.0);
    debug_assert!(other1_half_extent_b > 0.0);

    let other2_min_a = rng.gen_range (-10.0, 10.0);
    let other2_max_a = rng.gen_range (other2_min_a + 0.01, 10.01);
    let (other2_min_b, other2_max_b) = if rng.gen_bool (0.5) {
      // greater than
      let other2_min_b = rng.gen_range (other2_max_a + 0.01, 10.02);
      let other2_max_b = rng.gen_range (other2_min_b + 0.01, 10.03);
      debug_assert!(other2_max_a < other2_min_b);
      (other2_min_b, other2_max_b)
    } else {
      // less than
      let other2_min_b = rng.gen_range (-10.03, other2_min_a - 0.03);
      let other2_max_b = rng.gen_range (other2_min_b + 0.01, other2_min_a - 0.01);
      debug_assert!(other2_max_b < other2_min_a);
      (other2_min_b, other2_max_b)
    };
    let other2_pos_a = 0.5 * (other2_min_a + other2_max_a);
    let other2_half_extent_a = other2_max_a - other2_pos_a;
    let other2_pos_b = 0.5 * (other2_min_b + other2_max_b);
    let other2_half_extent_b = other2_max_b - other2_pos_b;
    debug_assert!(other2_half_extent_a > 0.0);
    debug_assert!(other2_half_extent_b > 0.0);

    let axis : usize = rng.gen_range (0, 3);

    if axis == 0 {
      // x
      (
        [base_pos_a, other1_pos_a, other2_pos_a].into(),
        shape::Cuboid::noisy ([
          base_half_extent_a, other1_half_extent_a, other2_half_extent_a]),
        [base_pos_b, other1_pos_b, other2_pos_b].into(),
        shape::Cuboid::noisy ([
          base_half_extent_b, other1_half_extent_b, other2_half_extent_b])
      )
    } else if axis == 1 {
      // y
      (
        [other1_pos_a, base_pos_a, other2_pos_a].into(),
        shape::Cuboid::noisy ([
          other1_half_extent_a, base_half_extent_a, other2_half_extent_a]),
        [other1_pos_b, base_pos_b, other2_pos_b].into(),
        shape::Cuboid::noisy ([
          other1_half_extent_b, base_half_extent_b, other2_half_extent_b])
      )
    } else {
      // z
      debug_assert_eq!(axis, 2);
      (
        [other1_pos_a, other2_pos_a, base_pos_a].into(),
        shape::Cuboid::noisy ([
          other1_half_extent_a, other2_half_extent_a, base_half_extent_a]),
        [other1_pos_b, other2_pos_b, base_pos_b].into(),
        shape::Cuboid::noisy ([
          other1_half_extent_b, other2_half_extent_b, base_half_extent_b])
      )
    }
  }

} // end tests