hemisphere visible faces only

#!BPY
import bpy
import math

#############
n = 5
#############
m = bpy.data.meshes.new('hemisphereMesh0')
o = bpy.data.objects.new('hemisphere0', m)
bpy.context.scene.objects.link(o)

def vert():
    yield (0.0, 0.0, 1.0)
    for k in range(1, n+1):
        theta = float(k)*math.pi/2.0/float(n)
        phiDelta = 2.0*math.pi/6.0/float(k)
        z = math.cos(theta)
        r = math.sin(theta)
        for lm in range(6*k):
            phi = phiDelta*float(lm)
            x = r*math.cos(phi)
            y = r*math.sin(phi)
            yield(x,y,z)
        
def edge():
#    for k in range(3*n*(n+1)):
#        yield(k, k+1)
    for p in range(n):
        if p == 0:
            u0 = 0
        else:
            u0 = 1+3*p*(p-1)
        l0 = 1+3*p*(p+1)
        for q in range(6):
            u = u0+p*q
            l = l0+(p+1)*q
            yield(u,l)
            yield(l,u)
            for r in range(p):
                uu = u+r
                ll = l+r
                yield(uu,ll+1)
                if r==p-1 and q==5:
                    yield(ll+1, u0)
                else:
                    yield(ll+1,uu+1)
            yield(u,l)

def face():
    for q in range(1,6):
        yield(q,0,q+1)
    yield(6,0,1)
    for p in range(1,n):
        u00 = 3*p*(p-1)+1
        l00 = 3*p*(p+1)+1
        u0 = u00
        l0 = l00
        for q in range(6):
            for r in range(p):
                yield(l0+r, u0+r, l0+r+1)
                if q==5 and r==p-1:
                    yield(u0+r, l0+r+1, u00)
                else:
                    yield(u0+r, l0+r+1, u0+r+1)                    
            if q<5:
                yield(l0+r+1, u0+r+1, l0+r+2)
            else:
                yield(l0+r+1, u00, l00)
            u0 += p
            l0 += p+1

#m.from_pydata(list(vert()), list(edge()), [])
list(face())
m.from_pydata(list(vert()), [], list(face()))
m.update()