import numpy

offset = numpy.zeros([29,3])
zeta = numpy.pi*numpy.array([0.25, 0.75, 1.25, 1.75])
theta = numpy.pi*numpy.array([0.25, 0.75])

r1 = numpy.cbrt(1.0/29.0)
r2 = numpy.cbrt(15.0/29.0)
r = numpy.array([0.75*(r2*r2*r2*r2 - r1*r1*r1*r1)/(r2*r2*r2-r1*r1*r1), 0.75*(1.0-r2*r2*r2*r2)/(1.0-r2*r2*r2)])

m = 1
for ir in range(2):
  for z in zeta:
    for t in theta:
      offset[m][0] = r[ir]*numpy.sin(t)*numpy.cos(z)
      offset[m][1] = r[ir]*numpy.sin(t)*numpy.sin(z)
      offset[m][2] = r[ir]*numpy.cos(t)
      m+=1
  for j in range(-1,3,2):
    offset[m][0] = r[ir]*float(j)
    offset[m][1] = 0.0
    offset[m][2] = 0.0
    m+=1
    offset[m][0] = 0.0
    offset[m][1] = r[ir]*float(j)
    offset[m][2] = 0.0
    m+=1
    offset[m][0] = 0.0
    offset[m][1] = 0.0
    offset[m][2] = r[ir]*float(j)
    m+=1

print(offset)