#!/usr/bin/env python
import sys, math


simple = True

points = {}
triangles = []

schwarz = "0x000000"
rot     = "0x0000ff"
gruen   = "0x00ff00"
gelb    = "0x0000ff"
blau    = "0x00ffff"

minc = None
maxc = None

class Vect:
   def __init__ (self, x, y, z):
      self.x = x
      self.y = y
      self.z = z


   def __add__ (self, other):
      return Vect (self.x + other.x,
                   self.y + other.y,
                   self.z + other.z)

   def __sub__ (self, other):
      return Vect (self.x - other.x,
                   self.y - other.y,
                   self.z - other.z)

   def __neg__ (self):
      return Vect (-self.x, -self.y, -self.z)

   def __mul__ (self, other):
      f = float (other)
      return Vect (f*self.x, f*self.y, f*self.z)

   def __div__ (self, other):
      f = float (other)
      return Vect (self.x/f, self.y/f, self.z/f)

   def __rmul__ (self, other):
      f = float (other)
      return Vect (f*self.x, f*self.y, f*self.z)

   def dot (self, other):
      return self.x * other.x + self.y * other.y + self.z * other.z

   def cross (self, other):
      return Vect (self.y * other.z - self.z * other.y,
                   self.z * other.x - self.x * other.z,
                   self.x * other.y - self.y * other.x)

   def len (self):
      return (self.x ** 2 + self.y ** 2 + self.z ** 2) ** 0.5

   def __repr__ (self):
      return "Vect (%f, %f, %f)" % (self.x, self.y, self.z)



class Quat:
   def __init__ (self, w, x, y, z):
      self.w = w
      self.x = x
      self.y = y
      self.z = z

   def __add__ (self, other):
      return Quat (self.w + other.w,
                   self.x + other.x,
                   self.y + other.y,
                   self.z + other.z)

   def qmult (s, o):
      return Quat (s.w * o.w - s.x * o.x - s.y * o.y - s.z * o.z,
                   s.w * o.x + s.x * o.w + s.y * o.z - s.z * o.y,
                   s.w * o.y - s.x * o.z + s.y * o.w + s.z * o.x,
                   s.w * o.z + s.x * o.y - s.y * o.x + s.z * o.w)


   def __mul__ (self, other):
      f = float (other)
      return Quat (f*self.w, f*self.x, f*self.y, f*self.z)

   def __neg__ (self):
      return Quat (-self.w, -self.x, -self.y, -self.z)

   def __invert__ (self):
      return Quat (self.w, -self.x, -self.y, -self.z)

   def __repr__ (self):
      return "Quat (%f, %f, %f, %f)" % (self.w, self.x, self.y, self.z)

   def len (self):
      return (self.w ** 2 + self.x ** 2 + self.y ** 2 + self.z ** 2) ** 0.5



def triangle_add (v1, v2, v3, color=None):
   global minc, maxc

   for x in [1, -1]:
      for y in [1, -1]:
         i = []
         if x * y > 0:
            vv1 = Vect (v1.x * x, v1.y * y, v1.z)
            vv2 = Vect (v2.x * x, v2.y * y, v2.z)
            vv3 = Vect (v3.x * x, v3.y * y, v3.z)
         else:
            vv1 = Vect (v1.x * x, v1.y * y, v1.z)
            vv2 = Vect (v3.x * x, v3.y * y, v3.z)
            vv3 = Vect (v2.x * x, v2.y * y, v2.z)
         for v in [vv1, vv2, vv3]:
            if minc == None:
               minc = [v.x, v.y, v.z]
               maxc = [v.x, v.y, v.z]
            else:
               minc = [min (minc[0], v.x),
                       min (minc[1], v.y),
                       min (minc[2], v.z)]
               maxc = [max (maxc[0], v.x),
                       max (maxc[1], v.y),
                       max (maxc[2], v.z)]
            k = "%10.5f %10.5f %10.5f" % (v.x, v.y, v.z)
            if not points.has_key (k):
               points[k] = len (points)
            i.append (points[k])

         n = (vv2 - vv1).cross (vv3 - vv1)
         l = n.len ()
         if l > 0.0001:
            n /= n.len ()
            if color:
               c = color
            elif n.z > 0:
               c = gruen
            else:
               c = rot

            triangles.append (tuple (i) + (c,))


def mkrot (angle, vec):
   sa = math.sin (angle/2)
   ca = math.cos (angle/2)
   l = (vec.x**2 + vec.y**2 + vec.z**2)**0.5

   return Quat (ca, sa*vec.x/l, sa*vec.y/l, sa*vec.z/l)


def find_center (v0, v1, v2):
   d1 = (v1 - v0) * 0.5
   d2 = (v2 - v0) * 0.5
   n = d1.cross (d2)

   x1 = v0 + d1
   x3 = v0 + d2

   c1 = n.cross (d1)
   c2 = n.cross (d2)

   x2 = x1 + c1
   x4 = x3 + c2

   a = x2 - x1
   b = x4 - x3
   c = x3 - x1

   s = (c.cross (b)).dot (a.cross (b)) / (a.cross(b)).dot (a.cross(b))

   return x1 + s * a


def interpol_circle (v1, v2, v3, steps=10):
   l = []
   c = find_center (v1, v2, v3)
   v1 = v1 - c
   v2 = v2 - c
   v3 = v3 - c

   angle = math.acos (v1.dot (v2) / v1.len () / v2.len ())
   n = v1.cross (v2)
   for i in range (steps):
      rot = mkrot (i * angle / (steps - 1), n)
      qv = Quat (0, v1.x, v1.y, v1.z)
      rv = rot.qmult (qv).qmult (~rot)
      l.append (Vect (rv.x, rv.y, rv.z) + c)

   return c, l


def gen_sixth (r = 20):
   a = (1./3) ** 0.5
   b = (2./3) ** 0.5
   w2 = 2**0.5

   v1 = Vect (a-a/w2 , 0, -a-a/w2)
   v2 = Vect (a, a, -a)
   v3 = Vect (0, 0, -1)
   v4 = Vect (0, a+a/w2, -a+a/w2)
  
   steps = 15
   
   # triangle_add (v2, Vect (a, -a, -a), v1)
   # triangle_add (v3, Vect (0, -a-a/w2, -a+a/w2), v4)
   c1, l1 = interpol_circle (v2, v4, Vect (-v2.x, v2.y, v2.z), steps)
   c2, l2 = interpol_circle (v1, v3, Vect (-v1.x, v1.y, v1.z), steps)

   face = []
   for i in range (len (l1)):
      c, l = interpol_circle (l1[i], l2[i],
                              Vect(l1[i].x, -l1[i].y, l1[i].z), steps)
      face.append (l)

   R = 30.0
   Ra = R * a

   if simple:
      Ri = R * 0.4
      Va = Ra
   else:
      Ri = R * 0.7
      Va = Ra - 0.2

   for j in range (1, len (face)):
      for i in range (1, len (face[j])):
         if i % 8 in [1, 4, 5, 0] and j % 8 in [1, 4, 5, 0]:
            c = schwarz
         elif j % 8 in [2, 3, 6, 7]:
            c = blau
         else:
            c = gelb
         if j % 8 in [2, 3] and i % 8 in [2, 3]:
            c = gelb
         if j % 8 in [6, 7] and i % 8 in [6, 7]:
            c = gelb
         triangle_add (face[j-1][i] * R, face[j-1][i-1] * R, face[j][i] * R, c)
         triangle_add (face[j-1][i-1] * R, face[j][i-1] * R, face[j][i] * R, c)


   # Fase Seitenflaeche
   l = face[0]
   v0 = Vect (l[0].x * R - (Ra - Va), l[0].y * R - (Ra - Va), l[0].z * R)
   v1 = Vect (l[-1].x * R, l[-1].y * R, l[-1].z * R + (Ra - Va))
   c1, l1 = interpol_circle (v0, v1, Vect (v0.x, -v0.y, v0.z), steps)
   for i in range (1, len (l)):
      triangle_add (l[i-1]*R, l[i]*R, l1[i-1], schwarz)
      triangle_add (l1[i-1], l[i]*R, l1[i], schwarz)

   # Seitenflaeche
      f0 = float (i-1) / (steps - 1)
      f1 = float (i) / (steps - 1)

      v2 = Vect (Va,  0, -Ra) * f0 + Vect (Va,  Va, -Ra) * (1. - f0)
      v3 = Vect (Va,  0, -Ra) * f1 + Vect (Va,  Va, -Ra) * (1. - f1)

      triangle_add (l1[i-1], l1[i], v3, gruen)
      triangle_add (l1[i-1], v3, v2, gruen)


   # Uebergang Innenflaeche
   l = [x[0] for x in face]

   # Innenflaeche
   for i in range (1, steps):
      f0 = float (i-1) / (steps - 1)
      f1 = float (i) / (steps - 1)

      v0 = R * l[i-1]
      v1 = R * l[i]
      v2 = Vect ( 0, Ra, -Ra) * f0 + Vect (Ra,  Ra, -Ra) * (1. - f0)
      v3 = Vect ( 0, Ra, -Ra) * f1 + Vect (Ra,  Ra, -Ra) * (1. - f1)

      triangle_add (v0, v3, v1, gruen)
      triangle_add (v0, v2, v3, gruen)

      v0 = v2 / Ra * Va
      v1 = v3 / Ra * Va
      v0.z = v1.z = -Ra
      triangle_add (v0, v1, v3, gruen)
      triangle_add (v0, v3, v2, gruen)



   # innere Kugelflaeche

   for x in range (1, steps):
      for y in range (1, steps):
         x0 = Va * (x-1) / (steps - 1)
         y0 = Va * (y-1) / (steps - 1)
         x1 = Va * x / (steps - 1)
         y1 = Va * y / (steps - 1)

         v0 = Vect (x0, y0, -Ra)
         v1 = Vect (x0, y1, -Ra)
         v2 = Vect (x1, y0, -Ra)
         v3 = Vect (x1, y1, -Ra)

         if simple:
            v0 = v0 / Ra * Ri
            v1 = v1 / Ra * Ri
            v2 = v2 / Ra * Ri
            v3 = v3 / Ra * Ri
         else:
            v0 = v0 / v0.len () * Ri
            v1 = v1 / v1.len () * Ri
            v2 = v2 / v2.len () * Ri
            v3 = v3 / v3.len () * Ri
               
         triangle_add (v0, v1, v3, rot)
         triangle_add (v0, v3, v2, rot)

   # innere Seitenwaende

   for i in range (1, steps):
      x0 = Va * (i-1) / (steps - 1)
      x1 = Va * i / (steps - 1)

      for v0, v1 in [(Vect (x0, Va, -Ra), Vect (x1, Va, -Ra)),
                     (Vect (Va, x1, -Ra), Vect (Va, x0, -Ra))]:
         if simple:
            v2 = v0 / Ra * Ri
            v3 = v1 / Ra * Ri
         else:
            v2 = v0 / v0.len () * Ri
            v3 = v1 / v1.len () * Ri

         triangle_add (v0, v1, v3, gruen)
         triangle_add (v0, v3, v2, gruen)


def gen_obj ():
   print "3DG1"
   print len (points)
   l = [(v, k) for k, v in points.items()]
   l.sort ()
   for i in l:
      print i[1]

   for i in triangles:
      print "3 %d %d %d %s" % i

def gen_stl ():
   print "solid puzzle"
   l = [(v, k) for k, v in points.items()]
   l.sort ()

   for j in triangles:
      i = [l[j[0]], l[j[1]], l[j[2]]]
      print "  facet normal 0 0 0"
      print "    outer loop"
      print "      vertex %s" % i[0][1]
      print "      vertex %s" % i[1][1]
      print "      vertex %s" % i[2][1]
      print "    endloop"
      print "  endfacet"
      
   print "endsolid puzzle"


if __name__ == '__main__':
   gen_sixth ()
   
   gen_stl ()

   print >> sys.stderr, minc
   print >> sys.stderr, maxc