//	NEBSpaceGeometry.swift
//
//	© 2026 by Jeff Weeks
//	See TermsOfUse.txt

import Darwin	//	for sin() and sinh()
import simd		//	for simd_double3x3


//	When we access the components of a SIMD3<Double> p,
//	as .x, .y and .z, the Swift compiler compiles our code quickly.
//	However, when working with a point p in the pool table's
//	XZW coordinates (with no Y !), it would be super confusing to let
//
//		p.x be the x coordinate
//		p.y be the z coordinate
//		p.z be the w coordinate
//
//	I tried writing p[0], p[1], p[2] instead, but they slowed
//	Swift's type inference system to a crawl at compile time.
//
//	Writing those as p[Int(0)], p[Int(1)], p[Int(2)] lets the code
//	compile faster, but leaves it ugly and almost unreadable.
//
//	So instead let's define a wrapper for a SIMD3<Double>
//	that lets us access the components as .x, .z and .w (no .y).
struct VectorXZW {

	var v: SIMD3<Double>
	
	init(
		_ x: Double,
		_ z: Double,
		_ w: Double
	) {
		self.v = SIMD3<Double>(x, z, w)
	}

	init(
		_ v: SIMD3<Double>
	) {
		self.v = v
	}
	
	var x: Double {
		get { return v.x }
		set { v.x = newValue }
	}
	var z: Double {
		get { return v.y }	//	v.y stores VectorXZW's z coordinate (it has no y coordinate)
		set { v.y = newValue }
	}
	var w: Double {
		get { return v.z }	//	v.z stores VectorXZW's w coordinate (not its z coordinate)
		set { v.z = newValue }
	}
}
prefix func - (vector: VectorXZW) -> VectorXZW {
	return VectorXZW(-vector.v)
}
func *(lhs: Double, rhs: VectorXZW) -> VectorXZW {
	return VectorXZW(lhs * rhs.v)
}
func -(lhs: VectorXZW, rhs: VectorXZW) -> VectorXZW {
	return VectorXZW(lhs.v - rhs.v)
}
func *(lhs: matrix_double3x3, rhs: VectorXZW) -> VectorXZW {
	return VectorXZW(lhs * rhs.v)
}

extension NEBRenderer {

	class NEBSpaceGeometry: SpaceGeometry {

		//	The spaceRadius and the tableEdgeLengthInMeters
		//	each determines the other, because a right-angled table's
		//	edge lengths in radians are fixed constants.
		//	We treat tableEdgeLengthInMeters as the "independent variable"
		//	and use it to compute the spaceRadius, because
		//	the tableEdgeLengthInMeters is the quantity that
		//	the player perceives directly.
		
		//	On the one hand, the table is ridiculously small
		//	when tableEdgeLengthInMeters = 0.5 or even 0.75,
		//	and in spherical geometry it's hard to slide
		//	the cue stick along itself for a successful shot.
		//	On the other hand, it feels great to play on a 1.5 m table,
		//	but alas because of visionOS's 1.5 m "safe circle",
		//	the player frequently sees the passthrough video,
		//	even with tableEdgeLengthInMeters = 1.25.
		//	So really tableEdgeLengthInMeters = 1.0 is pretty much
		//	the only practical choice.
		//
		//		Note: If Apple ever allows a larger radius
		//		for the safe circle, I can re-install my code
		//		for a user-adjustable tableEdgeLengthInMeters,
		//		which I've archives in the file
		//
		//			scraps/adjustable tableEdgeLengthInMeters.swift
		//
		let tableEdgeLengthInMeters: Double = 1.0
		
		//	The geometryType and the tableEdgeLengthInMeters
		//	fully determine all of the following variables,
		//	which get computed in the init() given below.
		
		let ballCircumference: Double		//	in meters
		let ballShadowRadius: Double		//	in meters
		let ballLift: simd_double4x4		//	dimensionless (lift in radians)
		let pocketPoses: [simd_double4x4]	//	in table coordinates, ignoring the table height; dimensionless
		let pocketCentersXZW: [VectorXZW]	//	in meters in XZW space
		let tableDistanceFunction:			//	a function
			(
				_ pointA_XZW: VectorXZW,	//	in meters
				_ pointB_XZW: VectorXZW		//	in meters
			) -> Double							//	returns the distance between points A and B in meters,
												//	as measured along a geodesic on the billiards table
		let cushionToPointTableDistanceFunction:	//	a function
			(
				_ cushionNormalXZW: VectorXZW,	//	dimensionless, unit length, inward pointing
												//	normal in XZW space to plane through origin
												//	that contains the cushion
				_ pointXZW: VectorXZW			//	in meters
			) -> Double		//	returns the distance from cushion to point in meters,
							//	as measured along a geodesic on the billiards table
		let northPoleXZW: VectorXZW	//	in meters
		
		init(
			geometryType: GeometryType
		) {

			//	How many meters correspond to 1 radian?
			//
			//	To find out, divide the table edge length in meters
			//	by the inherent table edge length in radians.
			//
			let theSpaceRadius = tableEdgeLengthInMeters
							   / geometryType.tableEdgeLengthInRadians

			//	The geometryType and tableEdgeLengthInMeters never change
			//	while the Billiards Immersive Space is running, so we may
			//	assign them to "let" constants and precompute the various
			//	other properties that depend on them.
					
			ballCircumference = NEBSpaceGeometry.computeCircleCircumference(
					radius: ballRadius,
					geometryType: geometryType,
					spaceRadius: theSpaceRadius)
			
			ballShadowRadius = NEBSpaceGeometry.computeBallShadowRadius(
					geometryType: geometryType,
					spaceRadius: theSpaceRadius)
			
			ballLift = NEBSpaceGeometry.computeBallLift(
					geometryType: geometryType,
					spaceRadius: theSpaceRadius)
			
			let thePocketPoses = NEBSpaceGeometry.computePocketPoses(
					geometryType: geometryType,
					spaceRadius: theSpaceRadius)

			pocketPoses = thePocketPoses
			
			pocketCentersXZW = NEBSpaceGeometry.computePocketCenters(
					spaceRadius: theSpaceRadius,
					pocketPoses: thePocketPoses)

			tableDistanceFunction = NEBSpaceGeometry.computeTableDistanceFunction(
					geometryType: geometryType,
					spaceRadius: theSpaceRadius)

			cushionToPointTableDistanceFunction = NEBSpaceGeometry.computeCushionToPointTableDistanceFunction(
					geometryType: geometryType,
					spaceRadius: theSpaceRadius)
			
			northPoleXZW = VectorXZW(0.0, 0.0, theSpaceRadius)
			
			super.init(
				geometryType: geometryType,
				spaceRadius: theSpaceRadius)
		}
		
		static func computeCircleCircumference(
			radius: Double,	//	in meters
			geometryType: GeometryType,
			spaceRadius: Double
		) -> Double {	//	in meters
		
			let theRadiusInSpaceRadians = radius / spaceRadius
			
			let theEffectiveRadiusInSpaceRadians: Double
			switch geometryType {
			
			case .spherical:
				theEffectiveRadiusInSpaceRadians = sin(theRadiusInSpaceRadians)
				
			case .euclidean:
				theEffectiveRadiusInSpaceRadians = theRadiusInSpaceRadians
				
			case .hyperbolic:
				theEffectiveRadiusInSpaceRadians = sinh(theRadiusInSpaceRadians)
			}
			
			let theEffectiveRadiusInMeters = theEffectiveRadiusInSpaceRadians * spaceRadius
			let theCircumferenceInMeters = 2.0 * Double.pi * theEffectiveRadiusInMeters
			
			return theCircumferenceInMeters
		}
		
		private static func computeBallShadowRadius(
			geometryType: GeometryType,
			spaceRadius: Double
		) -> Double {	//	in meters
		
			//	A few definitions:
			//
			//		Let "orthogonal projection" be the mapping
			//		that sends each point of a billiards ball
			//		to the nearest point on the table.
			//		(I'm not sure whether this term is commonly used
			//		in non-euclidean spaces or not, but this is
			//		how we'll use it.)
			//
			//		A ball's image under orthogonal projection
			//		is its "orthogonal shadow". The orthogonal shadow
			//		is, obviously, a 2-dimensional disk in the plane
			//		of the table.
			//
			//		The "shadow radius" is the radius
			//		of a ball's orthogonal shadow.
			//
			//	Note that
			//
			//		On a 3-sphere, a ball's shadow radius is slightly
			//			greater than the ball's true 3-dimensional radius.
			//
			//		In Euclidean space, a ball's shadow radius
			//			equals the ball's true 3-dimensional radius.
			//
			//		In hyperbolic space, a ball's shadow radius is slightly
			//			less than the ball's true 3-dimensional radius.
			//
			//
			//	The question is, how to compute a ball's shadow radius?
			//	Consider two balls touching each other:
			//
			//		    *   *         #   #
			//		 *         *   #         #
			//		*     C₀_   * #   __C₁    #
			//		*     |  `--*T#--´  |     #
			//		*     |     *|#     |     #
			//		 *    |    * | #    |    #
			//		    * P₀*____S____# P₁#
			//
			//	C₀ and C₁ are the ball's centers,
			//	P₀ and P₁ are their orthogonal projections onto the table.
			//	T is the point of tangency where the balls touch, and
			//	S is its orthogonal projection onto the table.
			//
			//	Take a careful look at the quadrilateral C₀P₀ST:
			//
			//		The angle at vertex P₀ must be a right angle
			//		because of the mirror symmetry (reflection
			//		across a vertical line through P₀).
			//
			//		The angle at S must also be a right angle
			//		because of mirror symmetry (reflection
			//		across a vertical line through S).
			//
			//		The angle at T must be a right angle
			//		because the radius of a circle
			//		(in this case the segment C₀T) is always
			//		orthogonal to a line tangent to the circle
			//		(in this case the line containing S and T).
			//
			//		The angle at C₀ is the sensitive one.
			//		In spherical geometry,  it will be greater than 90°.
			//		In Euclidean geometry,  it will exactly equal 90°.
			//		In hyperbolic geometry, it will be less than 90°.
			//
			//	To solve for the dimensions of this quadrilateral,
			//	the crucial observation is that it's symmetrical
			//	about the diagonal C₀S.
			//
			//		    *   *
			//		 *         *
			//		*     C₀_   *
			//		*     | \`--*T
			//		*     |  \_ *|
			//		 *    |    \ |
			//		    * P₀*____S
			//
			//	Proof: The sides C₀P₀ and C₀T have the same length
			//	because both are equal to the ball's radius ρ.
			//	The angles at P₀ and at T are equal because both
			//	are right angles. Therefore the quadrilateral
			//	is symmetrical about the angle bisector at C₀,
			//	which becomes the diagonal C₀S. QED
			//
			//	So really we need consider only the triangle
			//
			//		      C₀
			//		      | \_
			//		  R/ρ |   \_
			//		      |     \
			//		      P₀_____S
			//          90°  r/ρ 45°
			//
			//	where
			//
			//		 ρ  is the space radius in meters (which we know)
			//		 R  is the ball's radius in meters (which we know)
			//		R/ρ is the ball's radius in radians (of the 3D space)
			//		 r  is the ball's shadow radius in meters (which we want to know)
			//		r/ρ is the ball's shadow radius in radians (of the 3D space)
			//
			//	We know one edge length (namely R/ρ) and two angles (90° and 45°)
			//	and we want to solve for one of the other edge lengths (namely r/ρ).
			//
			//	Note:  For a typical space radius on the order of ~1 meter,
			//	the difference between the shadow radius and the ball radius
			//	is only ~10 µm, which is less than the width of a human hair.
			//	The difference would become significant only with a much
			//	smaller table or a much larger ball.
			
			let theBallRadiusInRadians = ballRadius / spaceRadius	//	in radians of the 3D space
			
			let theShadowRadiusInRadians: Double
			switch geometryType {
			
			case .spherical:
				
				//	One of Napier's rules for a right spherical triangle
				//	says that
				//
				//		tan(a) = tan(A) sin(b)
				//
				//	where a, b and c are the lengths
				//	of the spherical triangle's sides in radians,
				//	A, B and C are the opposite angles,
				//	and C is a right angle.
				//	In the present case this becomes
				//
				//		tan(R/ρ) = tan(45°) sin(r/ρ)
				//	⇒
				//		tan(R/ρ) = sin(r/ρ)
				//
				//	Note:
				//		R/ρ ≪ π/4, so tan(R/ρ) ≪ 1 and there's
				//		no risk of a domain error in asin().
				//
				theShadowRadiusInRadians = asin(tan(theBallRadiusInRadians))
				
			case .euclidean:
				
				//	The angle at C₀ is a 45⁰ angle
				//	in the triangle shown above,
				//	and the triangle is isoceles.
				//	Therefore r/ρ = R/ρ.
				//
				theShadowRadiusInRadians = theBallRadiusInRadians
				
			case .hyperbolic:
				
				//	The hyperbolic version of the relevant Napier rule is
				//
				//		tanh(a) = tan(A) sinh(b)
				//
				//	which in the present case becomes
				//
				//		tanh(R/ρ) = tan(45°) sinh(r/ρ)
				//	⇒
				//		tanh(R/ρ) = sinh(r/ρ)
				//
				theShadowRadiusInRadians = asinh(tanh(theBallRadiusInRadians))
			}

			let theShadowRadiusInMeters = spaceRadius * theShadowRadiusInRadians

			return theShadowRadiusInMeters
		}
		
		private static func computeBallLift(
			geometryType: GeometryType,
			spaceRadius: Double
		) -> simd_double4x4 {	//	dimensionless (lift in radians)

			let theBallLift = axisAlignedTranslation(
								axis: 1,
								distanceInRadians: ballRadius/spaceRadius,
								geometryType: geometryType)
			
			return theBallLift
		}
		
		private static func computePocketPoses(
			geometryType: GeometryType,
			spaceRadius: Double
		) -> [simd_double4x4] {	//	in table coordinates, ignoring the table height; dimensionless

			let r = geometryType.tableOutradius	//	in radians
					- pocketInset / spaceRadius

			let theRadialTranslation = axisAlignedTranslation(
										axis: 2,
										distanceInRadians: r,
										geometryType: geometryType)

			let n = geometryType.tableNumSides
			
			let thePocketPoses = (0 ..< n).map() { i in

				let φ = 2.0 * Double.pi * ( -0.5 + Double(i) ) / Double(n)

				let theRotation = axisAlignedRotation(
									from: 2,
									towards: 0,
									angle: φ)
				
				let thePose = theRotation	//	right-to-left matrix action
							* theRadialTranslation
				
				return thePose
			}
			
			return thePocketPoses
		}
		
		private static func computePocketCenters(
			spaceRadius: Double,
			pocketPoses:  [simd_double4x4]	//	in table coordinates, ignoring the table height;
											//		dimensionless
		) -> [VectorXZW] {	//	in XZW space, in meters
		
			let theNorthPole = SIMD4<Double>(0.0, 0.0, 0.0, spaceRadius)

			let thePocketCenters = pocketPoses.map() { pocketPose in
				
				let thePocketCenter4D = pocketPose * theNorthPole
				let thePocketCenter = VectorXZW(
										thePocketCenter4D.x,
										//	omit the y component, which should be zero
										thePocketCenter4D.z,
										thePocketCenter4D.w)
				return thePocketCenter
			}
			
			return thePocketCenters
		}
		
		private static func computeTableDistanceFunction(
			geometryType: GeometryType,
			spaceRadius: Double
		) -> (								//	a function:
				_ pointA_XZW: VectorXZW,	//	in meters
				_ pointB_XZW: VectorXZW		//	in meters
			 ) -> Double					//	returns the distance between points A and B in meters,
											//	as measured along a geodesic on the billiards table
		{
			switch geometryType {
			
			case .spherical:

				return {
					(_ pointA_XZW: VectorXZW, _ pointB_XZW: VectorXZW) -> Double in
					
					//	In the 4-dimensional Euclidean embedding space,
					//	pointA_XZW and pointB_XZW both have length
					//	equal to the spaceRadius. For conciseness, let
					//
					//		ρ denote the spaceRadius
					//		d denote the great-circle distance between the points, in meters.
					//
					//	Taking the dot product of A and B gives
					//
					//		A·B = |A| |B| cos(d/ρ)
					//			= ρ² cos(d/ρ)
					//	so
					//		d = ρ acos(A·B / ρ²)
					//

					let theDotProduct: Double
						= pointA_XZW.x * pointB_XZW.x
						+ pointA_XZW.z * pointB_XZW.z
						+ pointA_XZW.w * pointB_XZW.w
										
					let theNormalizedDotProduct = theDotProduct / (spaceRadius * spaceRadius)
					
					//	Guard against theNormalizedDotProduct = 1.0000001 or similar.
					let theSafeDotProduct = min(max(theNormalizedDotProduct, -1.0), +1.0)
						
					let theDistanceInRadians = acos(theSafeDotProduct)
					
					let theDistanceInMeters = spaceRadius * theDistanceInRadians
					
					return theDistanceInMeters
				}

			case .euclidean:

				return {
					(_ pointA_XZW: VectorXZW, _ pointB_XZW: VectorXZW) -> Double in
					
					let Δp = SIMD2<Double>(
						pointB_XZW.x - pointA_XZW.x,
						pointB_XZW.z - pointA_XZW.z
				// omit pointB_XZW.w - pointA_XZW.w = 0.0
					)
					
					return length(Δp)
				}
				
			case .hyperbolic:

				return {
					(_ pointA_XZW: VectorXZW, _ pointB_XZW: VectorXZW) -> Double in
					
					//	The hyperbolic case is similar to the spherical case,
					//	except that the hyperbolic 3-space is embedded in
					//	a (3+1)-dimensional Minkowski space.
					//	Taking the dot product of A and B gives
					//
					//		A·B = - |A| |B| cosh(d/ρ)
					//			= - ρ² cosh(d/ρ)
					//	so
					//		d = ρ acosh( - A·B / ρ² )
					//

					let theDotProduct: Double
						= pointA_XZW.x * pointB_XZW.x
						+ pointA_XZW.z * pointB_XZW.z
						- pointA_XZW.w * pointB_XZW.w
					
					let theNormalizedDotProduct = -theDotProduct / (spaceRadius * spaceRadius)	//	note the minus sign
					
					//	Guard against theNormalizedDotProduct = 0.9999999 or similar.
					let theSafeDotProduct = max(theNormalizedDotProduct, +1.0)
						
					let theDistanceInRadians = acosh(theSafeDotProduct)
					
					let theDistanceInMeters = spaceRadius * theDistanceInRadians
					
					return theDistanceInMeters
				}
			}
		}
		
		private static func computeCushionToPointTableDistanceFunction(
			geometryType: GeometryType,
			spaceRadius: Double
		) -> (								//	a function:
				_ cushionNormalXZW: VectorXZW,	//	dimensionless, unit length, inward pointing
												//	normal in XZW space to plane through origin
												//	that contains the cushion
				_ pointXZW: VectorXZW			//	in meters
			 ) -> Double	//	distance from cushion to point in meters,
							//	as measured along a geodesic on the billiards table
		{
			switch geometryType {
			
			case .spherical:
			
				return {
					(_ cushionNormalXZW: VectorXZW, _ pointXZW: VectorXZW) -> Double in
					
					let theDotProduct: Double
						= cushionNormalXZW.x * pointXZW.x
						+ cushionNormalXZW.z * pointXZW.z
						+ cushionNormalXZW.w * pointXZW.w

					let theNormalizedDotProduct = theDotProduct / spaceRadius

					//	Guard against theDotProduct = 1.0000001 or similar.
					let theSafeDotProduct = min(max(theNormalizedDotProduct, -1.0), +1.0)
						
					let theDistanceInRadians = asin(theSafeDotProduct)
					
					let theDistanceInMeters = spaceRadius * theDistanceInRadians
					
					return theDistanceInMeters
				}
				
			case .euclidean:
			
				return {
					(_ cushionNormalXZW: VectorXZW, _ pointXZW: VectorXZW) -> Double in
					
					//	The Euclidean case is most easily understood
					//	by looking at it in cross-section, focusing
					//	on the x and w coordinates while ignoring z.
					//	For brevity, let ρ denote the spaceRadius.
					//
					//		  w
					//		  ·   line=      P=
					//		  ·     (ℓ,ρ)   (p_x,ρ)
					//		--·-------•--------•-----
					//		  ·  ℓ   /
					//		  ·     /
					//		ρ ·    /
					//		  ·   /
					//		  ·  /
					//		  · /
					//		  ·/    1
					//		· + · · · · · · · · · · x
					//		  ·\___      |
					//		  ·    \__   | ℓ/ρ
					//		  ·       \__|
					//		  ·          L=(1,-ℓ/ρ)
					//		  ·
					//
					//	In this cross-sectional view, the line itself appears as a dot
					//	at the point (ℓ,ρ) while its normal vector is L = (L_x, L_z).
					//	For simplicity, normalize L to have L_x = 1. We can use
					//	the similar triangles in the above diagram to solve for L_z:
					//
					//		 ℓ    -L_z
					//		--- = -----
					//		 ρ      1
					//	⇒
					//		L_z = -ℓ/ρ
					//
					//	So L = (L_x, L_z) = (1,-ℓ/ρ).
					//
					//	Thus the dot product of L with any point P = (p_x, ρ)
					//	(using the usual Euclidean metric) gives
					//
					//		L·P = (1,-ℓ/ρ)·(p_x, ρ)
					//			= (1)(p_x) + (-ℓ/ρ)(ρ)
					//			= p_x - ℓ
					//
					//	which is exactly the Euclidean distance
					//	from the point P to the line, as desired.
					//

					let theDotProduct: Double
						= cushionNormalXZW.x * pointXZW.x
						+ cushionNormalXZW.z * pointXZW.z
						+ cushionNormalXZW.w * pointXZW.w

					let theDistanceInMeters = theDotProduct
					
					return theDistanceInMeters
				}
				
			case .hyperbolic:
			
				return {
					(_ cushionNormalXZW: VectorXZW, _ pointXZW: VectorXZW) -> Double in
					
					let theDotProduct: Double
						= cushionNormalXZW.x * pointXZW.x
						+ cushionNormalXZW.z * pointXZW.z
						- cushionNormalXZW.w * pointXZW.w

					let theNormalizedDotProduct = theDotProduct / spaceRadius
					
					let theDistanceInRadians = asinh(theNormalizedDotProduct)
					
					let theDistanceInMeters = spaceRadius * theDistanceInRadians
					
					return theDistanceInMeters
				}
				
			}
		}
	}
}

extension GeometryType {
	
	var playButtonTitle: String {
	
		switch self {
		case .spherical:	return "Play spherical pool"
		case .euclidean:	return "Play Euclidean pool"
		case .hyperbolic:	return "Play hyperbolic pool"
		}
	}

	var tableNumSides: Int {
	
		switch self {
		
		case .spherical:
			return 3	//	right equilateral triangle
			
		case .euclidean:
			return 4	//	square
			
		case .hyperbolic:
			return 5	//	right-angled regular pentagon
		}
	}

	//	Billiards table dimensions in radians
	//
	//	The billiards table bed will be a regular right-angled n-gon,
	//	for n = 3, 4 or 5. That is, it will be
	//
	//		- an equilateral right triangle in spherical geometry
	//		- a square in Euclidean geometry
	//		- a regular right-angled pentagon in hyperbolic geometry
	//
	//	In all cases, the regular right-angled n-gon may be
	//	tiled by six smaller triangles in the natural way,
	//	so each smaller triangle has angles π/2, π/4 and π/n.

	var tableInradius: Double {	//	in radians
	
		switch self {
		
		case .spherical:
			//	To compute the big triangle's inradius r,
			//	apply the "spherical law of cosines for angles"
			//	to one of those six smaller triangles:
			//
			//		cos(γ) = -cos(α) cos(β) + sin(α) sin(β) cos(c)
			//
			//		cos(π/4) = -cos(π/2) cos(π/3) + sin(π/2) sin(π/3) cos(r)
			//
			//		√½ = -(0)(½) + (1)(√¾) cos(r)
			//
			//		cos(r) = √½ / √¾
			//
			//		r = acos(√⅔)
			//
			return 0.61547970867038734107	//	acos(√⅔)

		case .euclidean:
			return 0.5 * tableEdgeLengthInRadians
			
		case .hyperbolic:
			//	To compute the pentagon's inradius r,
			//	apply the "hyperbolic law of cosines for angles"
			//	to one of those six smaller triangles:
			//
			//		cos(γ) = -cos(α) cos(β) + sin(α) sin(β) cosh(c)
			//
			//		cos(π/4) = -cos(π/2) cos(π/5) + sin(π/2) sin(π/5) cosh(r)
			//
			//		√½ = -(0)(cos(π/5)) + (1)(sin(π/5)) cosh(r)
			//
			//		cosh(r) = √½ / sin(π/5)
			//
			//		r = acosh( √½ / sin(π/5) )
			//
			return 0.62686966290617781414	//	acosh( √½ / sin(π/5) )
		}
	}
	
	var tableOutradius: Double {	//	in radians
	
		switch self {
		
		case .spherical:
			//	To compute the big triangle's outradius R, apply
			//	the "spherical law of cosines for angles" as for the tableInradius,
			//	but with a different assignment of the angles α, β and γ:
			//
			//		cos(γ) = -cos(α) cos(β) + sin(α) sin(β) cos(c)
			//
			//		cos(π/2) = -cos(π/3) cos(π/4) + sin(π/3) sin(π/4) cos(R)
			//
			//		0 = -(½)(√½) + (√¾)(√½) cos(R)
			//
			//		cos(R) = ( (½)(√½) ) / ( (√¾)(√½) )
			//
			//		cos(R) = 1/√3
			//
			//		R = acos(√⅓)
			//
			return 0.95531661812450927816	//	acos(√⅓)
			
		case .euclidean:
			return (0.5).squareRoot() * tableEdgeLengthInRadians
			
		case .hyperbolic:
			//	To compute the pentagon's outradius R, apply
			//	the "hyperbolic law of cosines for angles" as for the tableInradius,
			//	but with a different assignment of the angles α, β and γ:
			//
			//		cos(γ) = -cos(α) cos(β) + sin(α) sin(β) cosh(c)
			//
			//		cos(π/2) = -cos(π/5) cos(π/4) + sin(π/5) sin(π/4) cosh(R)
			//
			//		0 = -(cos(π/5))(√½) + (sin(π/5))(√½) cosh(R)
			//
			//		cosh(R) = ( (cos(π/5))(√½) ) / ( (sin(π/5))(√½) )
			//
			//		cosh(R) = cot(π/5)
			//
			//		R = acosh(cot(π/5))
			//
			return 0.84248208146200745911	//	acosh(cot(π/5))
		}
	}
	
	var tableEdgeLengthInRadians: Double {
	
		switch self {
		
		case .spherical:
			//	To compute half the big triangle's edge length s,
			//	apply the "spherical law of cosines for angles"
			//	as for the tableInradius and tableOutradius,
			//	with yet another assignment of the angles α, β and γ:
			//
			//		cos(γ) = -cos(α) cos(β) + sin(α) sin(β) cos(c)
			//
			//		cos(π/3) = -cos(π/2) cos(π/4) + sin(π/2) sin(π/4) cos(s/2)
			//
			//		½ = (0)(√½) + (1)(√½) cos(s/2)
			//
			//		cos(s/2) = ½ / √½
			//
			//		cos(s/2) = √½
			//
			//		s/2 = π/4
			//
			//		s = π/2
			//
			//	Note:  The answer s = π/2 agrees with the edge length
			//	one gets by examining an octahedral triangulation
			//	of a 2-sphere (and skipping the above computation!).
			//
			return 1.57079632679489661923	//	π/2
			
		case .euclidean:
			//	We may set the table's edge length in radians to any value
			//	we want. Our init() uses it to convert the player's
			//	requested tableEdgeLengthInMeters to the corresponding
			//	space radius
			//
			//		                edge length in meters
			//		space radius = ------------------------
			//		                edge length in radians
			//
			//	So the table and the lattice end up looking the same
			//	no matter what value we choose for the tableEdgeLengthInRadians.
			//
			//	(Note: This works in Euclidean geometry because Euclidean
			//	geometry has similarities. It would definitely not work
			//	in spherical or hyperbolic geometry.)
			//
			//	The simplest choice would be to set the edge length
			//	in radians to 1. The only risk there is that it tempts us
			//	(well, it tempted me!) into hard-coding that assumption
			//	into our constructions of the checkerboard poses' translations
			//	and the checkerboard cell's walls. To avoid that temptation,
			//	let's set the edge length in radians to some value
			//	other than 1.
			//
			return 0.75	//	radians
			
		case .hyperbolic:
			//	To compute half the pentagon's edge length s,
			//	apply the "hyperbolic law of cosines for angles"
			//	as for the tableInradius and tableOutradius,
			//	with yet another assignment of the angles α, β and γ:
			//
			//		cos(γ) = -cos(α) cos(β) + sin(α) sin(β) cosh(c)
			//
			//		cos(π/5) = -cos(π/2) cos(π/4) + sin(π/2) sin(π/4) cosh(s/2)
			//
			//		cos(π/5) = (0)(√½) + (1)(√½) cosh(s/2)
			//
			//		cosh(s/2) = cos(π/5) / √½
			//
			//		s/2 = acosh( cos(π/5) / √½ )
			//
			//		s = 2 acosh( cos(π/5) / √½ )
			//
			return 1.06127506190503565203	//	2 acosh( cos(π/5) / √½ )
		}
	}
	
	var cushionNormalVectorsXZW: [VectorXZW] {
		//	Returns an Array of dimensionless, unit length, inward pointing normals
		//	to the planes through origin in XZW space that contain the cushions.

		//	Note:  The following code and the code in tableMesh()
		//	share the assumptions that edge 0 sits with its midpoint
		//	on the positive z axis, and that the remaining edges
		//	progress around in the direction of the positive x axis.

		let θ = tableInradius	//	in radians
		let n = tableNumSides

		let theUnrotatedNormalXZW: VectorXZW
		switch self {
		
		case .spherical:
			theUnrotatedNormalXZW = VectorXZW(0.0, -cos(θ), sin(θ))
			
		case .euclidean:
			//	See comment in cushionToPointTableDistanceFunction()
			//	for an illustrated explanation of the normal vector
			//	that we're using in the Euclidean case.
			theUnrotatedNormalXZW = VectorXZW(0.0, -1.0, θ)
			
		case .hyperbolic:
			theUnrotatedNormalXZW = VectorXZW(0.0, -cosh(θ), -sinh(θ))
		}
		
		let Δφ = 2.0 * Double.pi / Double(n)
		let theRotatedNormalsXZW = (0 ..< n).map { i in
		
			let φ = Double(i) * Δφ
		
			let theRotation = simd_double3x3(
				SIMD3<Double>( cos(φ), -sin(φ), 0.0),
				SIMD3<Double>(+sin(φ),  cos(φ), 0.0),
				SIMD3<Double>(  0.0,     0.0,   1.0)
			)
			
			let theRotatedNormalXZW = theRotation	//	right-to-left matrix action
									* theUnrotatedNormalXZW
			
			return theRotatedNormalXZW
		}
		
		return theRotatedNormalsXZW
	}


// MARK: -
// MARK: 3D metric (XZW space)

	
	var xzwMetricForVertical: XZWMetric {	//	Euclidean needs  "vertical"  00+ metric
		switch self {
		case .spherical:	return .Euclidean	//	signature +++
		case .euclidean:	return .vertical	//	signature 00+
		case .hyperbolic:	return .MinkowskiT	//	signature --+
		}
	}
	
	var xzwMetricForHorizontal: XZWMetric {	//	Euclidean needs "horizontal" ++0 metric
		switch self {
		case .spherical:	return .Euclidean	//	signature +++
		case .euclidean:	return .horizontal	//	signature ++0
		case .hyperbolic:	return .MinkowskiS	//	signature ++-
		}
	}
}

enum XZWMetric {

	case Euclidean	//	signature +++
	case MinkowskiS	//	signature ++-  (spacelike vectors have positive squared length)
	case MinkowskiT	//	signature --+  ( timelike vectors have positive squared length)
	case horizontal	//	signature ++0
	case vertical	//	signature 00+
	
	func innerProduct(
		_ u: VectorXZW,
		_ v: VectorXZW
	) -> Double {

		switch self {
		case .Euclidean:	return   u.x * v.x  +  u.z * v.z  +  u.w * v.w
		case .MinkowskiS:	return   u.x * v.x  +  u.z * v.z  -  u.w * v.w
		case .MinkowskiT:	return -(u.x * v.x  +  u.z * v.z) +  u.w * v.w
		case .horizontal:	return   u.x * v.x  +  u.z * v.z
		case .vertical:		return                               u.w * v.w
		}
	}
	
	func length(
		_ u: VectorXZW
	) -> Double {
	
		return sqrt(abs(innerProduct(u, u)))
	}
	
	func normalize(
		_ u: VectorXZW
	) -> VectorXZW {
	
		let theLength = length(u)

		//	Don't try to normalize the unnormalizable.
		let ε = 1e-4
		if theLength < ε {
			assertionFailure("normalize() received a vector of length less than ε")
			return u
		}
		
		let theNormalizedVector = (1.0/theLength) * u
		
		return theNormalizedVector
	}
}
