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

import Metal


extension NEBRenderer {

	struct NEBVertex {
		let pos: SIMD4<Float32>	//	position (x,y,z,w)
		let tex: SIMD3<Float32>	//	texture coordinates (u,v,w) or (u,v,0)
	}

	class NEBMesh {
		
		//	To index the mesh vertices, use UInt16, which makes
		//	good use of the GPU's power-efficient 16-bit arithmetic
		//	and allows us to create meshes with up to 2¹⁶ = 65536 vertices,
		//	which is more than enough for our needs.
											
		let vertexCount: Int
		let indexCount: Int

		let vertexBuffer: MTLBuffer	//	array of NEBVertex
		let indexBuffer: MTLBuffer	//	array of UInt32
		
		init(
			device: MTLDevice,
			vertexCount: Int,
			indexCount: Int)
		{
			self.vertexCount = vertexCount
			self.indexCount = indexCount
			
			guard let theVertexBuffer = device.makeBuffer(
				length: vertexCount * MemoryLayout<NEBVertex>.stride,
				options: [MTLResourceOptions.storageModeShared]
			) else {
				fatalError("Couldn't make the vertexBuffer in NEBMesh()")
			}
			self.vertexBuffer = theVertexBuffer
			
			guard let theIndexBuffer = device.makeBuffer(
				length: indexCount * MemoryLayout<UInt16>.stride,
				options: [MTLResourceOptions.storageModeShared]
			) else {
				fatalError("Couldn't make the indexBuffer in NEBMesh()")
			}
			self.indexBuffer = theIndexBuffer
		}
	}


// MARK: -
// MARK: Table

	func makeTableMesh(
		device: MTLDevice
	) -> NEBMesh {

		let theGeometryType = spaceGeometry.geometryType
		let theSpaceRadius = spaceGeometry.spaceRadius

		//	The billiards table has 3, 4 or 5 corners in the
		//	spherical, Euclidean and hyperbolic cases respectively.
		//
		let n = theGeometryType.tableNumSides

		//	We triangulate it as 3, 4 or 5 large triangles,
		//	each of which is subdivided into smaller triangles like this:
		//
		//		        *						|
		//		       / \						|
		//		      *———* 					|
		//		     / \ / \					|
		//		    *———*———*					|
		//		   / \ / \ / \					|
		//		  *———*———*———*					|
		//		 / \ / \ / \ / \				|
		//		*———*———*———*———*				|
		//
		//	Let the number of mesh divisions per side be a power of two,
		//	so that we can use a recursive algorithm to position them
		//	more evenly on a curved surface.
		//
		let theNumRecursiveSubdivisions = 5
		let m = 1 << theNumRecursiveSubdivisions	// number of mesh edges per side

		//	All distances will be in radians until the very end.
		
		//	Determine the table's corners in the curved space itself.
		//
		let theNorthPoleInRadians = SIMD4<Double>(0.0, 0.0, 0.0, 1.0)
		let theOutradius = theGeometryType.tableOutradius
		let theRadialTranslation = axisAlignedTranslation(
									axis: 2,	// translate w axis towards z axis
									distanceInRadians: theOutradius,
									geometryType: theGeometryType)
		let theProtoCorner = theRadialTranslation * theNorthPoleInRadians
		let theTrailingCorners = Array(0 ..< n).map() { i in

			//	The following two lines and the code in cushionNormalVectorsXZW()
			//	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 φ = (Double(i) - 0.5) * (2.0 * Double.pi / Double(n))
			let theTrailingRotation = axisAlignedRotation(from: 2, towards: 0, angle: φ)
			let theTrailingCorner = theTrailingRotation * theProtoCorner

			return theTrailingCorner
		}

		//	Barycentric mesh indices
		//
		//	Each of the n large triangles keeps a map telling
		//	where to find the mesh vertices that comprise it.
		//	The map is organized in barycentric coordinates,
		//	for example with 4 mesh-divisions-per-side it looks like
		//
		//					   pivot vertex
		//					   (north pole)
		//		                 (4,0,0)
		//		             (3,1,0) (3,0,1)
		//		         (2,2,0) (2,1,1) (2,0,2)
		//		     (1,3,0) (1,2,1) (1,1,2) (1,0,3)
		//		 (0,4,0) (0,3,1) (0,2,2) (0,1,3) (0,0,4)
		//	"trailing"							   "leading"
		//	  vertex								 vertex
		//

		let theVertexCount =
			1	//	the pivot vertex where all 3, 4 or 5 large triangles meet
		+
			n * m	//	the vertices along the edges where adjacent large triangles meet
					//		(excluding the pivot vertex, which we've already counted)
		+
			n * ( (m - 1) * m / 2 )	//	all remaining vertices
		
		let theIndexCount = 3 * n * m * m	//	n * m * m facets, each with three vertices
		
		//	Allocate an Array to hold the mesh vertices' positions
		//	in the given space (a 3-sphere, a Euclidean 3-space,
		//	or a hyperbolic 3-space), measured in radians.
		//
		var theVertexPositionsInRadians = Array(
											repeating: SIMD4<Double>.zero,
											count: theVertexCount)

		//	Convert a vertex's barycentric coordinates
		//	to its index in the linear array theVertexPositionsInRadians.

		var vertexIndex: [ [ [ Int ] ] ]	//	index in linear array theVertexPositionsInRadians
			= Array(
				repeating:
					Array(
						repeating:
							Array(
								repeating:
									//	The three barycentric coordinates always sum
									//	to exactly m, so the last one may be left implicit,
									//	to reduce this array's memory requirements
									//	by a factor of (m + 1).
									Int.max,	//	dummy value
								count: m + 1),	//	second barycentric coordinate
						count: m + 1),	//	first barycentric coordinate
				count: n)	//	which of the n large triangles


			//	Keep a running count of how many indices we've assigned.
		var theRunningVertexCount = 0

			//	Assign an index to the pivot vertex at the north pole.
		for i in 0 ..< n {	//	which large triangle
			vertexIndex[i][m][0]/*[0]*/ = theRunningVertexCount
		}
		theRunningVertexCount += 1

			//	Assign indices to the vertices where exactly two large triangles meet.
		for i in 0 ..< n {	//	which large triangle
		
			let ii = (i + 1) % n
			
			for j in 0 ..< m {	//	weight relative to pivot vertex
								//	(the first coordinate in the diagram above)
								//	Note:  omits j = m
			
				//	Assign an index to the vertex on the
				//	leading edge of this large triangle.
				vertexIndex[i][j][0]/*[m - j]*/ = theRunningVertexCount

				//	Assign the same index to the corresponding vertex
				//	on the trailing edge of the next large triangle.
				vertexIndex[ii][j][m - j]/*[0]*/ = theRunningVertexCount

				theRunningVertexCount += 1
			}
		}

			//	Assign all remaining mesh vertices.
		for i in 0 ..< n {	//	which large triangle

			for j in 0 ..< m - 1 {	//	weight relative to pivot vertex
									//	(the first coordinate in the diagram above)
			
				for k in 1 ..< m - j {	//	weight relative to trailing vertex
										//	(the second coordinate in the diagram above)
										//	k > 0     excludes leading  edge at k = 0
										//	k < m - j excludes trailing edge at ℓ = 0
				
					//	ℓ = m - (j + k) is the weight relative to leading vertex
					//						(the third coordinate in the diagram above)

					//	Assign an index to this vertex.
					vertexIndex[i][j][k]/*[ℓ]*/ = theRunningVertexCount
					theRunningVertexCount += 1
				}
			}
		}
		
		precondition(
			theRunningVertexCount == theVertexCount,
			"Table mesh vertex indexing code found wrong number of vertices")


		//	Define the vertex positions.
		//
		//	First assign each vertex a position in the curved space
		//	measured in radians. In principle we could project
		//	the Euclidean diagram of each big triangle (shown above)
		//	onto the curved space via central projection,
		//	but that would give needlessly nonuniformly-spaced vertices
		//	in the curved space. So instead we'll
		//
		//	  - manually place the "primary" vetrices of each big triangle,
		//			whose indices are all multiples on m, then
		//
		//	  - bisect the "primary" edges to place the "secondary" vertices,
		//			shown in [ ] in the diagram below,
		//			whose indices are all multiples of m/2,
		//
		//		                 (4,0,0)
		//		                -       -
		//		         [2,2,0]         [2,0,2]
		//		        -                       -
		//		 (0,4,0) --------[0,2,2]-------- (0,0,4)
		//
		//	  - bisect the "secondary" edges to place the "tertiary" vertices,
		//			shown in { } in the diagram below,
		//			whose indices are all multiples of m/4,
		//
		//		                 (4,0,0)
		//		             {3,1,0} {3,0,1}
		//		         [2,2,0] {2,1,1} [2,0,2]
		//		     {1,3,0} {1,2,1} {1,1,2} {1,0,3}
		//		 (0,4,0) {0,3,1} [0,2,2] {0,1,3} (0,0,4)
		//
		//	  - and continue bisecting until we've placed even the vertices
		//			whose indices are multiples of 1
		//			(but not of any larger power of 2).
		//
		for i in 0 ..< n {	//	which large triangle
		
			//	Redundancy note:
			//	The pivot vertex will get set n times and
			//	the other vertices along the leading/trailing edges
			//	will each get set twice, but that's OK.

			//	Set the pivot vertex.
			theVertexPositionsInRadians[Int(vertexIndex[i][m][0]/*[0]*/)]
				= theNorthPoleInRadians
			
			//	Set the trailing vertex.
			theVertexPositionsInRadians[Int(vertexIndex[i][0][m]/*[0]*/)]
				= theTrailingCorners[i]
			
			//	Set the leading vertex.
			theVertexPositionsInRadians[Int(vertexIndex[i][0][0]/*[m]*/)]
				= theTrailingCorners[(i+1)%n]
			
			//	Bisect the "primary" edges to place the "secondary" vertices,
			//	and so on, as explained above.
			var theStride = m
			while theStride > 1 {
			
				let theHalfStride = theStride / 2

				for j in stride(from: 0, to: m, by: theStride) {			//	which row
					for k in stride(from: 0, to: m - j, by: theStride) {	//	how far along the row
					
						let ℓ = m - (j + k)
					
						theVertexPositionsInRadians[
								vertexIndex[i][j][k + theHalfStride]/*[ℓ - theHalfStride]*/ ]
							= theGeometryType.normalize(
								theVertexPositionsInRadians[Int(
									vertexIndex[i][j][      k      ]/*[      ℓ      ]*/)]
							  + theVertexPositionsInRadians[Int(
									vertexIndex[i][j][k + theStride]/*[ℓ - theStride]*/)] )
					
						theVertexPositionsInRadians[
								vertexIndex[i][ℓ - theHalfStride][j]/*[k + theHalfStride]*/ ]
							= theGeometryType.normalize(
								theVertexPositionsInRadians[
									vertexIndex[i][      ℓ      ][j]/*[      k      ]*/ ]
							  + theVertexPositionsInRadians[
									vertexIndex[i][ℓ - theStride][j]/*[k + theStride]*/ ] )
					
						theVertexPositionsInRadians[
								vertexIndex[i][k + theHalfStride][ℓ - theHalfStride]/*[j]*/ ]
							= theGeometryType.normalize(
								theVertexPositionsInRadians[
									vertexIndex[i][      k      ][      ℓ      ]/*[j]*/ ]
							  + theVertexPositionsInRadians[
									vertexIndex[i][k + theStride][ℓ - theStride]/*[j]*/ ] )
					}
				}

				//	Halve theStride and take another pass through the while-loop.
				theStride /= 2
			}
		}

		let theVertexPositionsInMeters = theVertexPositionsInRadians.map() { p in
			SIMD4<Float32>(theSpaceRadius * p)
		}


		//	Define the texture coordinates.
		
		//	"Unwrap" the mesh out onto the tangent space
		//	at the north pole (0,0,0,1), using the inverse of the
		//	azimuthal equidistant map, to get texture coordinates
		//	on a scale of radians. We'll factor in the space radius
		//	momentarily.
		let theTangentSpacePositionsInRadians = theVertexPositionsInRadians.map { p in

			let theXYZPart = SIMD3<Double>(p.x, p.y, p.z)
			let theXYZLength = length(theXYZPart)
			if theXYZLength == 0.0 {
				//	May occur only for p.w == +1.
				//	(Our spherical table never reaches p.w == -1.)
				return SIMD3<Double>.zero
			}
			let theDirection = theXYZPart / theXYZLength

			let theDistance: Double	//	in radians
			switch theGeometryType {
			
			case .spherical:
				//	asin(theXYZLength) may be more numerically accurate
				//	than acos(p.w) for small angles.
				theDistance = asin(theXYZLength)	//	= acos(p.w)
				
			case .euclidean:
				theDistance = theXYZLength
				
			case .hyperbolic:
				//	asinh(theXYZLength) may be more numerically accurate
				//	than acosh(p.w) for small angles.
				theDistance = asinh(theXYZLength)	//	= acosh(p.w])
			}
			
			return theDistance * theDirection
		}

		//	Let's let the green felt have random "features" on the scale of 1 mm,
		//	even though such a fine feature scale would, if unmodified, cause
		//	the table texture to flicker in the more distant parts of the table.
		//	To avoid the flicker, the GPU vertex function calls textureContrast()
		//	and the GPU fragment function applies the result to damp the contrast
		//	in the affected regions.
		//
		//		Note: If we ever do any psychology experiment to test
		//		how people perceive the table, we can increase the feature scale
		//		to 2 mm here while increasing rThreshold and rSaturation
		//		in the textureContrast() function to reduce the damping.
		//		But for general use, 2 mm features look unnaturally coarse.
		//		1 mm features look more like a real fabric.
		//
		let theFeatureSizeInMeters = 0.001
		let theFeaturesPerRadian = theSpaceRadius / theFeatureSizeInMeters
		let theTextureCoordinates = theTangentSpacePositionsInRadians.map { p in
		
			//	The point p has y coordinate 0 because it came from a point on the table,
			//	and it also has w coordinate 0 because we just unwrapped it
			//	onto the tangent space at the north pole.
			//	So we may now use (x,z) as texture coordinates.
			//	The radial scale will be constant, but of course there'll
			//	be some variation in the azimuthal scale to accommodate
			//	the space's curvature.
			
			SIMD2<Float32>(theFeaturesPerRadian * SIMD2<Double>(p.x, p.z))
		}

		
		//	Define the mesh facets.
		
		var theFacets = Array(repeating: 0, count: theIndexCount)

		var c = 0	//	the running facet index count

		for i in 0 ..< n {	//	which large triangle

			//	We may create each rightside-up (Δ) facet
			//	by starting at its vertex closest to the pivot vertex
			//	and tracing clockwise around the facet.
			for j in 1 ... m {	//	which row
				for k in 0 ... m - j {		//	how far along the row

					//	ℓ = m - (j + k)
					
					theFacets[c+0] = vertexIndex[i][  j  ][  k  ]/*[  ℓ  ]*/
					theFacets[c+1] = vertexIndex[i][j - 1][  k  ]/*[ℓ + 1]*/
					theFacets[c+2] = vertexIndex[i][j - 1][k + 1]/*[  ℓ  ]*/

					c += 3
				}
			}

			//	We may create each upside-down (∇) facet
			//	by starting at its vertex furthest to the pivot vertex
			//	and tracing clockwise around the facet.
			for j in 0 ..< m - 1 {			//	which row
				for k in 1 ... m - j - 1 {	//	how far along the row

					//	ℓ = m - (j + k)
					
					theFacets[c+0] = vertexIndex[i][  j  ][  k  ]/*[  ℓ  ]*/
					theFacets[c+1] = vertexIndex[i][j + 1][  k  ]/*[ℓ - 1]*/
					theFacets[c+2] = vertexIndex[i][j + 1][k - 1]/*[  ℓ  ]*/

					c += 3
				}
			}
		}
		
		precondition(
			c == theIndexCount,
			"Table mesh facet code found wrong number of facets")


		//	Package up the mesh.
		
		let theMesh = NEBMesh(
						device: device,
						vertexCount: theVertexCount,
						indexCount: theIndexCount)
		
		precondition(theMesh.vertexCount == theVertexCount, "Wrong vertex count in tableMesh()")
		let theVertices = theMesh.vertexBuffer.contents()
			.bindMemory(to: NEBVertex.self, capacity: theVertexCount)
		for i in 0 ..< theVertexCount {
		
			theVertices[i] = NEBVertex(
				pos: theVertexPositionsInMeters[i],
				tex: SIMD3<Float32>(theTextureCoordinates[i], 0.0) )
		}

		precondition(theMesh.indexCount == theIndexCount, "Wrong index count in tableMesh()")
		let theIndices = theMesh.indexBuffer.contents()
			.bindMemory(to: UInt16.self, capacity: theIndexCount)
		for i in 0 ..< theIndexCount {
			theIndices[i] = UInt16(theFacets[i])
		}

		return theMesh
	}


// MARK: -
// MARK: Disk

	func makeDiskMesh(
		device: MTLDevice,
		r: Double,	//	disk radius in meters
		n: Int		//	number of triangles
	) -> NEBMesh {

		let ρ = spaceGeometry.spaceRadius	//	in meters
		let φ = r / ρ	//	disk radius in spatial radians
		
		let c: Double
		let s: Double
		switch spaceGeometry.geometryType {

		case .spherical:
			c = cos(φ)
			s = sin(φ)
			
		case .euclidean:
			c = 1.0
			s = φ

		case .hyperbolic:
			c = cosh(φ)
			s = sinh(φ)
		}

		//	Triangulate the disk using n triangles that all meet
		//	at the disk's center. Even though it would be easy
		//	to triangulate a disk using only vertices on its circumference,
		//	we intentionally include the vertex at the disk's center,
		//	so that the grounding shadow mesh may stay closer to the table's
		//	curved representation when adjustPerceivedDistances is enabled
		//	and we're drawing the disk in the azimuthal equidistant map.
			
		//	Define the vertex positions.
		var thePositions = Array(repeating: SIMD4<Float32>.zero, count: n + 1)
		thePositions[0] = SIMD4<Float32>(0.0, 0.0, 0.0, Float32(ρ))
		for i in 1 ... n {
			let θ = 2.0 * Double.pi * Double(i - 1)/Double(n)	//	2π (i - 1)/n
			thePositions[i] = SIMD4<Float32>( ρ * SIMD4<Double>(s * cos(θ), 0.0, s * sin(θ), c) )
		}
		
		//	Define the texture coordinates.
		//
		//	The first coordinate tells the relative radial distance:
		//		0.0 for the center vertex and
		//		1.0 for each circumferential vertex.
		//
		//	The second coordinate is unused.
		//
		var theTextureCoordinates = Array(repeating: SIMD2<Float32>(1.0, 0.0), count: n + 1)
		theTextureCoordinates[0] = SIMD2<Float32>(0.0, 0.0)	//	center vertex
		
		//	Define the mesh facets.
		var theFacetIndices = Array(repeating: 0, count: 3*n)
		for i in 0 ... n - 1 {
			theFacetIndices[3*i + 0] = i + 1
			theFacetIndices[3*i + 1] = ( i == n - 1 ? 1 : i + 2 )
			theFacetIndices[3*i + 2] = 0
		}

		//	Package up the mesh.
		
		let theMesh = NEBMesh(
						device: device,
						vertexCount: n + 1,
						indexCount: 3*n)
		
		precondition(theMesh.vertexCount == n + 1, "Wrong vertex count in makeDiskMesh()")
		let theVertices = theMesh.vertexBuffer.contents()
			.bindMemory(to: NEBVertex.self, capacity: n + 1)
		for i in 0 ..< n + 1 {
		
			theVertices[i] = NEBVertex(
				pos: thePositions[i],
				tex: SIMD3<Float32>(theTextureCoordinates[i], 0.0) )
		}

		precondition(theMesh.indexCount == 3*n, "Wrong index count in makeDiskMesh()")
		let theIndices = theMesh.indexBuffer.contents()
			.bindMemory(to: UInt16.self, capacity: 3*n)
		for i in 0 ..< 3*n {
			theIndices[i] = UInt16(theFacetIndices[i])
		}

		return theMesh
	}


// MARK: -
// MARK: Ball

	//	What polyhedron should makeBallMesh() start with?
	enum SphereBase {

		//	An icosahedron works best for most purposes
		//	because it distributes the curvature more evenly.
		case icosahedron
		
		//	For the iOS/macOS versions of 4D Maze and the Torus Games'
		//	3D mazes we started with an octahedron instead,
		//	because its intersections with the coordinate planes
		//	are regular 2ⁿ-gons, which mate perfectly with cylindrical tubes
		//	running parallel to the x, y and z axes.
		//
		//		Note:  The octahedron-based mesh also
		//		offers the possibility of drawing only
		//		a half, a quarter or an eighth of the sphere.
		//
		//	I'm unlikely to ever use the octahedron-based mesh defined here,
		//	because I'm doing my Euclidean visionOS games with RealityKit,
		//	and these Compositor Services-based non-Euclidean games
		//	are unlikely to have tubes meeting at right angles to
		//	an octahedron-based sphere.
		case octahedron
	}

	//	To facilitate creating a sphere in a non-Euclidean space,
	//	let's first create a reference unit sphere in Euclidean 3-space.
	private struct Mesh3D {
		let vertexBuffer: [SIMD3<Double>]
		let indexBuffer: [UInt16]
	}

	func makeBallMesh(
		device: MTLDevice,
		r: Double = ballRadius,	//	in meters
		baseShape: SphereBase = .icosahedron,
		refinementLevel: Int = 4	//	refinementLevel = 3 is also fine
	) -> NEBMesh {

		//	Starting with an icosahedron or octahedron, we call
		//	subdivideUnitRadiusPolyhedronMesh() several times
		//	to subdivide the mesh, each time replacing each old facet
		//	with four new ones:
		//
		//		  /\		//
		//		 /__\		//
		//		/_\/_\		//
		//
		//	So if we start with an icosahedron, the refinementLevel
		//	determines the number of faces as follows:
		//
		//		0	  20 faces
		//		1	  80 faces
		//		2	 320 faces
		//		3	1280 faces
		//		etc.
		//
		//	If we instead start with an octahedron,
		//	the number of faces will be
		//
		//		0	   8 faces
		//		1	  32 faces
		//		2	 128 faces
		//		3	 512 faces
		//		etc.
		//
		var theUnitSphere: Mesh3D
		switch baseShape {
		
		case .icosahedron:
			theUnitSphere = unitIcosahedronMesh()
			
		case .octahedron:
			theUnitSphere = unitOctahedronMesh()
		}
		for _ in 0 ..< refinementLevel {
			theUnitSphere = subdivideUnitRadiusPolyhedronMesh(mesh: theUnitSphere)
		}
		
		//	Use theUnitSphere to create a sphere
		//	of the desired radius and spaceGeometry.

		let ρ = spaceGeometry.spaceRadius	//	in meters
		let φ = r / ρ	//	sphere radius in spatial radians

		let c: Double
		let s: Double
		switch spaceGeometry.geometryType {

		case .spherical:
			c = cos(φ)
			s = sin(φ)
			
		case .euclidean:
			c = 1.0
			s = φ

		case .hyperbolic:
			c = cosh(φ)
			s = sinh(φ)
		}

		let theMesh = NEBMesh(
						device: device,
						vertexCount: theUnitSphere.vertexBuffer.count,
						indexCount: theUnitSphere.indexBuffer.count)

		precondition(
			theMesh.vertexCount == theUnitSphere.vertexBuffer.count,
			"Wrong vertex count in makeBallMesh()")
		let theVertices = theMesh.vertexBuffer.contents()
			.bindMemory(to: NEBVertex.self, capacity: theMesh.vertexCount)
		for i in 0 ..< theMesh.vertexCount {
		
			let v = theUnitSphere.vertexBuffer[i]
		
			theVertices[i] = NEBVertex(
				pos: SIMD4<Float32>( ρ * SIMD4<Double>(s * v, c) ),
				tex: SIMD3<Float32>(v) )
		}

		precondition(
			theMesh.indexCount == theUnitSphere.indexBuffer.count,
			"Wrong index count in makeBallMesh()")
		let theIndices = theMesh.indexBuffer.contents()
			.bindMemory(to: UInt16.self, capacity: theMesh.indexCount)
		for i in 0 ..< theMesh.indexCount {
			theIndices[i] = theUnitSphere.indexBuffer[i]
		}

		return theMesh
	}

	private func unitIcosahedronMesh() -> Mesh3D {

		//	The icosahedron's (unnormalized) vertices sit at
		//
		//			( 0, ±1, ±φ)
		//			(±1, ±φ,  0)
		//			(±φ,  0, ±1)
		//
		//	where φ is the golden ratio. The golden ratio is a root
		//	of the irreducible polynomial φ² - φ - 1,
		//	with numerical value φ = (1 + √5)/2 ≈ 1.6180339887…

		let φ = 0.5 * (1.0 + sqrt(5.0))
			//	golden ratio
			//	≈ 1.61803398874989484820
		
		let theNormalizationFactor = 1.0 / sqrt(φ*φ + 1.0)
			//	inverse length of unnormalized vertex
			//	 = 1/√(φ² + 1) ≈ 0.52573111211913360603
		
		let a = 1.0 * theNormalizationFactor
		let b =  φ  * theNormalizationFactor

		//	The icosahedron's 12 vertices, normalized to unit length
		let v: [SIMD3<Double>] = [
		
			SIMD3<Double>( 0.0, -a,  -b  ),
			SIMD3<Double>( 0.0, +a,  -b  ),
			SIMD3<Double>( 0.0, -a,  +b  ),
			SIMD3<Double>( 0.0, +a,  +b  ),

			SIMD3<Double>( -a,  -b,  0.0 ),
			SIMD3<Double>( +a,  -b,  0.0 ),
			SIMD3<Double>( -a,  +b,  0.0 ),
			SIMD3<Double>( +a,  +b,  0.0 ),

			SIMD3<Double>( -b,  0.0, -a  ),
			SIMD3<Double>( -b,  0.0, +a  ),
			SIMD3<Double>( +b,  0.0, -a  ),
			SIMD3<Double>( +b,  0.0, +a  )
		]
		
		//	The icosahedron's 20 facets
		let f: [ UInt16 ] =
		[
			//	Winding order is clockwise when viewed
			//	from outside the icosahedron in a right-handed coordinate system.

			//	side-based faces

			 0,  1,  8,
			 1,  0, 10,
			 2,  3, 11,
			 3,  2,  9,

			 4,  5,  0,
			 5,  4,  2,
			 6,  7,  3,
			 7,  6,  1,

			 8,  9,  4,
			 9,  8,  6,
			10, 11,  7,
			11, 10,  5,

			//	corner-based faces
			
			 0,  8,  4,
			 2,  4,  9,
			 1,  6,  8,
			 3,  9,  6,
			 0,  5, 10,
			 2, 11,  5,
			 1, 10,  7,
			 3,  7, 11
		]

		let theIcosahedronMesh = Mesh3D(
									vertexBuffer: v,
									indexBuffer: f)
		
		return theIcosahedronMesh
	}

	private func unitOctahedronMesh() -> Mesh3D {
		
		//	The octahedron's 6 vertices sit at
		//
		//			(±1,  0,  0)
		//			( 0, ±1,  0)
		//			( 0,  0, ±1)
		//
		let v = [
		
			SIMD3<Double>(-1.0,  0.0,  0.0),
			SIMD3<Double>( 1.0,  0.0,  0.0),
			SIMD3<Double>( 0.0, -1.0,  0.0),
			SIMD3<Double>( 0.0,  1.0,  0.0),
			SIMD3<Double>( 0.0,  0.0, -1.0),
			SIMD3<Double>( 0.0,  0.0,  1.0)
		]

		//	The octahedron's 8 facets
		let f: [ UInt16 ] =
		[
			//	Winding order is clockwise when viewed
			//	from outside the octahedron in a right-handed coordinate system.

			5, 3, 1,
			4, 1, 3,
			4, 2, 1,
			5, 1, 2,
			4, 3, 0,
			4, 0, 2,
			5, 0, 3,
			5, 2, 0
		 ]

		let theOctahedronMesh = Mesh3D(
									vertexBuffer: v,
									indexBuffer: f)
		
		return theOctahedronMesh
	}

	private func subdivideUnitRadiusPolyhedronMesh(
		mesh: Mesh3D
	) -> Mesh3D {

		precondition(
			mesh.indexBuffer.count % 3 == 0,
			"The index count must be a multiple of 3")
		
		let theSrcVertexCount = mesh.vertexBuffer.count
		let theSrcFacetCount = mesh.indexBuffer.count / 3

		let theSrcVertices = mesh.vertexBuffer
		let theSrcIndices = mesh.indexBuffer	//	treat as bare indices, not as facets

		//	We'll subdivide the mesh, replacing each old facet with four new ones.
		//
		//		  /\		//
		//		 /__\		//
		//		/_\/_\		//
		//
		
		let theSubdivisionFacetCount = 4 * theSrcFacetCount
		let theSubdivisionIndexCount = 3 * theSubdivisionFacetCount
		
		//	Euler's formula says that for any cell decomposition of a sphere
		//
		//		v - e + f = 2
		//
		//	Our meshes are comprised exclusively of triangles, so
		//
		//		e = 3f/2
		//
		//	Combining those two formulas gives
		//
		//		v = f/2 + 2
		//
		//	The diagram above shows that
		//
		//		subdivisionF = 4 * originalF
		//
		//	so
		//
		//		subdivisionV = subdivisionF/2 + 2
		//					 = (4*originalF)/2 + 2
		//					 = 2*originalF + 2
		//
		let theSubdivisionVertexCount =  2 * theSrcFacetCount  +  2

		//	Allocate theSubdivisionMesh
		var theSubdivisionVertices = Array(
										repeating: SIMD3<Double>.zero,
										count: theSubdivisionVertexCount)
		var theSubdivisionFacetIndices = Array(
											repeating: UInt16(0),
											count: theSubdivisionIndexCount)

		//	Keep a running count of the number of vertices
		//	we've written into theSubdivisionVertices,
		//	so we always know where to write the next one.
		//
		var theVertexCount = 0

		//	First copy the source mesh vertex positions to the destination mesh...
		for i in 0..<theSrcVertexCount {
			theSubdivisionVertices[i] = theSrcVertices[i]
		}
		theVertexCount += theSrcVertexCount

		//	...and then create one new vertex on each edge.
		//
		//	Use theNewVertexTable[][] to index them,
		//	so two facets sharing an edge can share the same vertex.
		//
		//	theNewVertexTable[v0][v1] takes the indices v0 and v1 of two old vertices,
		//	and gives the index of the new vertex that sits at the midpoint of the edge
		//	that connects v0 and v1.
		//
		//	The size of theNewVertexTable grows as the square of the number of source vertices.
		//	For modest meshes this won't be a problem. For larger meshes a fancier algorithm,
		//	with linear rather than quadratic memory demands, could be used.
		//
		//	Initialize all of theNewVertexTable[][]'s entries to
		//
		//		notPresent = Int.max
		//
		//	to indicate that no new vertices have yet been created.
		//	This seems easier than working with optional values
		//	(which would be of type "Int?").
		//
		let notPresent = Int.max
		var theNewVertexTable: [ [ Int ] ]
			= Array(
				repeating: Array(
						repeating: notPresent,
						count: theSrcVertexCount),
				count: theSrcVertexCount
			)

		//	For each edge in the source mesh, create a new vertex at its midpoint.
		for i in 0..<theSrcFacetCount {
			for j in 0..<3 {
			
				let v0 = Int(theSrcIndices[3*i +    j   ])	//	Cast from UInt16 to Int
				let v1 = Int(theSrcIndices[3*i + (j+1)%3])
				
				if theNewVertexTable[v0][v1] == notPresent {
				
					precondition(
						theVertexCount < theSubdivisionVertexCount,
						"Internal error #1 in subdivideUnitRadiusPolyhedronMesh()")

					//	The new vertex will sit midway between vertices v0 and v1.
					theSubdivisionVertices[theVertexCount]
						= normalize(Double(0.5) * ( theSubdivisionVertices[v0]
												  + theSubdivisionVertices[v1]))
						//	normalize() should be safe here, because
						//		we'll never pass it a zero vector

					//	Record the new vertex at [v1][v0] as well as [v0][v1],
					//	so only one new vertex will get created for each edge.
					theNewVertexTable[v0][v1] = theVertexCount
					theNewVertexTable[v1][v0] = theVertexCount

					theVertexCount += 1
				}
			}
		}
		precondition(
			theVertexCount == theSubdivisionVertexCount,
			"Internal error #2 in subdivideUnitRadiusPolyhedronMesh()")

		//	For each facet in the source mesh,
		//	create four smaller facets in the subdivision.
		for i in 0 ..< theSrcFacetCount {
		
			//	Call the old vertices incident to this facet v0, v1 and v2
			let v0 = Int(theSrcIndices[3*i + 0])	//	Cast from UInt16 to Int
			let v1 = Int(theSrcIndices[3*i + 1])
			let v2 = Int(theSrcIndices[3*i + 2])

			//	Call the new vertices -- which sit at the midpoints
			//	of the old edges -- vv0, vv1, and vv2.
			//	Each vv_j sits opposite the corresponding v_j .
			let vv0 = theNewVertexTable[v1][v2]
			let vv1 = theNewVertexTable[v2][v0]
			let vv2 = theNewVertexTable[v0][v1]

			//	Create the new facets.

			theSubdivisionFacetIndices[ 3*(4*i + 0) + 0 ] = UInt16(vv0)
			theSubdivisionFacetIndices[ 3*(4*i + 0) + 1 ] = UInt16(vv1)
			theSubdivisionFacetIndices[ 3*(4*i + 0) + 2 ] = UInt16(vv2)

			theSubdivisionFacetIndices[ 3*(4*i + 1) + 0 ] = UInt16(v0)
			theSubdivisionFacetIndices[ 3*(4*i + 1) + 1 ] = UInt16(vv2)
			theSubdivisionFacetIndices[ 3*(4*i + 1) + 2 ] = UInt16(vv1)

			theSubdivisionFacetIndices[ 3*(4*i + 2) + 0 ] = UInt16(v1)
			theSubdivisionFacetIndices[ 3*(4*i + 2) + 1 ] = UInt16(vv0)
			theSubdivisionFacetIndices[ 3*(4*i + 2) + 2 ] = UInt16(vv2)

			theSubdivisionFacetIndices[ 3*(4*i + 3) + 0 ] = UInt16(v2)
			theSubdivisionFacetIndices[ 3*(4*i + 3) + 1 ] = UInt16(vv1)
			theSubdivisionFacetIndices[ 3*(4*i + 3) + 2 ] = UInt16(vv0)
		}
		

		let theSubdivisionMesh = Mesh3D(
									vertexBuffer: theSubdivisionVertices,
									indexBuffer: theSubdivisionFacetIndices)
		
		return theSubdivisionMesh
	}


// MARK: -
// MARK: Cue stick

	func makeCueStickMesh(
		device: MTLDevice,
		buttRadius: Double = cueStickButtRadius,
		tipRadius: Double = cueStickTipRadius,
		length: Double = cueStickLength
	) -> NEBMesh {

		//	We subdivide the mesh not only in the meridional direction
		//	but in the longitudinal direction as well, so that players
		//	may see how the cue stick seems to bend towards or away from
		//	them. This is definitely an effect when the scene is rendered
		//	in the azimuthal equidistant map to simulate native perception
		//	of distances, but it's slightly noticeable in un-adjusted
		//	rendering or not, especially in spherical geometry
		//	with a cue stick that's at least 1 meter long.

		let nL = 32	//	number of divisons in longitudinal direction
		let nM = 16	//	number of divisons in  meridional  direction

		//	Note:  Going around in the meridional direction,
		//	the first and last vertices will share the same position
		//	but with different texture coordinates.
		
		let theVertexCount = (nL + 1) * (nM + 1)
		
		let theSideFacetCount = 2 * nL * nM
		let theButtFacetCount = nM - 2
		let theTipFacetCount  = nM - 2
		let theFacetCount     = theSideFacetCount + theButtFacetCount + theTipFacetCount
		
		let theIndexCount = 3 * theFacetCount

		//	Define the vertices.
		//
		//	Place the cue stick's tip at the north pole,
		//	with the stick extending in the positive z direction
		//	(towards the observer).
		
		//	This code gets run only once, when we create the cue stick mesh.
		//	So I've written it to favor clarity over runtime efficiency.

		let ρ = spaceGeometry.spaceRadius

		//	Define the vertex positions.
		var thePositions = Array(repeating: SIMD4<Double>.zero, count: theVertexCount)
		var theTexCoords = Array(repeating: SIMD3<Double>.zero, count: theVertexCount)
		for ell in 0 ... nL {
		
			//	Define a meridional ring of vertices wrapping once around the cue stick.
			
			//	What fraction of way from the tip to the butt is this meridional ring?
			let s = Double(ell) / Double(nL)	//	0.0 ≤ s ≤ 1.0
			
			//	How far, in radians, is this ring from the cue stick's tip?
			let φ = (s * length) / ρ
			
			//	The cue stick's radius increases linearly from tipRadius to buttRadius.
			//	Divide by the spaceRadius to convert it to radians.
			let r = ( (1.0 - s)*tipRadius + s*buttRadius ) / ρ
			
			//	To distribute the curvature more evenly,
			//	let the odd-numbered rings be a half-step out of phase
			//	with the even-numbered rings.
			let thePhaseShift = (ell % 2 == 1 ? 0.5 : 0.0)
				
			let theRadialTranslation = axisAlignedTranslation(
										axis: 0,
										distanceInRadians: r,
										geometryType: spaceGeometry.geometryType)
		
			let theLongitudinalTranslation = axisAlignedTranslation(
												axis: 2,
												distanceInRadians: φ,
												geometryType: spaceGeometry.geometryType)
		
			let theBasePosition = theRadialTranslation * SIMD4<Double>(0.0, 0.0, 0.0,  ρ )

			for m in 0 ... nM {

				//	What fraction of way around the meridional ring is this vertex?
				//
				//		Note:  t may exceed 1.0 because of the phase shift.
				//		Be sure to set the MTLSamplerAddressMode to .repeat
				//		in the meridional direction (but .clamp in the
				//		longitudinal direction)
				//
				let t = (Double(m) + thePhaseShift) / Double(nM)

				let θ = 2.0 * Double.pi * t
				
				let theMeridionalRotation = axisAlignedRotation(
												from: 0,	//	zero, not theta
												towards: 1,
												angle: θ)	//	theta, not zero

				let theRotatedPosition = theMeridionalRotation * theBasePosition

				let theTranslatedPosition = theLongitudinalTranslation * theRotatedPosition

				thePositions[(nM + 1)*ell + m] = theTranslatedPosition
				theTexCoords[(nM + 1)*ell + m] = SIMD3<Double>(s, t, 0)
			}
		}
		
		//	Define the mesh facets.
		//
		//	Let the facets wind clockwise when viewed from the outside
		//	in a right-handed coordinate system.
		
		var theFacetIndices = Array(repeating: 0, count: theIndexCount)
		var c = 0	//	running count of number of facets already defined
		
		//	Let's define the facets in more-or-less the same
		//	order that the vertices appear in the vertex array,
		//	just in case that makes it easier for the GPU
		//	to re-use transformed vertices.

		//	Tip facets
		let theTipBaseVertex = 0
		for m in 1 ... nM - 2 {
		
			theFacetIndices[3*c + 0] = theTipBaseVertex +  0
			theFacetIndices[3*c + 1] = theTipBaseVertex +  m
			theFacetIndices[3*c + 2] = theTipBaseVertex + m+1
			c += 1
		}
		
		//	Side facets
		for ell in 0 ..< nL {	//	excludes nL
			for m in 0 ..< nM {	//	excludes nM
				
				theFacetIndices[3*c + 0] = (nM + 1) *   ell    +    m
				theFacetIndices[3*c + 1] = (nM + 1) * (ell+1)  +    m
				theFacetIndices[3*c + 2] = (nM + 1) *   ell    +  (m+1)
				c += 1

				theFacetIndices[3*c + 0] = (nM + 1) *   ell    +  (m+1)
				theFacetIndices[3*c + 1] = (nM + 1) * (ell+1)  +    m
				theFacetIndices[3*c + 2] = (nM + 1) * (ell+1)  +  (m+1)
				c += 1
			}
		}

		//	Butt facets
		let theButtBaseVertex = (nM + 1) * nL
		for m in 1 ... nM - 2 {
		
			theFacetIndices[3*c + 0] = theButtBaseVertex +  0
			theFacetIndices[3*c + 1] = theButtBaseVertex + m+1
			theFacetIndices[3*c + 2] = theButtBaseVertex +  m
			c += 1
		}
		
		precondition(
			c == theFacetCount,
			"Internal error: bad facet count in makeCueStickMesh()")

		//	Package up the mesh.
		
		let theMesh = NEBMesh(
						device: device,
						vertexCount: theVertexCount,
						indexCount: theIndexCount)
		
		precondition(
			theMesh.vertexCount == theVertexCount,
			"Wrong vertex count in makeCueStickMesh()")
		let theVertices = theMesh.vertexBuffer.contents()
			.bindMemory(to: NEBVertex.self, capacity: theVertexCount)
		for i in 0 ..< theVertexCount {
		
			theVertices[i] = NEBVertex(
				pos: SIMD4<Float32>(thePositions[i]),
				tex: SIMD3<Float32>(theTexCoords[i]) )
		}

		precondition(
			theMesh.indexCount == theIndexCount,
			"Wrong index count in makeCueStickMesh()")
		let theIndices = theMesh.indexBuffer.contents()
			.bindMemory(to: UInt16.self, capacity: theIndexCount)
		for i in 0 ..< theIndexCount {
			theIndices[i] = UInt16(theFacetIndices[i])
		}

		return theMesh
	}
}

extension VirtualHand {
	
	var cueStickGripPoint: SIMD4<Double> {

		let theGripPointInHandCoordinates
			= joints.cueStickGripPointInHandCoordinates()
			
		let theGripPointInGameCoordinates
			= pose * theGripPointInHandCoordinates
		
		return theGripPointInGameCoordinates
	}
	
	var cueStickBridgePoint: SIMD4<Double> {

		let theBridgePointInHandCoordinates
			= joints.cueStickBridgePointInHandCoordinates()
			
		let theBridgePointInGameCoordinates
			= pose * theBridgePointInHandCoordinates
		
		return theBridgePointInGameCoordinates
	}
	
	var cueBallInHandPoint: SIMD4<Double> {	//	in game coordinates

		let theBallInHandPointInHandCoordinates
			= joints.cueBallInHandPointInHandCoordinates()
			
		let theBallInHandPointInGameCoordinates
			= pose * theBallInHandPointInHandCoordinates
		
		return theBallInHandPointInGameCoordinates
	}
}
