curveopt: add a commented-out call to synch everything

[moebius3.git] / symbolic.py

diff --git a/symbolic.py b/symbolic.py
index d505e8c64189614276df37199524829a0149134b..8680378ea06d906cbf58b4d230f09451e5ec281d 100644 (file)
--- a/symbolic.py
+++ b/symbolic.py
@@ -47,15 +47,15 @@ E, H = vector_symbols('E H')
 F0, G0 = vector_symbols('F0 G0')
 En, Hn = vector_symbols('En Hn')
 
-EFl, HGl = symbols('EFl HGl')
+EFlq, HGlq = symbols('EFlq HGlq')
 
 def vector_component(v, ix):
   return v.components[N.base_vectors()[ix]]
 
 # x array in numerical algorithm has:
 #    N x 3    coordinates of points 0..N-3
-#    1        EFl = length parameter |EF| for point 1
-#    1        HGl = length parameter |HG| for point N-2
+#    1        EFlq = sqrt of length parameter |EF| for point 1
+#    1        HGlq = sqrt of length parameter |HG| for point N-2
 
 # fixed array in numerical algorithm has:
 #    4 x 3    E, En, H, Hn
@@ -124,8 +124,8 @@ def calculate():
   # ---------- actual cost computation formulae ----------
 
   global F, G
-  F = E + En * EFl
-  G = H + Hn * HGl
+  F = E + En * pow(EFlq, 2)
+  G = H + Hn * pow(HGlq, 2)
 
   global a,b, al,bl, au,bu
   a,  b  = CoordArray ('a_ b_',   'NP-1', B-A           ).s() # [mm]
@@ -133,15 +133,14 @@ def calculate():
   au, bu = CoordArray ('au bu',   'NP-1', a / al        ).s() # [1]
 
   tan_theta = (au ^ bu) / (au & bu)     # [1]     bending
-  curvature = tan_theta / sqrt(al * bl) # [1/mm]  bending per unit length
 
   global mu, nu
-  mu, nu = CoordArray ('mu nu', 'NP-2', curvature ).s() # [1/mm]
+  mu, nu = CoordArray ('mu nu', 'NP-2', tan_theta ).s() # [1]
 
-  CostComponent('NP-3', sqnorm(mu - nu)) # [1/mm^2]
+  CostComponent('NP-3', sqnorm(mu - nu)) # [1]
 
-  d_density = 1/al - 1/bl # [1/mm]
-  CostComponent('NP-2', pow(d_density, 2)) # [1/mm^2]
+  dl2 = pow(al - bl, 2) # [mm^2]
+  CostComponent('NP-2', dl2 / (al*bl)) # [1]
 
   # ---------- end of cost computation formulae ----------
 
@@ -241,9 +240,9 @@ def gen_point_references():
   cprint(' &X[3*((PP)-2)]')
   cprintraw(')')
 
-  cprintraw('#define EFl X[ NX_DIRECT + 0 ]')
-  cprintraw('#define HGl X[ NX_DIRECT + 1 ]')
-  cprintraw('#define NX   ( NX_DIRECT + 2 )')
+  cprintraw('#define EFlq X[ NX_DIRECT + 0 ]')
+  cprintraw('#define HGlq X[ NX_DIRECT + 1 ]')
+  cprintraw('#define NX    ( NX_DIRECT + 2 )')
 
   for ai in range(0, len(abcd)):
     cprintraw('#define %s POINT(P%+d)' % (abcd[ai], ai))
@@ -261,8 +260,9 @@ def gen_prepare():
                          (H,'H', G0,'G')):
     EFHGv = FG0 - EH
     EFHGl = EFHGv.magnitude()
+    EFHGlq = sqrt(EFHGl)
     cassign_vector(EFHGv/EFHGl, EHs+'n', 'tmp_'+EHs)
-    cassign(EFHGl, EHs+FGs+'l', 'tmp_l'+EHs)
+    cassign(EFHGlq, EHs+FGs+'lq', 'tmp_l'+EHs)
   cprintraw('')
 
 def gen_calculate_FG():