Nested Angular Grids (for adaptive quadrature)

We may need to manually generated nested grids. One can generate icosahedral nested grids by subdividing an icosahedral triangulation. The sequence would be $12, 48=4 \times 12, 192 = 4 \times 48, 768 = 4 \times 192, 3072 = 4 \times 768, \ldots$. Once the sequence is found, one has to solve for the weights that ensure that as many spherical harmonics as possible are integrated correctly. There are $(\ell+1)^2$ spherical harmonics with angular momentum less than or equal to $\ell$, so (at a minimum) the aforementioned grids have degree $3,6,13,27,55,\ldots$. In general, due to symmetry, one should be able to do a bit better than this.

I asked ChatGPT to make a script to generate the nested grid points, and got the following (not checked); the example at the end is for 192 points. The angular degrees would probably be updated based on these grids. We'd have angular grids:
(48, 192, 768) fine
(192, 768, 3072) crazy big
We'd need to figure out which radial grids are a good match for these angular grids.

import numpy as np

def icosahedron_vertices():
    # Golden ratio
    phi = (1 + np.sqrt(5)) / 2

    # Normalize the length to 1
    a = 1 / np.sqrt(1 + phi**2)
    b = phi * a

    # Vertices of an icosahedron
    vertices = [
        (-a,  b,  0), ( a,  b,  0), (-a, -b,  0), ( a, -b,  0),
        ( 0, -a,  b), ( 0,  a,  b), ( 0, -a, -b), ( 0,  a, -b),
        ( b,  0, -a), ( b,  0,  a), (-b,  0, -a), (-b,  0,  a)
    ]
    return np.array(vertices)

def subdivide(vertices, faces):
    # New vertex set, including old vertices
    new_vertices = vertices.tolist()
    vertex_map = {}

    # Helper function to find the midpoint and index
    def midpoint_index(v1, v2):
        # Check if we have already calculated this
        edge = tuple(sorted([v1, v2]))
        if edge in vertex_map:
            return vertex_map[edge]

        # Calculate midpoint and normalize it to the sphere
        midpoint = (vertices[v1] + vertices[v2]) / 2
        midpoint /= np.linalg.norm(midpoint)
        
        # Store the new vertex index
        index = len(new_vertices)
        new_vertices.append(midpoint)
        vertex_map[edge] = index
        return index

    new_faces = []

    for v1, v2, v3 in faces:
        # Get the midpoints
        a = midpoint_index(v1, v2)
        b = midpoint_index(v2, v3)
        c = midpoint_index(v3, v1)

        # Create four new faces
        new_faces.extend([
            [v1, a, c],
            [v2, b, a],
            [v3, c, b],
            [a, b, c]
        ])

    return np.array(new_vertices), np.array(new_faces)

def generate_icosphere(subdivisions):
    # Initial icosahedron vertices and faces
    vertices = icosahedron_vertices()
    faces = [
        [0, 11, 5], [0, 5, 1], [0, 1, 7], [0, 7, 10], [0, 10, 11],
        [1, 5, 9], [5, 11, 4], [11, 10, 2], [10, 7, 6], [7, 1, 8],
        [3, 9, 4], [3, 4, 2], [3, 2, 6], [3, 6, 8], [3, 8, 9],
        [4, 9, 5], [2, 4, 11], [6, 2, 10], [8, 6, 7], [9, 8, 1]
    ]

    # Subdivide the triangles the number of times specified
    for _ in range(subdivisions):
        vertices, faces = subdivide(vertices, faces)

    return vertices

# Generate vertices for a subdivided icosahedron with 2 subdivisions
icosphere_vertices = generate_icosphere(2)
icosphere_vertices