#!/usr/bin/env python

solid_name  = 'ASCII_STL_of_a_sphericon_by_CMG_Lee'
n_face_side = 60
r           = 1000

import re, io, math
def fmt(string): ## string.format(**vars()) using tags {expression!format} by CMG Lee
 def f(tag): i_sep = tag.rfind('!'); return (re.sub('\.0+$', '', str(eval(tag[1:-1])))
  if (i_sep < 0) else ('{:%s}' % tag[i_sep + 1:-1]).format(eval(tag[1:i_sep])))
 return (re.sub(r'(?<!{){[^{}]+}', lambda m:f(m.group()), string)
         .replace('{{', '{').replace('}}', '}'))
def append(obj, string): return obj.append(fmt(string))
def tabbify(cellss, separator='|'):
 cellpadss = [list(rows) + [''] * (len(max(cellss, key=len)) - len(rows)) for rows in cellss]
 fmts = ['%%%ds' % (max([len(str(cell)) for cell in cols])) for cols in zip(*cellpadss)]
 return '\n'.join([separator.join(fmts) % tuple(rows) for rows in cellpadss])
def roundm(x, multiple=1):
 if   (isinstance(x, tuple)): return tuple(roundm(list(x), multiple))
 elif (isinstance(x, list )): return [roundm(x_i, multiple) for x_i in x]
 else: return int(math.floor(float(x) / multiple + 0.5)) * multiple

facetss = []
## Find facets
for i_face_side in range(n_face_side + 1):
 rad = math.pi * i_face_side / n_face_side
 (r_sin,r_cos) = roundm([r * function(rad) for function in [math.sin,math.cos]])
 if (i_face_side > 0):
  facetss.append([1000 + i_face_side,  r, 0,0, 0,-r_cos_old,-r_sin_old, 0,-r_cos,-r_sin])
  facetss.append([2000 + i_face_side, -r, 0,0, 0, r_cos_old,-r_sin_old, 0, r_cos,-r_sin])
  facetss.append([3000 + i_face_side,  0, r,0, -r_cos_old,0, r_sin_old, -r_cos,0, r_sin])
  facetss.append([4000 + i_face_side,  0,-r,0,  r_cos_old,0, r_sin_old,  r_cos,0, r_sin])
 (r_sin_old,r_cos_old) = (r_sin,r_cos)
## Calculate normals
for facets in facetss:
 us = [facets[i_xyz + 3] - facets[i_xyz] for i_xyz in range(3)]
 vs = [facets[i_xyz] - facets[i_xyz + 6] for i_xyz in range(3)]
 normals = [us[1]*vs[2] - us[2]*vs[1], us[2]*vs[0] - us[0]*vs[2], us[0]*vs[1] - us[1]*vs[0]]
 normal_length = sum([component * component      for component in normals]) ** 0.5
 facets += ['%.5f' % (component / normal_length) for component in normals]
facetss = sorted(facetss)
print(tabbify([['s.f'] + ['%s%d' % (xyz, n) for n in range(3) for xyz in list('XYZ')] +
               ['N%s' % (xyz) for xyz in list('xyz')]] + facetss))
## Compile STL
outs = [fmt('''\
 facet normal {facets[10]} {facets[11]} {facets[12]}
  outer loop
   vertex {facets[1]} {facets[2]} {facets[3]}
   vertex {facets[4]} {facets[5]} {facets[6]}
   vertex {facets[7]} {facets[8]} {facets[9]}
  endloop
 endfacet''') for facets in facetss]
with io.open(__file__[:__file__.rfind('.')] + '.stl', 'w', newline='\n') as f_out:
 f_out.write('solid %s\n%s\nendsolid %s\n\n## Please keep Python script below\n%s' %
             (solid_name, '\n'.join(outs), solid_name, io.open(__file__).read()))