from openscad import *
from math import *

def mean(p1,p2,r):
    pm = [(p1[i]+p2[i])/2.0 for i in range(3) ]
    f=r/sqrt(pm[0]*pm[0]+pm[1]*pm[1]+pm[2]*pm[2])
    return [pm[i]*f for i in range(3) ]
   
def geodesic(r,n):
    if n == 0:
        yield [[ 0,-r,0],[r, 0, 0],[0,0, r]]
        yield [[ r, 0,0],[0, r, 0],[0,0, r]]
        yield [[ 0, r,0],[-r, 0,0],[0,0, r]]
        yield [[-r, 0,0],[0, -r,0],[0,0, r]]
        yield [[ 0, r,0],[r, 0, 0],[0,0,-r]]
        yield [[-r, 0,0],[0, r, 0],[0,0,-r]]
        yield [[ 0,-r,0],[-r, 0,0],[0,0,-r]]
        yield [[ r, 0,0],[0, -r,0],[0,0,-r]]
    else:
        for p1,p2,p3 in geodesic(r,n-1):
            p12=mean(p1,p2,r)
            p23=mean(p2,p3,r)
            p31=mean(p3,p1,r)
            yield [p1, p12, p31]
            yield [p2, p23, p12]
            yield [p3, p31, p23]
            yield [p12,p23, p31]
        
vertices = []
triangles = []
cols = []
ptmap = {}
def find_index(pt):
    key="%g/%g/%g"%(pt[0],pt[1],pt[2])
    if not key in ptmap:
        ptmap[key]=len(vertices)        
        vertices.append(pt)            
    return ptmap[key]
    
r=10    
for tri in geodesic(r,3):
    triangles.append([find_index(tri[0]),find_index(tri[1]),find_index(tri[2]) ])
    cols.append([abs(tri[0][i]+tri[1][i]+tri[2][i])/(3*r) for i in range(3) ])
    
polyhedron(vertices, triangles, colors=cols).show()