       x=xy[tubeIndex]
       y=xy[tubeIndex+1]
 
-      #vectr from the previous point to this
+      #(dx,dy,dz) is an approximate tangent
       dx = xyNext[tubeIndex] - xyPrev[tubeIndex]
       dy = xyNext[tubeIndex+1] - xyPrev[tubeIndex+1] #y-y0
       #dz = @stepZ
+      iqxyz2 = 1.0 / ( dx*dx+dy*dy+dz*dz)
+      iqxyz = 1.0 / Math.sqrt( iqxyz2 ) #normalizing K
 
       #x0 = x
       #y0 = y
+
+      #main normal is calculated by projecting the last normal 
+      # pn X tangent:
+      tan = (x_pn*dx + y_pn*dy + z_pn*dz) * iqxyz2 #tangent component, scaled by iqxyz
+      xn1 = x_pn - tan * dx
+      yn1 = y_pn - tan * dy
+      zn1 = z_pn - tan * dz
+      #now normalize it
+      inorm1 = 1.0/Math.sqrt(xn1*xn1+yn1*yn1+zn1*zn1)
+      if isFinite inorm1
+        xn1 *= inorm1
+        yn1 *= inorm1
+        zn1 *= inorm1
+      else
+        xn1=1.0
+        yn1=0.0
+        zn1=0.0
+      
 
       #compute normals
       # noraml of form

      # xn1 =   dx / sqrt(dx^2+dy^2)
      # yn1 = - dy / sqrt(dx^2+dy^2)

      qxy = Math.sqrt( dx*dx+dy*dy)
      #zn1 = 0
      if qxy < 1e-3
        xn1 = 1.0
        yn1 = 0.0
      else
        iqxy = 1 / qxy
        xn1 = dy * iqxy
        yn1 = -dx * iqxy
      #now calculate the third vector, as a X-product of
      # (xn1, yn1, 0) X (dx, dy, dz)
      # -dz*yn1*i +dz*xn1*j+ k*(dx*yn1-dy*xn1)
      iqxyz = 1.0 / Math.sqrt( dx*dx+dy*dy+dz*dz)
      xn2 = -dz * iqxyz * yn1
      yn2 =  dz * iqxyz * xn1
      zn2 = (dx*yn1-dy*xn1) * iqxyz

      vindex = curV / 3 | 0

      for i in [0 ... shape.length] by 2
        vx = shape[i]
        vy = shape[i+1]
        pushXYZ x+xn1*vx+xn2*vy, y+yn1*vx+yn2*vy, z+zn2*vy+zn1*vx

      x_pn = xn1
      y_pn = yn1
      z_pn = zn1
            
      if iz > 1 and (Math.abs(xyPrev[tubeIndex] - x)+Math.abs(xyPrev[tubeIndex+1] - y) < jumpTreshold)
        for j in [0...tubeEdges] by 1
          j1 = (j+1)%tubeEdges
          pushQuad vindex-tubeEdges+j, vindex+j, vindex-tubeEdges+j1, vindex+j1

    @prevNormals[tubeIndex/2] = x_pn
    @prevNormals[tubeIndex/2+1] = y_pn
    @prevNormals[tubeIndex/2+2] = z_pn
    if curV isnt vs.length then throw new Error "Not all vertices filled"

    if curIx isnt ixs.length
