Fix bug in skin, add some complex stuff to math and tests

arrays.scad

@@ -1682,4 +1682,17 @@ function transpose(arr, reverse=false) =
            arr;
 
 
+// Function: is_matrix_symmetric()
+// Usage:
+//   b = is_matrix_symmetric(A,<eps>)
+// Description:
+//   Returns true if the input matrix is symmetric, meaning it equals its transpose.
+//   Matrix should have numerical entries.
+// Arguments:
+//   A = matrix to test
+//   eps = epsilon for comparing equality.  Default: 1e-12
+function is_matrix_symmetric(A,eps=1e-12) =
+    approx(A,transpose(A));
+
+
 // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap

math.scad

@@ -1465,19 +1465,52 @@ function deriv3(data, h=1, closed=false) =
 
 // Section: Complex Numbers
 
+
+// Function: complex()
+// Usage:
+//   z = complex(list)
+// Description:
+//   Converts a real valued number, vector or matrix into its complex analog
+//   by replacing all entries with a 2-vector that has zero imaginary part.
+function complex(list) =
+   is_num(list) ? [list,0] :
+   [for(entry=list) is_num(entry) ? [entry,0] : complex(entry)];
+
+
 // Function: c_mul()
 // Usage:
 //   c = c_mul(z1,z2)
 // Description:
-//   Multiplies two complex numbers represented by 2D vectors.  
-//   Returns a complex number as a 2D vector [REAL, IMAGINARY].
+//   Multiplies two complex numbers, vectors or matrices, where complex numbers
+//   or entries are represented as vectors: [REAL, IMAGINARY].  Note that all
+//   entries in both arguments must be complex.  
 // Arguments:
-//   z1 = First complex number, given as a 2D vector [REAL, IMAGINARY]
-//   z2 = Second complex number, given as a 2D vector [REAL, IMAGINARY]
-function c_mul(z1,z2) = 
-    assert( is_matrix([z1,z2],2,2), "Complex numbers should be represented by 2D vectors" )
+//   z1 = First complex number, vector or matrix
+//   z2 = Second complex number, vector or matrix
+
+function _split_complex(data) =
+   is_vector(data,2) ? data
+ : is_num(data[0][0]) ? [data*[1,0], data*[0,1]]
+ : [
+    [for(vec=data) vec * [1,0]],
+    [for(vec=data) vec * [0,1]]
+   ];
+
+function _combine_complex(data) =
+    is_vector(data,2) ? data
+  : is_num(data[0][0]) ? [for(i=[0:len(data[0])-1]) [data[0][i],data[1][i]]]
+  : [for(i=[0:1:len(data[0])-1])
+        [for(j=[0:1:len(data[0][0])-1])  
+            [data[0][i][j], data[1][i][j]]]];
+
+function _c_mul(z1,z2) = 
     [ z1.x*z2.x - z1.y*z2.y, z1.x*z2.y + z1.y*z2.x ];
 
+function c_mul(z1,z2) =
+  is_matrix([z1,z2],2,2) ? _c_mul(z1,z2) :
+  _combine_complex(_c_mul(_split_complex(z1), _split_complex(z2)));
+
+
 // Function: c_div()
 // Usage:
 //   x = c_div(z1,z2)
@@ -1493,7 +1526,52 @@ function c_div(z1,z2) =
     let(den = z2.x*z2.x + z2.y*z2.y)
     [(z1.x*z2.x + z1.y*z2.y)/den, (z1.y*z2.x - z1.x*z2.y)/den];
 
-// For the sake of consistence with Q_mul and vmul, c_mul should be called C_mul
+
+// Function: c_conj()
+// Usage:
+//   w = c_conj(z)
+// Description:
+//   Computes the complex conjugate of the input, which can be a complex number,
+//   complex vector or complex matrix.  
+function c_conj(z) =
+   is_vector(z,2) ? [z.x,-z.y] :
+   [for(entry=z) c_conj(entry)];
+
+// Function: c_real()
+// Usage:
+//   x = c_real(z)
+// Description:
+//   Returns real part of a complex number, vector or matrix.
+function c_real(z) = 
+     is_vector(z,2) ? z.x
+   : is_num(z[0][0]) ? z*[1,0]
+   : [for(vec=z) vec * [1,0]];
+
+// Function: c_imag()
+// Usage:
+//   x = c_imag(z)
+// Description:
+//   Returns imaginary part of a complex number, vector or matrix.
+function c_imag(z) = 
+     is_vector(z,2) ? z.y
+   : is_num(z[0][0]) ? z*[0,1]
+   : [for(vec=z) vec * [0,1]];
+
+
+// Function: c_ident()
+// Usage:
+//   I = c_ident(n)
+// Description:
+//   Produce an n by n complex identity matrix
+function c_ident(n) = [for (i = [0:1:n-1]) [for (j = [0:1:n-1]) (i==j)?[1,0]:[0,0]]];
+
+// Function: c_norm()
+// Usage:
+//   n = c_norm(z)
+// Description:
+//   Compute the norm of a complex number or vector. 
+function c_norm(z) = norm_fro(z);
+
 
 // Section: Polynomials
 
@@ -1539,12 +1617,12 @@ function quadratic_roots(a,b,c,real=false) =
 //   where a_n is the z^n coefficient.  Polynomial coefficients are real.
 //   The result is a number if `z` is a number and a complex number otherwise.
 function polynomial(p,z,k,total) =
-    is_undef(k)
-    ?   assert( is_vector(p) , "Input polynomial coefficients must be a vector." )
-        assert( is_finite(z) || is_vector(z,2), "The value of `z` must be a real or a complex number." )
-        polynomial( _poly_trim(p), z, 0, is_num(z) ? 0 : [0,0])
-    : k==len(p) ? total
-    : polynomial(p,z,k+1, is_num(z) ? total*z+p[k] : c_mul(total,z)+[p[k],0]);
+  is_undef(k)
+  ? assert( is_vector(p) , "Input polynomial coefficients must be a vector." )
+    assert( is_finite(z) || is_vector(z,2), "The value of `z` must be a real or a complex number." )
+    polynomial( _poly_trim(p), z, 0, is_num(z) ? 0 : [0,0])
+  : k==len(p) ? total
+  : polynomial(p,z,k+1, is_num(z) ? total*z+p[k] : c_mul(total,z)+[p[k],0]);
 
 // Function: poly_mult()
 // Usage:
@@ -1554,12 +1632,12 @@ function polynomial(p,z,k,total) =
 //   Given a list of polynomials represented as real coefficient lists, with the highest degree coefficient first, 
 //   computes the coefficient list of the product polynomial.  
 function poly_mult(p,q) = 
-    is_undef(q) ?
-        len(p)==2 
+  is_undef(q) ?
+    len(p)==2 
         ? poly_mult(p[0],p[1]) 
-        : poly_mult(p[0], poly_mult(select(p,1,-1)))
-    :
-    assert( is_vector(p) && is_vector(q),"Invalid arguments to poly_mult")
+    : poly_mult(p[0], poly_mult(select(p,1,-1)))
+  :
+  assert( is_vector(p) && is_vector(q),"Invalid arguments to poly_mult")
     p*p==0 || q*q==0
     ? [0]
     : _poly_trim(convolve(p,q));

skin.scad

@@ -407,7 +407,7 @@ function skin(profiles, slices, refine=1, method="direct", sampling, caps, close
   assert(methodlistok==[], str("method list contains invalid method at ",methodlistok))
   assert(len(method) == profcount,"Method list is the wrong length")
   assert(in_list(sampling,["length","segment"]), "sampling must be set to \"length\" or \"segment\"")
-  assert(sampling=="segment" || (!in_list("distance",method) && !in_list("fast_distance") && !in_list("tangent",method)), "sampling is set to \"length\" which is only allowed with methods \"direct\" and \"reindex\"")
+  assert(sampling=="segment" || (!in_list("distance",method) && !in_list("fast_distance",method) && !in_list("tangent",method)), "sampling is set to \"length\" which is only allowed with methods \"direct\" and \"reindex\"")
   assert(capsOK, "caps must be boolean or a list of two booleans")
   assert(!closed || !caps, "Cannot make closed shape with caps")
   let(

tests/test_math.scad

@@ -824,19 +824,78 @@ module test_lcm() {
 test_lcm();
 
 
-module test_C_times() {
-    assert_equal(C_times([4,5],[9,-4]), [56,29]);
-    assert_equal(C_times([-7,2],[24,3]), [-174, 27]);
+module test_c_mul() {
+    assert_equal(c_mul([4,5],[9,-4]), [56,29]);
+    assert_equal(c_mul([-7,2],[24,3]), [-174, 27]);
+    assert_equal(c_mul([3,4], [[3,-7], [4,9], [4,8]]), [[37,-9],[-24,43], [-20,40]]);
+    assert_equal(c_mul([[3,-7], [4,9], [4,8]], [[1,1],[3,4],[-3,4]]), [-58,31]);
+    M = [
+           [ [3,4], [9,-1], [4,3] ],
+           [ [2,9], [4,9], [3,-1] ]
+        ];
+    assert_equal(c_mul(M, [ [3,4], [4,4],[5,5]]), [[38,91], [-30, 97]]);
+    assert_equal(c_mul([[4,4],[9,1]], M), [[5,111],[67,117], [32,22]]);
+    assert_equal(c_mul(M,transpose(M)), [  [[80,30], [30, 117]], [[30,117], [-134, 102]]]);
+    assert_equal(c_mul(transpose(M),M), [  [[-84,60],[-42,87],[15,50]], [[-42,87],[15,54],[60,46]], [[15,50],[60,46],[15,18]]]);
 }
-test_C_times();
+test_c_mul();
 
 
-module test_C_div() {
-    assert_equal(C_div([56,29],[9,-4]), [4,5]);
-    assert_equal(C_div([-174,27],[-7,2]), [24,3]);
+module test_c_div() {
+    assert_equal(c_div([56,29],[9,-4]), [4,5]);
+    assert_equal(c_div([-174,27],[-7,2]), [24,3]);
 }    
-test_C_div();
+test_c_div();
 
+module test_c_conj(){
+    assert_equal(c_conj([3,4]), [3,-4]);
+    assert_equal(c_conj(           [ [2,9], [4,9], [3,-1] ]),            [ [2,-9], [4,-9], [3,1] ]);
+    M = [
+           [ [3,4], [9,-1], [4,3] ],
+           [ [2,9], [4,9], [3,-1] ]
+        ];
+    Mc = [
+           [ [3,-4], [9,1], [4,-3] ],
+           [ [2,-9], [4,-9], [3,1] ]
+        ];
+    assert_equal(c_conj(M), Mc);
+}
+test_c_conj();
+
+module test_c_real(){
+    M = [
+           [ [3,4], [9,-1], [4,3] ],
+           [ [2,9], [4,9], [3,-1] ]
+        ];
+    assert_equal(c_real(M), [[3,9,4],[2,4,3]]);
+    assert_equal(c_real(           [ [3,4], [9,-1], [4,3] ]), [3,9,4]);
+    assert_equal(c_real([3,4]),3);
+}
+test_c_real();
+
+
+module test_c_imag(){
+    M = [
+           [ [3,4], [9,-1], [4,3] ],
+           [ [2,9], [4,9], [3,-1] ]
+        ];
+    assert_equal(c_imag(M), [[4,-1,3],[9,9,-1]]);
+    assert_equal(c_imag(           [ [3,4], [9,-1], [4,3] ]), [4,-1,3]);
+    assert_equal(c_imag([3,4]),4);
+}
+test_c_imag();
+
+
+module test_c_ident(){
+  assert_equal(c_ident(3), [[[1, 0], [0, 0], [0, 0]], [[0, 0], [1, 0], [0, 0]], [[0, 0], [0, 0], [1, 0]]]);
+}
+test_c_ident();
+
+module test_c_norm(){
+  assert_equal(c_norm([3,4]), 5);
+  assert_approx(c_norm([[3,4],[5,6]]), 9.273618495495704);
+}
+test_c_norm();
 
 module test_back_substitute(){
    R = [[12,4,3,2],