from bpy import context, ops
from math import cos, sin
from mathutils import Vector
from random import TWOPI
import cmath

def stereo(z,w):
    "Return the stereographic projection of a point (z,w) in C^2"
    return Vector((z.real / (1-w.imag), z.imag / (1-w.imag), w.real / (1-w.imag)))

def add_curve(startpt):
    ops.curve.primitive_bezier_circle_add(enter_editmode=True)
    curve = context.active_object
    ops.curve.subdivide(number_cuts=36)
    bez_points = curve.data.splines[0].bezier_points
    sz = len(bez_points)
    i_to_theta = TWOPI / sz
    j = complex(0,1)
    for i in range(0, sz, 1):
        newz = startpt[0]*cmath.exp(i * i_to_theta*j)
        neww = startpt[1]*cmath.exp(i * i_to_theta*j)
        bez_points[i].co = stereo(newz,neww)
    ops.object.mode_set(mode='OBJECT')
    return curve

phi = TWOPI/4+.01
N=33
curves = list()
for k in range(N-3):
    theta = TWOPI*k/N
    for i in range(k,k+1):
        startpt = (complex(sin(theta)*sin(phi)*cos(TWOPI*i/(4*N)),cos(theta)*sin(phi)*cos(TWOPI*i/(4*N))),complex(sin(phi)*sin(TWOPI*i/(4*N)),cos(phi)))
        (x,y,z) = stereo(startpt[0],startpt[1])        
        curve = add_curve(startpt)
        curves.append(curve)