import numpy as np

# Number of circumferential faces (should be even number)
n_faces = 40

# Outer Radius 
ro = 27. # mm

# Width
w = 4.2 # mm

# Ratio
ratio = 2/(1+np.sqrt(5))# 0 ~ 1
ratio = np.max([0., ratio])
ratio = np.min([1., ratio])

offset = 2 * ratio - 1.

pdo =  self.GetPolyDataOutput()

dtheta = 2 * np.pi / n_faces

multiplier = np.arange(n_faces, dtype=np.double)
multiplier[1::2] = multiplier[1::2] + offset

theta = multiplier * dtheta

x = ro * np.sin(theta)
y = ro * np.cos(theta)

torsion = dtheta + 0.2*dtheta

sin_dtheta = np.sin(torsion)
cos_dtheta = np.cos(torsion)

# rotation_mat = np.array([cos_dtheta, - sin_dtheta, sin_dtheta, cos_dtheta]).reshape((2,2))

x_ = cos_dtheta * x - sin_dtheta * y
y_ = sin_dtheta * x + cos_dtheta * y

x1 = np.concatenate((x_[1:], x_[:1]))
y1 = np.concatenate((y_[1:], y_[:1]))

nodes = vtk.vtkPoints()
for i, (xi, yi) in enumerate(zip(x, y)):
    nodes.InsertNextPoint(xi, yi, 0.0)

for i, (xi, yi) in enumerate(zip(x1, y1)):
    nodes.InsertNextPoint(xi, yi, w)

elements = vtk.vtkCellArray()
for i in range(n_faces):

    # Create the object that stores the list of nodeIDs
    elementIdList0 = vtk.vtkIdList()

    # Insert the nodeIDs in the appropriate order.
    i_next = (i + 1) % n_faces
    elementIdList0.InsertNextId(i)
    elementIdList0.InsertNextId(i+n_faces)
    elementIdList0.InsertNextId(i_next+n_faces)
    elementIdList0.InsertNextId(i_next)

    elements.InsertNextCell(elementIdList0)

# Specify that the variable 'nodes' as the points in the polydata
pdo.SetPoints(nodes)
# Similarly, specify that the variable 'elements' as the polys in polydata
pdo.SetPolys(elements)