simplex wip: use gsl_vector_get for X, for abandonment

[moebius3.git] / meshscad

diff --git a/meshscad b/meshscad
index 84cc314e81b8926d302db6fea9a57c82ca599872..ad3194929243001d5942498650ceee7216fe159d 100755 (executable)
--- a/meshscad
+++ b/meshscad
@@ -1,4 +1,4 @@
-#!/usr/bin/python
+#!/usr/bin/python3
 
 from __future__ import print_function
 
@@ -9,54 +9,95 @@ from moebius import *
 from scad import *
 
 nomsize = 30;
-wire = 3.0;
+wire = 1.500;
 
 # number of wires
-kv = 10
-kw = 10
+kv = 24
+kw = 16
 
 # resolution
-nv = 40
-nw = 40
-ns = 10 # around tube
+nv = 200
+nw = 200
+ns = 4 # around tube, should be even
+
+from moedebug import *
+dbg_file(sys.stderr)
 
 nv += -nv % kv
-nw += -hw % kw
+nw += -nw % kw
 
 each_v = nv / kv
 each_w = nw / kw
 
 m = Moebius(nv, nw)
 
-def round_wire(p, vec_surfacenormal, vec_acrosswire, sigma):
-  return (p +
-          wire * (vec_surfacenormal * sin(sigma) +
-                  vec_acrosswire * cos*(sigma)))
+def dpr(v): return '%+3f %+3f %+3f' % tuple(v)
+
+def points_round_wire(p, norm, acrs, sigmas):
+  for sigma in sigmas:
+    delta = norm * sin(sigma) + acrs * cos(sigma)
+    r = p + wire/nomsize * delta
+    yield r
+
+def calc_sigmas(ss):
+  return [ (s + 0.5)/ns * tau for s in ss ]
 
 def make_moebius(objname):
   print('module %s(){' % objname)
-  extents = ScadObject() # wires along extents
-  travers = ScadObject(): # wires along traverses
+  # wires:
+  extents = [ ScadObject() for w in range(0,nw) ] # along extents
+  travers = [ ScadObject() for v in range(0,nv) ] # along traverses
+
+  def qc(v, w, sigmas, is_trav):
+    #print(' QCv,w,T',v,w,is_trav, file=sys.stderr)
+    for ab in 0,1:
+      p, norm, extt, trav = m.details(v + ab*(not is_trav), w + ab*is_trav)
+      if is_trav: acrs = extt
+      else: acrs = trav
+      #print(' RW,ab,sx',ab,sx,
+      #      'r=',dpr(r),
+      #      'p=',dpr(p),
+      #      'norm=',dpr(norm),
+      #      'extt=',dpr(extt),
+      #      'acrs=',dpr(acrs),
+      #      'delta=',dpr(delta),
+      #      's=','%4f' % sigma,
+      #      file=sys.stderr)
+      for r in points_round_wire(p, norm, acrs, sigmas):
+        yield r
+
   for v in range(0, nv):
-    for w in range(0, nw):
+    for w in range(0, nw+1):
+      extw = w
+      if extw > nw/2: extw = nw-w
+      # each extent wire has to go round twice to meet itself
+      # (except the middle one, if there is one, where nw/2 = nw - nw/2
+      # and it gets to be by itself
       for s in range(0, ns):
-        sigmas = [ (s + sx)/ns * tau for sigma in 0,1 ]
-        eq = []
+        sigmas = calc_sigmas([s + sx for sx in (0,1)])
+        #print('VWS',v,w,s, sigmas, file=sys.stderr)
         if not w % each_w:
-          extents[w].quad([ round_wire(p, norm, trav, sigmas[sx])
-                            for p, norm, extt, trav in m.details(v+a, w)
-                            for sx in 0,1
-                            for a in 0,1 ])
-        if not v % each_v:
-          travers[v].quad([ round_wire(p, norm, extt, sigmas[sx])
-                            for p, norm, extt, trav in m.details(v, w+b)
-                            for sx in 0,1
-                            for b in 0,1 ])
-  for v in range(0, nv):
-    travers[v].writeout_core()
+          extents[extw].quad([ cnr for cnr in qc(v,w,sigmas,False) ])
+        if not v % each_v and w < nw:
+          travers[v].rquad([ cnr for cnr in qc(v,w,sigmas,True) ])
+    if not v % each_v:
+      for w in 0, nw:
+        p, norm, extt, trav = m.details(v, w)
+        cnrs = points_round_wire(p, norm, extt, calc_sigmas(range(0,ns)))
+        cnrs = list(cnrs)
+        if w: cnrs.reverse()
+        for s in range(0, ns-1):
+          travers[v].triangle(cnrs[s],
+                              cnrs[s+1],
+                              cnrs[ns-1])
+
   for w in range(0, nw):
-    travers[w].writeout_core()
-  print '}'
+    print('// extent w=', w)
+    extents[w].writeout_core(nomsize)
+  for v in range(0, nv):
+    print('// travers v=', v)
+    travers[v].writeout_core(nomsize)
+  print('}')
 
 make_moebius('MoebiusMesh')
 print('moebiuscore_nomsize=%s;' % repr(nomsize))