diff --git a/affine.scad b/affine.scad index 4a1a7ffc..b39529fb 100644 --- a/affine.scad +++ b/affine.scad @@ -277,9 +277,9 @@ function affine_frame_map(x,y,z, reverse=false) = assert(yvalid,"Input y must be a length 3 vector") assert(zvalid,"Input z must be a length 3 vector") let( - x = is_undef(x)? undef : unit(x), - y = is_undef(y)? undef : unit(y), - z = is_undef(z)? undef : unit(z), + x = is_undef(x)? undef : unit(x,RIGHT), + y = is_undef(y)? undef : unit(y,BACK), + z = is_undef(z)? undef : unit(z,UP), map = is_undef(x)? [cross(y,z), y, z] : is_undef(y)? [x, cross(z,x), z] : is_undef(z)? [x, y, cross(x,y)] : diff --git a/arrays.scad b/arrays.scad index 9d9463a0..ea47d480 100644 --- a/arrays.scad +++ b/arrays.scad @@ -91,7 +91,11 @@ function slice(arr,st,end) = let( // in_list("bar", ["foo", "bar", "baz"]); // Returns true. // in_list("bee", ["foo", "bar", "baz"]); // Returns false. // in_list("bar", [[2,"foo"], [4,"bar"], [3,"baz"]], idx=1); // Returns true. -function in_list(x,l,idx=undef) = search([x], l, num_returns_per_match=1, index_col_num=idx) != [[]]; +function in_list(val,list,idx=undef) = + let( s = search([val], list, num_returns_per_match=1, index_col_num=idx)[0] ) + s==[] ? false + : is_undef(idx) ? val==list[s] + : val==list[s][idx]; // Function: min_index() diff --git a/attachments.scad b/attachments.scad index 323a2ef2..1e7a9fb6 100644 --- a/attachments.scad +++ b/attachments.scad @@ -420,12 +420,12 @@ function find_anchor(anchor, geom) = bot = point3d(vmul(point2d(size)/2,axy),-h/2), top = point3d(vmul(point2d(size2)/2,axy)+shift,h/2), pos = point3d(cp) + lerp(bot,top,u) + offset, - sidevec = unit(rot(from=UP, to=top-bot, p=point3d(axy))), - vvec = anchor==CENTER? UP : unit([0,0,anchor.z]), + sidevec = unit(rot(from=UP, to=top-bot, p=point3d(axy)),UP), + vvec = anchor==CENTER? UP : unit([0,0,anchor.z],UP), vec = anchor==CENTER? UP : - approx(axy,[0,0])? unit(anchor) : + approx(axy,[0,0])? unit(anchor,UP) : approx(anchor.z,0)? sidevec : - unit((sidevec+vvec)/2) + unit((sidevec+vvec)/2,UP) ) [anchor, pos, vec, oang] ) : type == "cyl"? ( //r1, r2, l, shift let( @@ -433,24 +433,24 @@ function find_anchor(anchor, geom) = r1 = is_num(rr1)? [rr1,rr1] : rr1, r2 = is_num(rr2)? [rr2,rr2] : rr2, u = (anchor.z+1)/2, - axy = unit(point2d(anchor)), + axy = unit(point2d(anchor),[0,0]), bot = point3d(vmul(r1,axy), -l/2), top = point3d(vmul(r2,axy)+shift, l/2), pos = point3d(cp) + lerp(bot,top,u) + offset, sidevec = rot(from=UP, to=top-bot, p=point3d(axy)), - vvec = anchor==CENTER? UP : unit([0,0,anchor.z]), + vvec = anchor==CENTER? UP : unit([0,0,anchor.z],UP), vec = anchor==CENTER? UP : - approx(axy,[0,0])? unit(anchor) : + approx(axy,[0,0])? unit(anchor,UP) : approx(anchor.z,0)? sidevec : - unit((sidevec+vvec)/2) + unit((sidevec+vvec)/2,UP) ) [anchor, pos, vec, oang] ) : type == "spheroid"? ( //r let( rr = geom[1], r = is_num(rr)? [rr,rr,rr] : rr, - anchor = unit(point3d(anchor)), + anchor = unit(point3d(anchor),CENTER), pos = point3d(cp) + vmul(r,anchor) + offset, - vec = unit(vmul(r,anchor)) + vec = unit(vmul(r,anchor),UP) ) [anchor, pos, vec, oang] ) : type == "vnf_isect"? ( //vnf let( @@ -490,7 +490,8 @@ function find_anchor(anchor, geom) = faceplane = plane_from_points(faceverts), nrm = plane_normal(faceplane) ) nrm - ]) / len(nfaces) + ]) / len(nfaces), + UP ) ) [anchor, pos, n, oang] @@ -513,15 +514,15 @@ function find_anchor(anchor, geom) = frpt = [size.x/2*anchor.x, -size.y/2], bkpt = [size2/2*anchor.x, size.y/2], pos = point2d(cp) + lerp(frpt, bkpt, u) + offset, - vec = unit(rot(from=BACK, to=bkpt-frpt, p=anchor)) + vec = unit(rot(from=BACK, to=bkpt-frpt, p=anchor),[0,1]) ) [anchor, pos, vec, 0] ) : type == "circle"? ( //r let( rr = geom[1], r = is_num(rr)? [rr,rr] : rr, pos = point2d(cp) + vmul(r,anchor) + offset, - anchor = unit(point2d(anchor)), - vec = unit(vmul([r.y,r.x],anchor)) + anchor = unit(point2d(anchor),[0,0]), + vec = unit(vmul([r.y,r.x],anchor),[0,1]) ) [anchor, pos, vec, 0] ) : type == "path_isect"? ( //path let( @@ -534,7 +535,7 @@ function find_anchor(anchor, geom) = isect = ray_segment_intersection([[0,0],anchor], seg1), n = is_undef(isect)? [0,1] : !approx(isect, t[1])? line_normal(seg1) : - unit((line_normal(seg1)+line_normal(seg2))/2), + unit((line_normal(seg1)+line_normal(seg2))/2,[0,1]), n2 = vector_angle(anchor,n)>90? -n : n ) if(!is_undef(isect) && !approx(isect,t[0])) [norm(isect), isect, n2] @@ -542,7 +543,7 @@ function find_anchor(anchor, geom) = maxidx = max_index(subindex(isects,0)), isect = isects[maxidx], pos = point2d(cp) + isect[1], - vec = unit(isect[2]) + vec = unit(isect[2],[0,1]) ) [anchor, pos, vec, 0] ) : type == "path_extent"? ( //path let( diff --git a/coords.scad b/coords.scad index 97729beb..e45eab12 100644 --- a/coords.scad +++ b/coords.scad @@ -28,7 +28,6 @@ function point2d(p, fill=0) = [for (i=[0:1]) (p[i]==undef)? fill : p[i]]; // every vector has the same length. // Arguments: // points = A list of 2D or 3D points/vectors. -// fill = Value to fill missing values in vectors with. function path2d(points) = assert(is_path(points,dim=undef,fast=true),"Input to path2d is not a path") let (result = points * concat(ident(2), repeat([0,0], len(points[0])-2))) diff --git a/debug.scad b/debug.scad index 3ff5704c..ada55e8b 100644 --- a/debug.scad +++ b/debug.scad @@ -191,16 +191,17 @@ module debug_faces(vertices, faces, size=1, disabled=false) { if (face[0] < 0 || face[1] < 0 || face[2] < 0 || face[0] >= vlen || face[1] >= vlen || face[2] >= vlen) { echo("BAD FACE: ", vlen=vlen, face=face); } else { - v0 = vertices[face[0]]; - v1 = vertices[face[1]]; - v2 = vertices[face[2]]; - c = mean(select(vertices,face)); + verts = select(vertices,face); + c = mean(verts); + v0 = verts[0]; + v1 = verts[1]; + v2 = verts[2]; dv0 = unit(v1 - v0); dv1 = unit(v2 - v0); - nrm0 = unit(cross(dv0, dv1)); - nrm1 = [0, 0, 1]; - axis = unit(cross(nrm0, nrm1)); - ang = vector_angle(nrm0, nrm1); + nrm0 = cross(dv0, dv1); + nrm1 = UP; + axis = vector_axis(nrm0, nrm1); + ang = vector_angle(nrm0, nrm1); theta = atan2(nrm0[1], nrm0[0]); translate(c) { rotate(a=180-ang, v=axis) { diff --git a/geometry.scad b/geometry.scad index b5f4c29c..dd35ddc3 100644 --- a/geometry.scad +++ b/geometry.scad @@ -133,8 +133,9 @@ function distance_from_line(line, pt) = // color("green") stroke([p1,p1+10*n], endcap2="arrow2"); // color("blue") move_copies([p1,p2]) circle(d=2, $fn=12); function line_normal(p1,p2) = - is_undef(p2)? line_normal(p1[0],p1[1]) : - unit([p1.y-p2.y,p2.x-p1.x]); + is_undef(p2)? + assert(is_path(p1,2)) line_normal(p1[0],p1[1]) : + assert(is_vector(p1,2)&&is_vector(p2,2)) unit([p1.y-p2.y,p2.x-p1.x]); // 2D Line intersection from two segments. @@ -252,34 +253,192 @@ function segment_intersection(s1,s2,eps=EPSILON) = // Usage: // line_closest_point(line,pt); // Description: -// Returns the point on the given `line` that is closest to the given point `pt`. +// Returns the point on the given 2D or 3D `line` that is closest to the given point `pt`. +// The `line` and `pt` args should either both be 2D or both 3D. // Arguments: // line = A list of two points that are on the unbounded line. // pt = The point to find the closest point on the line to. +// Example(2D): +// line = [[-30,0],[30,30]]; +// pt = [-32,-10]; +// p2 = line_closest_point(line,pt); +// stroke(line, endcaps="arrow2"); +// color("blue") translate(pt) circle(r=1,$fn=12); +// color("red") translate(p2) circle(r=1,$fn=12); +// Example(2D): +// line = [[-30,0],[30,30]]; +// pt = [-5,0]; +// p2 = line_closest_point(line,pt); +// stroke(line, endcaps="arrow2"); +// color("blue") translate(pt) circle(r=1,$fn=12); +// color("red") translate(p2) circle(r=1,$fn=12); +// Example(2D): +// line = [[-30,0],[30,30]]; +// pt = [40,25]; +// p2 = line_closest_point(line,pt); +// stroke(line, endcaps="arrow2"); +// color("blue") translate(pt) circle(r=1,$fn=12); +// color("red") translate(p2) circle(r=1,$fn=12); +// Example(FlatSpin): +// line = [[-30,-15,0],[30,15,30]]; +// pt = [5,5,5]; +// p2 = line_closest_point(line,pt); +// stroke(line, endcaps="arrow2"); +// color("blue") translate(pt) sphere(r=1,$fn=12); +// color("red") translate(p2) sphere(r=1,$fn=12); +// Example(FlatSpin): +// line = [[-30,-15,0],[30,15,30]]; +// pt = [-35,-15,0]; +// p2 = line_closest_point(line,pt); +// stroke(line, endcaps="arrow2"); +// color("blue") translate(pt) sphere(r=1,$fn=12); +// color("red") translate(p2) sphere(r=1,$fn=12); +// Example(FlatSpin): +// line = [[-30,-15,0],[30,15,30]]; +// pt = [40,15,25]; +// p2 = line_closest_point(line,pt); +// stroke(line, endcaps="arrow2"); +// color("blue") translate(pt) sphere(r=1,$fn=12); +// color("red") translate(p2) sphere(r=1,$fn=12); function line_closest_point(line,pt) = + assert(is_path(line)&&len(line)==2) + assert(same_shape(pt,line[0])) + assert(!approx(line[0],line[1])) let( - n = line_normal(line), - isect = _general_line_intersection(line,[pt,pt+n]) - ) isect[0]; + seglen = norm(line[1]-line[0]), + segvec = (line[1]-line[0])/seglen, + projection = (pt-line[0]) * segvec + ) + line[0] + projection*segvec; + + +// Function: ray_closest_point() +// Usage: +// ray_closest_point(seg,pt); +// Description: +// Returns the point on the given 2D or 3D ray `ray` that is closest to the given point `pt`. +// The `ray` and `pt` args should either both be 2D or both 3D. +// Arguments: +// ray = The ray, given as a list `[START,POINT]` of the start-point START, and a point POINT on the ray. +// pt = The point to find the closest point on the ray to. +// Example(2D): +// ray = [[-30,0],[30,30]]; +// pt = [-32,-10]; +// p2 = ray_closest_point(ray,pt); +// stroke(ray, endcap2="arrow2"); +// color("blue") translate(pt) circle(r=1,$fn=12); +// color("red") translate(p2) circle(r=1,$fn=12); +// Example(2D): +// ray = [[-30,0],[30,30]]; +// pt = [-5,0]; +// p2 = ray_closest_point(ray,pt); +// stroke(ray, endcap2="arrow2"); +// color("blue") translate(pt) circle(r=1,$fn=12); +// color("red") translate(p2) circle(r=1,$fn=12); +// Example(2D): +// ray = [[-30,0],[30,30]]; +// pt = [40,25]; +// p2 = ray_closest_point(ray,pt); +// stroke(ray, endcap2="arrow2"); +// color("blue") translate(pt) circle(r=1,$fn=12); +// color("red") translate(p2) circle(r=1,$fn=12); +// Example(FlatSpin): +// ray = [[-30,-15,0],[30,15,30]]; +// pt = [5,5,5]; +// p2 = ray_closest_point(ray,pt); +// stroke(ray, endcap2="arrow2"); +// color("blue") translate(pt) sphere(r=1,$fn=12); +// color("red") translate(p2) sphere(r=1,$fn=12); +// Example(FlatSpin): +// ray = [[-30,-15,0],[30,15,30]]; +// pt = [-35,-15,0]; +// p2 = ray_closest_point(ray,pt); +// stroke(ray, endcap2="arrow2"); +// color("blue") translate(pt) sphere(r=1,$fn=12); +// color("red") translate(p2) sphere(r=1,$fn=12); +// Example(FlatSpin): +// ray = [[-30,-15,0],[30,15,30]]; +// pt = [40,15,25]; +// p2 = ray_closest_point(ray,pt); +// stroke(ray, endcap2="arrow2"); +// color("blue") translate(pt) sphere(r=1,$fn=12); +// color("red") translate(p2) sphere(r=1,$fn=12); +function ray_closest_point(ray,pt) = + assert(is_path(ray)&&len(ray)==2) + assert(same_shape(pt,ray[0])) + assert(!approx(ray[0],ray[1])) + let( + seglen = norm(ray[1]-ray[0]), + segvec = (ray[1]-ray[0])/seglen, + projection = (pt-ray[0]) * segvec + ) + projection<=0 ? ray[0] : + ray[0] + projection*segvec; // Function: segment_closest_point() // Usage: // segment_closest_point(seg,pt); // Description: -// Returns the point on the given line segment `seg` that is closest to the given point `pt`. +// Returns the point on the given 2D or 3D line segment `seg` that is closest to the given point `pt`. +// The `seg` and `pt` args should either both be 2D or both 3D. // Arguments: // seg = A list of two points that are the endpoints of the bounded line segment. // pt = The point to find the closest point on the segment to. +// Example(2D): +// seg = [[-30,0],[30,30]]; +// pt = [-32,-10]; +// p2 = segment_closest_point(seg,pt); +// stroke(seg); +// color("blue") translate(pt) circle(r=1,$fn=12); +// color("red") translate(p2) circle(r=1,$fn=12); +// Example(2D): +// seg = [[-30,0],[30,30]]; +// pt = [-5,0]; +// p2 = segment_closest_point(seg,pt); +// stroke(seg); +// color("blue") translate(pt) circle(r=1,$fn=12); +// color("red") translate(p2) circle(r=1,$fn=12); +// Example(2D): +// seg = [[-30,0],[30,30]]; +// pt = [40,25]; +// p2 = segment_closest_point(seg,pt); +// stroke(seg); +// color("blue") translate(pt) circle(r=1,$fn=12); +// color("red") translate(p2) circle(r=1,$fn=12); +// Example(FlatSpin): +// seg = [[-30,-15,0],[30,15,30]]; +// pt = [5,5,5]; +// p2 = segment_closest_point(seg,pt); +// stroke(seg); +// color("blue") translate(pt) sphere(r=1,$fn=12); +// color("red") translate(p2) sphere(r=1,$fn=12); +// Example(FlatSpin): +// seg = [[-30,-15,0],[30,15,30]]; +// pt = [-35,-15,0]; +// p2 = segment_closest_point(seg,pt); +// stroke(seg); +// color("blue") translate(pt) sphere(r=1,$fn=12); +// color("red") translate(p2) sphere(r=1,$fn=12); +// Example(FlatSpin): +// seg = [[-30,-15,0],[30,15,30]]; +// pt = [40,15,25]; +// p2 = segment_closest_point(seg,pt); +// stroke(seg); +// color("blue") translate(pt) sphere(r=1,$fn=12); +// color("red") translate(p2) sphere(r=1,$fn=12); function segment_closest_point(seg,pt) = + assert(is_path(seg)&&len(seg)==2) + assert(same_shape(pt,seg[0])) + approx(seg[0],seg[1])? seg[0] : let( - n = line_normal(seg), - isect = _general_line_intersection(seg,[pt,pt+n]) + seglen = norm(seg[1]-seg[0]), + segvec = (seg[1]-seg[0])/seglen, + projection = (pt-seg[0]) * segvec ) - norm(n)==0? seg[0] : - isect[1]<=0? seg[0] : - isect[1]>=1? seg[1] : - isect[0]; + projection<=0 ? seg[0] : + projection>=seglen ? seg[1] : + seg[0] + projection*segvec; // Section: 2D Triangles @@ -1221,15 +1380,18 @@ function find_noncollinear_points(points) = // Usage: // pointlist_bounds(pts); // Description: -// Finds the bounds containing all the 2D or 3D points in `pts`. -// Returns `[[MINX, MINY, MINZ], [MAXX, MAXY, MAXZ]]` +// Finds the bounds containing all the points in `pts` which can be a list of points in any dimension. +// Returns a list of two items: a list of the minimums and a list of the maximums. For example, with +// 3d points `[[MINX, MINY, MINZ], [MAXX, MAXY, MAXZ]]` // Arguments: // pts = List of points. -function pointlist_bounds(pts) = [ - [for (a=[0:2]) min([ for (x=pts) point3d(x)[a] ]) ], - [for (a=[0:2]) max([ for (x=pts) point3d(x)[a] ]) ] -]; - +function pointlist_bounds(pts) = + assert(is_matrix(pts)) + let(ptsT = transpose(pts)) + [ + [for(row=ptsT) min(row)], + [for(row=ptsT) max(row)] + ]; // Function: closest_point() // Usage: diff --git a/math.scad b/math.scad index 09e3b9ce..7629cb18 100644 --- a/math.scad +++ b/math.scad @@ -67,7 +67,7 @@ function hypot(x,y,z=0) = norm([x,y,z]); // Usage: // x = factorial(n,[d]); // Description: -// Returns the factorial of the given integer value. +// Returns the factorial of the given integer value, or n!/d! if d is given. // Arguments: // n = The integer number to get the factorial of. (n!) // d = If given, the returned value will be (n! / d!) @@ -75,7 +75,10 @@ function hypot(x,y,z=0) = norm([x,y,z]); // x = factorial(4); // Returns: 24 // y = factorial(6); // Returns: 720 // z = factorial(9); // Returns: 362880 -function factorial(n,d=1) = product([for (i=[n:-1:d]) i]); +function factorial(n,d=0) = + assert(n>=0 && d>=0, "Factorial is not defined for negative numbers") + assert(d<=n, "d cannot be larger than n") + product([1,for (i=[n:-1:d+1]) i]); // Function: lerp() @@ -525,6 +528,17 @@ function deltas(v) = [for (p=pair(v)) p.y-p.x]; function product(v, i=0, tot=undef) = i>=len(v)? tot : product(v, i+1, ((tot==undef)? v[i] : is_vector(v[i])? vmul(tot,v[i]) : tot*v[i])); +// Function: outer_product() +// Description: +// Compute the outer product of two vectors, a matrix. +// Usage: +// M = outer_product(u,v); +function outer_product(u,v) = + assert(is_vector(u) && is_vector(v)) + assert(len(u)==len(v)) + [for(i=[0:len(u)-1]) [for(j=[0:len(u)-1]) u[i]*v[j]]]; + + // Function: mean() // Description: // Returns the arithmatic mean/average of all entries in the given array. @@ -563,7 +577,7 @@ function median(v) = // Description: // Solves the linear system Ax=b. If A is square and non-singular the unique solution is returned. If A is overdetermined // the least squares solution is returned. If A is underdetermined, the minimal norm solution is returned. -// If A is rank deficient or singular then linear_solve returns `undef`. If b is a matrix that is compatible with A +// If A is rank deficient or singular then linear_solve returns []. If b is a matrix that is compatible with A // then the problem is solved for the matrix valued right hand side and a matrix is returned. Note that if you // want to solve Ax=b1 and Ax=b2 that you need to form the matrix transpose([b1,b2]) for the right hand side and then // transpose the returned value. @@ -582,7 +596,7 @@ function linear_solve(A,b) = R = submatrix(qr[1],[0:mindim-1], [0:mindim-1]), zeros = [for(i=[0:mindim-1]) if (approx(R[i][i],0)) i] ) - zeros != [] ? undef : + zeros != [] ? [] : m=0 && j0 && d[0]!=0 , "Denominator is zero or has leading zero coefficient") + len(n)qlen ? p : q, + short = plen>qlen ? q : p + ) + _poly_trim(long + concat(repeat(0,len(long)-len(short)),short)); + + // Function: poly_roots() // Usage: // poly_roots(p,[tol]) @@ -1062,14 +1151,18 @@ function polynomial(p, z, k, zk, total) = // Dario Bini. "Numerical computation of polynomial zeros by means of Aberth's Method", Numerical Algorithms, Feb 1996. // https://www.researchgate.net/publication/225654837_Numerical_computation_of_polynomial_zeros_by_means_of_Aberth's_method -function poly_roots(p,tol=1e-14) = +function poly_roots(p,tol=1e-14,error_bound=false) = assert(p!=[], "Input polynomial must have a nonzero coefficient") assert(is_vector(p), "Input must be a vector") - p[0] == 0 ? poly_roots(slice(p,1,-1)) : // Strip leading zero coefficients - p[len(p)-1] == 0 ? [[0,0], // Strip trailing zero coefficients - each poly_roots(select(p,0,-2))] : - len(p)==1 ? [] : // Nonzero constant case has no solutions - len(p)==2 ? [[-p[1]/p[0],0]] : // Linear case needs special handling + p[0] == 0 ? poly_roots(slice(p,1,-1),tol=tol,error_bound=error_bound) : // Strip leading zero coefficients + p[len(p)-1] == 0 ? // Strip trailing zero coefficients + let( solutions = poly_roots(select(p,0,-2),tol=tol, error_bound=error_bound)) + (error_bound ? [ [[0,0], each solutions[0]], [0, each solutions[1]]] + : [[0,0], each solutions]) : + len(p)==1 ? (error_bound ? [[],[]] : []) : // Nonzero constant case has no solutions + len(p)==2 ? let( solution = [[-p[1]/p[0],0]]) // Linear case needs special handling + (error_bound ? [solution,[0]] : solution) + : let( n = len(p)-1, // polynomial degree pderiv = [for(i=[0:n-1]) p[i]*(n-i)], @@ -1082,9 +1175,12 @@ function poly_roots(p,tol=1e-14) = init = [for(i=[0:1:n-1]) // Initial guess for roots let(angle = 360*i/n+270/n/PI) [beta,0]+r*[cos(angle),sin(angle)] - ] + ], + roots = _poly_roots(p,pderiv,s,init,tol=tol), + error = error_bound ? [for(xi=roots) n * (norm(polynomial(p,xi))+tol*polynomial(s,norm(xi))) / + abs(norm(polynomial(pderiv,xi))-tol*polynomial(s,norm(xi)))] : 0 ) - _poly_roots(p,pderiv,s,init,tol=tol); + error_bound ? [roots, error] : roots; // p = polynomial // pderiv = derivative polynomial of p @@ -1114,7 +1210,10 @@ function _poly_roots(p, pderiv, s, z, tol, i=0) = // The polynomial is specified as p=[a_n, a_{n-1},...,a_1,a_0] // where a_n is the x^n coefficient. This function works by // computing the complex roots and returning those roots where -// the imaginary part is closed to zero, specifically: z.y/(1+norm(z)) < eps. Because +// the imaginary part is closed to zero. By default it uses a computed +// error bound from the polynomial solver to decide whether imaginary +// parts are zero. You can specify eps, in which case the test is +// z.y/(1+norm(z)) < eps. Because // of poor convergence and higher error for repeated roots, such roots may // be missed by the algorithm because their imaginary part is large. // Arguments: @@ -1122,11 +1221,13 @@ function _poly_roots(p, pderiv, s, z, tol, i=0) = // eps = used to determine whether imaginary parts of roots are zero // tol = tolerance for the complex polynomial root finder -function real_roots(p,eps=EPSILON,tol=1e-14) = +function real_roots(p,eps=undef,tol=1e-14) = let( - roots = poly_roots(p) + roots_err = poly_roots(p,error_bound=true), + roots = roots_err[0], + err = roots_err[1] ) - [for(z=roots) if (abs(z.y)/(1+norm(z))0 ? Q_Ident() : Quat([ v1.y, -v1.x, 0], 180) + : Quat(ax, atan2( n , v1*v2 )); + + // Function: Q_Ident() // Description: Returns the "Identity" zero-rotation Quaternion. function Q_Ident() = [0, 0, 0, 1]; @@ -81,55 +129,85 @@ function Q_Ident() = [0, 0, 0, 1]; // Function: Q_Add_S() // Usage: // Q_Add_S(q, s) -// Description: Adds a scalar value `s` to the W part of a quaternion `q`. -function Q_Add_S(q, s) = q+[0,0,0,s]; +// Description: +// Adds a scalar value `s` to the W part of a quaternion `q`. +// The returned quaternion is usually not normalized. +function Q_Add_S(q, s) = + assert( is_finite(s), "Invalid scalar" ) + q+[0,0,0,s]; // Function: Q_Sub_S() // Usage: // Q_Sub_S(q, s) -// Description: Subtracts a scalar value `s` from the W part of a quaternion `q`. -function Q_Sub_S(q, s) = q-[0,0,0,s]; +// Description: +// Subtracts a scalar value `s` from the W part of a quaternion `q`. +// The returned quaternion is usually not normalized. +function Q_Sub_S(q, s) = + assert( is_finite(s), "Invalid scalar" ) + q-[0,0,0,s]; // Function: Q_Mul_S() // Usage: // Q_Mul_S(q, s) -// Description: Multiplies each part of a quaternion `q` by a scalar value `s`. -function Q_Mul_S(q, s) = q*s; +// Description: +// Multiplies each part of a quaternion `q` by a scalar value `s`. +// The returned quaternion is usually not normalized. +function Q_Mul_S(q, s) = + assert( is_finite(s), "Invalid scalar" ) + q*s; // Function: Q_Div_S() // Usage: // Q_Div_S(q, s) -// Description: Divides each part of a quaternion `q` by a scalar value `s`. -function Q_Div_S(q, s) = q/s; +// Description: +// Divides each part of a quaternion `q` by a scalar value `s`. +// The returned quaternion is usually not normalized. +function Q_Div_S(q, s) = + assert( is_finite(s) && ! approx(s,0) , "Invalid scalar" ) + q/s; // Function: Q_Add() // Usage: // Q_Add(a, b) -// Description: Adds each part of two quaternions together. -function Q_Add(a, b) = a+b; +// Description: +// Adds each part of two quaternions together. +// The returned quaternion is usually not normalized. +function Q_Add(a, b) = + assert( Q_is_quat(a) && Q_is_quat(a), "Invalid quaternion(s)") + assert( ! approx(norm(a+b),0), "Quaternions cannot be opposed" ) + a+b; // Function: Q_Sub() // Usage: // Q_Sub(a, b) -// Description: Subtracts each part of quaternion `b` from quaternion `a`. -function Q_Sub(a, b) = a-b; +// Description: +// Subtracts each part of quaternion `b` from quaternion `a`. +// The returned quaternion is usually not normalized. +function Q_Sub(a, b) = + assert( Q_is_quat(a) && Q_is_quat(a), "Invalid quaternion(s)") + assert( ! approx(a,b), "Quaternions cannot be equal" ) + a-b; // Function: Q_Mul() // Usage: // Q_Mul(a, b) -// Description: Multiplies quaternion `a` by quaternion `b`. -function Q_Mul(a, b) = [ - a[3]*b.x + a.x*b[3] + a.y*b.z - a.z*b.y, - a[3]*b.y - a.x*b.z + a.y*b[3] + a.z*b.x, - a[3]*b.z + a.x*b.y - a.y*b.x + a.z*b[3], - a[3]*b[3] - a.x*b.x - a.y*b.y - a.z*b.z, -]; +// Description: +// Multiplies quaternion `a` by quaternion `b`. +// The returned quaternion is normalized if both `a` and `b` are normalized +function Q_Mul(a, b) = + assert( Q_is_quat(a) && Q_is_quat(b), "Invalid quaternion(s)") + [ + a[3]*b.x + a.x*b[3] + a.y*b.z - a.z*b.y, + a[3]*b.y - a.x*b.z + a.y*b[3] + a.z*b.x, + a[3]*b.z + a.x*b.y - a.y*b.x + a.z*b[3], + a[3]*b[3] - a.x*b.x - a.y*b.y - a.z*b.z, + ]; // Function: Q_Cumulative() @@ -139,6 +217,8 @@ function Q_Mul(a, b) = [ // Given a list of Quaternions, cumulatively multiplies them, returning a list // of each cumulative Quaternion product. It starts with the first quaternion // given in the list, and applies successive quaternion rotations in list order. +// The quaternion in the returned list are normalized if each quaternion in v +// is normalized. function Q_Cumulative(v, _i=0, _acc=[]) = _i==len(v) ? _acc : Q_Cumulative( @@ -154,42 +234,65 @@ function Q_Cumulative(v, _i=0, _acc=[]) = // Usage: // Q_Dot(a, b) // Description: Calculates the dot product between quaternions `a` and `b`. -function Q_Dot(a, b) = a[0]*b[0] + a[1]*b[1] + a[2]*b[2] + a[3]*b[3]; - +function Q_Dot(a, b) = + assert( Q_is_quat(a) && Q_is_quat(b), "Invalid quaternion(s)" ) + a*b; // Function: Q_Neg() // Usage: // Q_Neg(q) // Description: Returns the negative of quaternion `q`. -function Q_Neg(q) = -q; +function Q_Neg(q) = + assert( Q_is_quat(q), "Invalid quaternion" ) + -q; // Function: Q_Conj() // Usage: // Q_Conj(q) // Description: Returns the conjugate of quaternion `q`. -function Q_Conj(q) = [-q.x, -q.y, -q.z, q[3]]; +function Q_Conj(q) = + assert( Q_is_quat(q), "Invalid quaternion" ) + [-q.x, -q.y, -q.z, q[3]]; + + +// Function: Q_Inverse() +// Usage: +// qc = Q_Inverse(q) +// Description: Returns the multiplication inverse of quaternion `q` that is normalized only if `q` is normalized. +function Q_Inverse(q) = + assert( Q_is_quat(q), "Invalid quaternion" ) + let(q = _Qnorm(q) ) + [-q.x, -q.y, -q.z, q[3]]; // Function: Q_Norm() // Usage: // Q_Norm(q) -// Description: Returns the `norm()` "length" of quaternion `q`. -function Q_Norm(q) = norm(q); +// Description: +// Returns the `norm()` "length" of quaternion `q`. +// Normalized quaternions have unitary norm. +function Q_Norm(q) = + assert( Q_is_quat(q), "Invalid quaternion" ) + norm(q); // Function: Q_Normalize() // Usage: // Q_Normalize(q) // Description: Normalizes quaternion `q`, so that norm([W,X,Y,Z]) == 1. -function Q_Normalize(q) = q/norm(q); +function Q_Normalize(q) = + assert( Q_is_quat(q) , "Invalid quaternion" ) + q/norm(q); // Function: Q_Dist() // Usage: // Q_Dist(q1, q2) // Description: Returns the "distance" between two quaternions. -function Q_Dist(q1, q2) = norm(q2-q1); +function Q_Dist(q1, q2) = + assert( Q_is_quat(q1) && Q_is_quat(q2), "Invalid quaternion(s)" ) + norm(q2-q1); // Function: Q_Slerp() @@ -214,24 +317,24 @@ function Q_Dist(q1, q2) = norm(q2-q1); // for (q = Q_Slerp(a, b, [0:0.1:1])) // Qrot(q) right(80) cube([10,10,1]); // #sphere(r=80); -function Q_Slerp(q1, q2, u) = - assert(is_num(u) || is_num(u[0])) - !is_num(u)? [for (uu=u) Q_Slerp(q1,q2,uu)] : - let( - q1 = Q_Normalize(q1), - q2 = Q_Normalize(q2), - dot = Q_Dot(q1, q2) - ) let( - q2 = dot<0? Q_Neg(q2) : q2, - dot = dot<0? -dot : dot - ) (dot>0.9995)? Q_Normalize(q1 + (u * (q2-q1))) : - let( - dot = constrain(dot,-1,1), - theta_0 = acos(dot), - theta = theta_0 * u, - q3 = Q_Normalize(q2 - q1*dot), - out = q1*cos(theta) + q3*sin(theta) - ) out; +function Q_Slerp(q1, q2, u, _dot) = + is_undef(_dot) + ? assert(is_finite(u) || is_range(u) || is_vector(u), "Invalid interpolation coefficient(s)") + assert(Q_is_quat(q1) && Q_is_quat(q2), "Invalid quaternion(s)" ) + let( + _dot = q1*q2, + q1 = q1/norm(q1), + q2 = _dot<0 ? -q2/norm(q2) : q2/norm(q2), + dot = abs(_dot) + ) + ! is_finite(u) ? [for (uu=u) Q_Slerp(q1, q2, uu, dot)] : + Q_Slerp(q1, q2, u, dot) + : _dot>0.9995 + ? _Qnorm(q1 + u*(q2-q1)) + : let( theta = u*acos(_dot), + q3 = _Qnorm(q2 - _dot*q1) + ) + _Qnorm(q1*cos(theta) + q3*sin(theta)); // Function: Q_Matrix3() @@ -239,11 +342,13 @@ function Q_Slerp(q1, q2, u) = // Q_Matrix3(q); // Description: // Returns the 3x3 rotation matrix for the given normalized quaternion q. -function Q_Matrix3(q) = [ - [1-2*q[1]*q[1]-2*q[2]*q[2], 2*q[0]*q[1]-2*q[2]*q[3], 2*q[0]*q[2]+2*q[1]*q[3]], - [ 2*q[0]*q[1]+2*q[2]*q[3], 1-2*q[0]*q[0]-2*q[2]*q[2], 2*q[1]*q[2]-2*q[0]*q[3]], - [ 2*q[0]*q[2]-2*q[1]*q[3], 2*q[1]*q[2]+2*q[0]*q[3], 1-2*q[0]*q[0]-2*q[1]*q[1]] -]; +function Q_Matrix3(q) = + let( q = Q_Normalize(q) ) + [ + [1-2*q[1]*q[1]-2*q[2]*q[2], 2*q[0]*q[1]-2*q[2]*q[3], 2*q[0]*q[2]+2*q[1]*q[3]], + [ 2*q[0]*q[1]+2*q[2]*q[3], 1-2*q[0]*q[0]-2*q[2]*q[2], 2*q[1]*q[2]-2*q[0]*q[3]], + [ 2*q[0]*q[2]-2*q[1]*q[3], 2*q[1]*q[2]+2*q[0]*q[3], 1-2*q[0]*q[0]-2*q[1]*q[1]] + ]; // Function: Q_Matrix4() @@ -251,12 +356,14 @@ function Q_Matrix3(q) = [ // Q_Matrix4(q); // Description: // Returns the 4x4 rotation matrix for the given normalized quaternion q. -function Q_Matrix4(q) = [ - [1-2*q[1]*q[1]-2*q[2]*q[2], 2*q[0]*q[1]-2*q[2]*q[3], 2*q[0]*q[2]+2*q[1]*q[3], 0], - [ 2*q[0]*q[1]+2*q[2]*q[3], 1-2*q[0]*q[0]-2*q[2]*q[2], 2*q[1]*q[2]-2*q[0]*q[3], 0], - [ 2*q[0]*q[2]-2*q[1]*q[3], 2*q[1]*q[2]+2*q[0]*q[3], 1-2*q[0]*q[0]-2*q[1]*q[1], 0], - [ 0, 0, 0, 1] -]; +function Q_Matrix4(q) = + let( q = Q_Normalize(q) ) + [ + [1-2*q[1]*q[1]-2*q[2]*q[2], 2*q[0]*q[1]-2*q[2]*q[3], 2*q[0]*q[2]+2*q[1]*q[3], 0], + [ 2*q[0]*q[1]+2*q[2]*q[3], 1-2*q[0]*q[0]-2*q[2]*q[2], 2*q[1]*q[2]-2*q[0]*q[3], 0], + [ 2*q[0]*q[2]-2*q[1]*q[3], 2*q[1]*q[2]+2*q[0]*q[3], 1-2*q[0]*q[0]-2*q[1]*q[1], 0], + [ 0, 0, 0, 1] + ]; // Function: Q_Axis() @@ -264,16 +371,28 @@ function Q_Matrix4(q) = [ // Q_Axis(q) // Description: // Returns the axis of rotation of a normalized quaternion `q`. -function Q_Axis(q) = let(d = sqrt(1-(q[3]*q[3]))) (d==0)? [0,0,1] : [q[0]/d, q[1]/d, q[2]/d]; - +// The input doesn't need to be normalized. +function Q_Axis(q) = + assert( Q_is_quat(q) , "Invalid quaternion" ) + let( d = norm(_Qvec(q)) ) + approx(d,0)? [0,0,1] : _Qvec(q)/d; // Function: Q_Angle() // Usage: -// Q_Angle(q) +// a = Q_Angle(q) +// a12 = Q_Angle(q1,q2); // Description: -// Returns the angle of rotation (in degrees) of a normalized quaternion `q`. -function Q_Angle(q) = 2 * acos(q[3]); - +// If only q1 is given, returns the angle of rotation (in degrees) of that quaternion. +// If both q1 and q2 are given, returns the angle (in degrees) between them. +// The input quaternions don't need to be normalized. +function Q_Angle(q1,q2) = + assert(Q_is_quat(q1) && (is_undef(q2) || Q_is_quat(q2)), "Invalid quaternion(s)" ) + let( n1 = is_undef(q2)? norm(_Qvec(q1)): norm(q1) ) + is_undef(q2) + ? 2 * atan2(n1,_Qreal(q1)) + : let( q1 = q1/norm(q1), + q2 = q2/norm(q2) ) + 4 * atan2(norm(q1 - q2), norm(q1 + q2)); // Function&Module: Qrot() // Usage: As Module @@ -305,9 +424,9 @@ module Qrot(q) { } function Qrot(q,p) = - is_undef(p)? Q_Matrix4(q) : - is_vector(p)? Qrot(q,[p])[0] : - apply(Q_Matrix4(q), p); + is_undef(p)? Q_Matrix4(q) : + is_vector(p)? Qrot(q,[p])[0] : + apply(Q_Matrix4(q), p); // Module: Qrot_copies() @@ -327,4 +446,214 @@ function Qrot(q,p) = module Qrot_copies(quats) for (q=quats) Qrot(q) children(); +// Function: Q_Rotation() +// Usage: +// Q_Rotation(R) +// Description: +// Returns a normalized quaternion corresponding to the rotation matrix R. +// R may be a 3x3 rotation matrix or a homogeneous 4x4 rotation matrix. +// The last row and last column of R are ignored for 4x4 matrices. +// It doesn't check whether R is in fact a rotation matrix. +// If R is not a rotation, the returned quaternion is an unpredictable quaternion . +function Q_Rotation(R) = + assert( is_matrix(R,3,3) || is_matrix(R,4,4) , + "Matrix is neither 3x3 nor 4x4") + let( tr = R[0][0]+R[1][1]+R[2][2] ) // R trace + tr>0 + ? let( r = 1+tr ) + _Qnorm( _Qset([ R[1][2]-R[2][1], R[2][0]-R[0][2], R[0][1]-R[1][0] ], -r ) ) + : let( i = max_index([ R[0][0], R[1][1], R[2][2] ]), + r = 1 + 2*R[i][i] -R[0][0] -R[1][1] -R[2][2] ) + i==0 ? _Qnorm( _Qset( [ 4*r, (R[1][0]+R[0][1]), (R[0][2]+R[2][0]) ], (R[2][1]-R[1][2])) ): + i==1 ? _Qnorm( _Qset( [ (R[1][0]+R[0][1]), 4*r, (R[2][1]+R[1][2]) ], (R[0][2]-R[2][0])) ): + _Qnorm( _Qset( [ (R[2][0]+R[0][2]), (R[1][2]+R[2][1]), 4*r ], (R[1][0]-R[0][1])) ) ; + + +// Function&Module: Q_Rotation_path(q1, n, [q2]) +// Usage: As a function +// path = Q_Rotation_path(q1, n, q2); +// path = Q_Rotation_path(q1, n); +// Usage: As a module +// Q_Rotation_path(q1, n, q2) ... +// Description: +// If q2 is undef and it is called as a function, the path, with length n+1 (n>=1), will be the +// cumulative multiplications of the matrix rotation of q1 by itself. +// If q2 is defined and it is called as a function, returns a rotation matrix path of length n+1 (n>=1) +// that interpolates two given rotation quaternions. The first matrix of the sequence is the +// matrix rotation of q1 and the last one, the matrix rotation of q2. The intermediary matrix +// rotations are an uniform interpolation of the path extreme matrices. +// When called as a module, applies to its children() each rotation of the sequence computed +// by the function. +// The input quaternions don't need to be normalized. +// Arguments: +// q1 = The quaternion of the first rotation. +// q2 = The quaternion of the last rotation. +// n = An integer defining the path length ( path length = n+1). +// Example(3D): as a function +// a = QuatY(-135); +// b = QuatXYZ([0,-30,30]); +// for (M=Q_Rotation_path(a, 10, b)) +// multmatrix(M) +// right(80) cube([10,10,1]); +// #sphere(r=80); +// Example(3D): as a module +// a = QuatY(-135); +// b = QuatXYZ([0,-30,30]); +// Q_Rotation_path(a, 10, b) +// right(80) cube([10,10,1]); +// #sphere(r=80); +// Example(3D): as a function +// a = QuatY(5); +// for (M=Q_Rotation_path(a, 10)) +// multmatrix(M) +// right(80) cube([10,10,1]); +// #sphere(r=80); +// Example(3D): as a module +// a = QuatY(5); +// Q_Rotation_path(a, 10) +// right(80) cube([10,10,1]); +// #sphere(r=80); +function Q_Rotation_path(q1, n=1, q2) = + assert( Q_is_quat(q1) && (is_undef(q2) || Q_is_quat(q2) ), "Invalid quaternion(s)" ) + assert( is_finite(n) && n>=1 && n==floor(n), "Invalid integer" ) + assert( is_undef(q2) || ! approx(norm(q1+q2),0), "Quaternions cannot be opposed" ) + is_undef(q2) + ? [for( i=0, dR=Q_Matrix4(q1), R=dR; i<=n; i=i+1, R=dR*R ) R] + : let( q2 = Q_Normalize( q1*q2<0 ? -q2: q2 ), + dq = Q_pow( Q_Mul( q2, Q_Inverse(q1) ), 1/n ), + dR = Q_Matrix4(dq) ) + [for( i=0, R=Q_Matrix4(q1); i<=n; i=i+1, R=dR*R ) R]; + +module Q_Rotation_path(q1, n=1, q2) { + for(Mi=Q_Rotation_path(q1, n, q2)) + multmatrix(Mi) + children(); +} + + +// Function: Q_Nlerp() +// Usage: +// q = Q_Nlerp(q1, q2, u); +// Description: +// Returns a quaternion that is a normalized linear interpolation between two quaternions +// when u is a number. +// If u is a list of numbers, computes the interpolations for each value in the +// list and returns the interpolated quaternions in a list. +// The input quaternions don't need to be normalized. +// Arguments: +// q1 = The first quaternion. (u=0) +// q2 = The second quaternion. (u=1) +// u = A value (or a list of values), between 0 and 1, of the proportion(s) of each quaternion in the interpolation. +// Example(3D): Giving `u` as a Scalar +// a = QuatY(-135); +// b = QuatXYZ([0,-30,30]); +// for (u=[0:0.1:1]) +// Qrot(Q_Nlerp(a, b, u)) +// right(80) cube([10,10,1]); +// #sphere(r=80); +// Example(3D): Giving `u` as a Range +// a = QuatZ(-135); +// b = QuatXYZ([90,0,-45]); +// for (q = Q_Nlerp(a, b, [0:0.1:1])) +// Qrot(q) right(80) cube([10,10,1]); +// #sphere(r=80); +function Q_Nlerp(q1,q2,u) = + assert(is_finite(u) || is_range(u) || is_vector(u) , + "Invalid interpolation coefficient(s)" ) + assert(Q_is_quat(q1) && Q_is_quat(q2), "Invalid quaternion(s)" ) + assert( ! approx(norm(q1+q2),0), "Quaternions cannot be opposed" ) + let( q1 = Q_Normalize(q1), + q2 = Q_Normalize(q2) ) + is_num(u) + ? _Qnorm((1-u)*q1 + u*q2 ) + : [for (ui=u) _Qnorm((1-ui)*q1 + ui*q2 ) ]; + + +// Function: Q_Squad() +// Usage: +// qn = Q_Squad(q1,q2,q3,q4,u); +// Description: +// Returns a quaternion that is a cubic spherical interpolation of the quaternions +// q1 and q4 taking the other two quaternions, q2 and q3, as parameter of a cubic +// on the sphere similar to the control points of a Bezier curve. +// If u is a number, usually between 0 and 1, returns the quaternion that results +// from the interpolation. +// If u is a list of numbers, computes the interpolations for each value in the +// list and returns the interpolated quaternions in a list. +// The input quaternions don't need to be normalized. +// Arguments: +// q1 = The start quaternion. (u=0) +// q1 = The first intermediate quaternion. +// q2 = The second intermediate quaternion. +// q4 = The end quaternion. (u=1) +// u = A value (or a list of values), of the proportion(s) of each quaternion in the cubic interpolation. +// Example(3D): Giving `u` as a Scalar +// a = QuatY(-135); +// b = QuatXYZ([-50,-50,120]); +// c = QuatXYZ([-50,-40,30]); +// d = QuatY(-45); +// color("red"){ +// Qrot(b) right(80) cube([10,10,1]); +// Qrot(c) right(80) cube([10,10,1]); +// } +// for (u=[0:0.05:1]) +// Qrot(Q_Squad(a, b, c, d, u)) +// right(80) cube([10,10,1]); +// #sphere(r=80); +// Example(3D): Giving `u` as a Range +// a = QuatY(-135); +// b = QuatXYZ([-50,-50,120]); +// c = QuatXYZ([-50,-40,30]); +// d = QuatY(-45); +// for (q = Q_Squad(a, b, c, d, [0:0.05:1])) +// Qrot(q) right(80) cube([10,10,1]); +// #sphere(r=80); +function Q_Squad(q1,q2,q3,q4,u) = + assert(is_finite(u) || is_range(u) || is_vector(u) , + "Invalid interpolation coefficient(s)" ) + is_num(u) + ? Q_Slerp( Q_Slerp(q1,q4,u), Q_Slerp(q2,q3,u), 2*u*(1-u)) + : [for(ui=u) Q_Slerp( Q_Slerp(q1,q4,ui), Q_Slerp(q2,q3,ui), 2*ui*(1-ui) ) ]; + + +// Function: Q_exp() +// Usage: +// q2 = Q_exp(q); +// Description: +// Returns the quaternion that is the exponential of the quaternion q in base e +// The returned quaternion is usually not normalized. +function Q_exp(q) = + assert( is_vector(q,4), "Input is not a valid quaternion") + let( nv = norm(_Qvec(q)) ) // q may be equal to zero here! + exp(_Qreal(q))*Quat(_Qvec(q),2*nv); + + +// Function: Q_ln() +// Usage: +// q2 = Q_ln(q); +// Description: +// Returns the quaternion that is the natural logarithm of the quaternion q. +// The returned quaternion is usually not normalized and may be zero. +function Q_ln(q) = + assert(Q_is_quat(q), "Input is not a valid quaternion") + let( nq = norm(q), + nv = norm(_Qvec(q)) ) + approx(nv,0) ? _Qset([0,0,0] , ln(nq) ) : + _Qset(_Qvec(q)*atan2(nv,_Qreal(q))/nv, ln(nq)); + + +// Function: Q_pow() +// Usage: +// q2 = Q_pow(q, r); +// Description: +// Returns the quaternion that is the power of the quaternion q to the real exponent r. +// The returned quaternion is normalized if `q` is normalized. +function Q_pow(q,r=1) = + assert( Q_is_quat(q) && is_finite(r), + "Invalid inputs") + let( theta = 2*atan2(norm(_Qvec(q)),_Qreal(q)) ) + Quat(_Qvec(q), r*theta); // Q_exp(r*Q_ln(q)); + + + // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap diff --git a/rounding.scad b/rounding.scad index 6f2a065b..8d6e2c6a 100644 --- a/rounding.scad +++ b/rounding.scad @@ -159,6 +159,7 @@ include // fwd(60) // Note how the different points are cut back by different amounts // stroke(round_corners(zig,radius=1.5,closed=false),width=1); // Example(FlatSpin): Rounding some random 3D paths +// $fn=36; // list1= [ // [2.887360, 4.03497, 6.372090], // [5.682210, 9.37103, 0.783548], @@ -926,10 +927,10 @@ function os_profile(points, extra,check_valid, quality, offset_maxstep, offset) // // Example: Chamfered elliptical prism. If you stretch a chamfered cylinder the chamfer will be uneven. // convex_offset_extrude(bottom = os_chamfer(height=-2), top=os_chamfer(height=1), height=7) -// xscale(4)circle(r=6); +// xscale(4)circle(r=6,$fn=64); // Example: Elliptical prism with circular roundovers. // convex_offset_extrude(bottom=os_circle(r=-2), top=os_circle(r=1), height=7,steps=10) -// xscale(4)circle(r=6); +// xscale(4)circle(r=6,$fn=64); // Example: If you give a non-convex input you get a convex hull output // right(50) linear_extrude(height=7) star(5,r=22,ir=13); // convex_offset_extrude(bottom = os_chamfer(height=-2), top=os_chamfer(height=1), height=7) @@ -1122,6 +1123,7 @@ function _remove_undefined_vals(list) = // fwd(7) // offset_stroke(arc, width=2, start=os_pointed(loc=2,dist=2),end=os_pointed(loc=.5,dist=-1)); // Example(2D): The os_round() end treatment adds roundovers to the end corners by specifying the `cut` parameter. In the first example, the cut parameter is the same at each corner. The bezier smoothness parameter `k` is given to allow a larger cut. In the second example, each corner is given a different roundover, including zero for no rounding at all. The red shows the same strokes without the roundover. +// $fn=36; // arc = arc(points=[[1,1],[3,4],[6,3]],N=50); // path = [[0,0],[6,2],[9,7],[8,10]]; // offset_stroke(path, width=2, rounded=false,start=os_round(angle=-20, cut=0.4,k=.9), end=os_round(angle=-35, cut=0.4,k=.9)); @@ -1133,9 +1135,9 @@ function _remove_undefined_vals(list) = // Example(2D): Negative cut values produce a flaring end. Note how the absolute angle aligns the ends of the first example withi the axes. In the second example positive and negative cut values are combined. Note also that very different cuts are needed at the start end to produce a similar looking flare. // arc = arc(points=[[1,1],[3,4],[6,3]],N=50); // path = [[0,0],[6,2],[9,7],[8,10]]; -// offset_stroke(path, width=2, rounded=false,start=os_round(cut=-1, abs_angle=90), end=os_round(cut=-0.5, abs_angle=0)); +// offset_stroke(path, width=2, rounded=false,start=os_round(cut=-1, abs_angle=90), end=os_round(cut=-0.5, abs_angle=0),$fn=36); // right(10) -// offset_stroke(arc, width=2, rounded=false, start=os_round(cut=[-.75,-.2], angle=-45), end=os_round(cut=[-.2,.2], angle=20)); +// offset_stroke(arc, width=2, rounded=false, start=os_round(cut=[-.75,-.2], angle=-45), end=os_round(cut=[-.2,.2], angle=20),$fn=36); // Example(2D): Setting the width to a vector allows generation of a set of parallel strokes // path = [[0,0],[4,4],[8,4],[2,9],[10,10]]; // for(i=[0:.25:2]) @@ -1146,9 +1148,9 @@ function _remove_undefined_vals(list) = // offset_stroke(path, rounded=true,width = [i,i+.08]); // Example(2D): In this example a spurious triangle appears. This results from overly enthusiastic validity checking. Turning validity checking off fixes it in this case. // path = [[0,0],[4,4],[8,4],[2,9],[10,10]]; -// offset_stroke(path, check_valid=true,rounded=false,width = [1.4, 1.45]); +// offset_stroke(path, check_valid=true,rounded=false,width = [1.4, 1.5]); // right(2) -// offset_stroke(path, check_valid=false,rounded=false,width = [1.4, 1.45]); +// offset_stroke(path, check_valid=false,rounded=false,width = [1.4, 1.5]); // Example(2D): But in this case, disabling the validity check produces an invalid result. // path = [[0,0],[4,4],[8,4],[2,9],[10,10]]; // offset_stroke(path, check_valid=true,rounded=false,width = [1.9, 2]); @@ -1786,7 +1788,7 @@ function _circle_mask(r) = // less than 180 degrees, and the path shouldn't have closely spaced points at concave points of high curvature because // this will cause self-intersection in the mask polyhedron, resulting in CGAL failures. // Arguments: -// r|radius = center radius of the cylindrical shell to cut a hole in +// r / radius = center radius of the cylindrical shell to cut a hole in // thickness = thickness of cylindrical shell (may need to be slighly oversized) // path = 2d path that defines the hole to cut // Example: The mask as long pointed ends because this was the most efficient way to close off those ends. diff --git a/scripts/docs_gen.py b/scripts/docs_gen.py index d2f819c5..b1131d3d 100755 --- a/scripts/docs_gen.py +++ b/scripts/docs_gen.py @@ -97,7 +97,7 @@ def image_compare(file1, file2): sq = (value * (i % 256) ** 2 for i, value in enumerate(diff)) sum_squares = sum(sq) rms = math.sqrt(sum_squares / float(img1.size[0] * img1.size[1])) - return rms<10 + return rms<2 def image_resize(infile, outfile, newsize=(320,240)): diff --git a/shapes2d.scad b/shapes2d.scad index 9a711243..9622ad4f 100644 --- a/shapes2d.scad +++ b/shapes2d.scad @@ -237,8 +237,8 @@ module stroke( } else { quatsums = Q_Cumulative([ for (i = idx(path2,end=-2)) let( - vec1 = i==0? UP : unit(path2[i]-path2[i-1]), - vec2 = unit(path2[i+1]-path2[i]), + vec1 = i==0? UP : unit(path2[i]-path2[i-1], UP), + vec2 = unit(path2[i+1]-path2[i], UP), axis = vector_axis(vec1,vec2), ang = vector_angle(vec1,vec2) ) Quat(axis,ang) @@ -961,8 +961,8 @@ function regular_ngon(n=6, r, d, or, od, ir, id, side, rounding=0, realign=false tipp = polar_to_xy(r-inset+rounding,a1), pos = (p1+p2)/2 ) each [ - anchorpt(str("tip",i), tipp, unit(tipp), 0), - anchorpt(str("side",i), pos, unit(pos), 0), + anchorpt(str("tip",i), tipp, unit(tipp,BACK), 0), + anchorpt(str("side",i), pos, unit(pos,BACK), 0), ] ] ) reorient(anchor,spin, two_d=true, path=path, extent=false, p=path, anchors=anchors); @@ -983,8 +983,8 @@ module regular_ngon(n=6, r, d, or, od, ir, id, side, rounding=0, realign=false, tipp = polar_to_xy(r-inset+rounding,a1), pos = (p1+p2)/2 ) each [ - anchorpt(str("tip",i), tipp, unit(tipp), 0), - anchorpt(str("side",i), pos, unit(pos), 0), + anchorpt(str("tip",i), tipp, unit(tipp,BACK), 0), + anchorpt(str("side",i), pos, unit(pos,BACK), 0), ] ]; attachable(anchor,spin, two_d=true, path=path, extent=false, anchors=anchors) { @@ -1332,9 +1332,9 @@ function star(n, r, d, or, od, ir, id, step, realign=false, anchor=CENTER, spin= p3 = polar_to_xy(r,a3), pos = (p1+p3)/2 ) each [ - anchorpt(str("tip",i), p1, unit(p1), 0), - anchorpt(str("corner",i), p2, unit(p2), 0), - anchorpt(str("midpt",i), pos, unit(pos), 0), + anchorpt(str("tip",i), p1, unit(p1,BACK), 0), + anchorpt(str("corner",i), p2, unit(p2,BACK), 0), + anchorpt(str("midpt",i), pos, unit(pos,BACK), 0), ] ] ) reorient(anchor,spin, two_d=true, path=path, p=path, anchors=anchors); @@ -1355,9 +1355,9 @@ module star(n, r, d, or, od, ir, id, step, realign=false, anchor=CENTER, spin=0) p3 = polar_to_xy(r,a3), pos = (p1+p3)/2 ) each [ - anchorpt(str("tip",i), p1, unit(p1), 0), - anchorpt(str("corner",i), p2, unit(p2), 0), - anchorpt(str("midpt",i), pos, unit(pos), 0), + anchorpt(str("tip",i), p1, unit(p1,BACK), 0), + anchorpt(str("corner",i), p2, unit(p2,BACK), 0), + anchorpt(str("midpt",i), pos, unit(pos,BACK), 0), ] ]; attachable(anchor,spin, two_d=true, path=path, anchors=anchors) { diff --git a/tests/hull.scad b/tests/hull.scad index 98bdb3f5..50366d29 100644 --- a/tests/hull.scad +++ b/tests/hull.scad @@ -1,5 +1,5 @@ -include -include +include <../std.scad> +include <../hull.scad> testpoints_on_sphere = [ for(p = diff --git a/tests/polyhedra.scad b/tests/polyhedra.scad index 5f2bedb2..3a77fbcb 100644 --- a/tests/polyhedra.scad +++ b/tests/polyhedra.scad @@ -1,5 +1,5 @@ -include -include +include<../std.scad> +include<../polyhedra.scad> $fn=96; diff --git a/tests/test_affine.scad b/tests/test_affine.scad index a37153c9..a759b568 100644 --- a/tests/test_affine.scad +++ b/tests/test_affine.scad @@ -1,4 +1,4 @@ -include +include <../std.scad> module test_ident() { diff --git a/tests/test_arrays.scad b/tests/test_arrays.scad index 0ed2bd6d..756e433d 100644 --- a/tests/test_arrays.scad +++ b/tests/test_arrays.scad @@ -1,5 +1,4 @@ -include - +include <../std.scad> // List/Array Ops @@ -16,6 +15,11 @@ module test_in_list() { assert(in_list("bar", ["foo", "bar", "baz"])); assert(!in_list("bee", ["foo", "bar", "baz"])); assert(in_list("bar", [[2,"foo"], [4,"bar"], [3,"baz"]], idx=1)); + assert(!in_list(undef, [3,4,5])); + assert(in_list(undef,[3,4,undef,5])); + assert(!in_list(3,[])); + assert(!in_list(3,[4,5,[3]])); + } test_in_list(); diff --git a/tests/test_common.scad b/tests/test_common.scad index ed97cae7..fa6bdbdc 100644 --- a/tests/test_common.scad +++ b/tests/test_common.scad @@ -1,4 +1,4 @@ -include +include <../std.scad> module test_typeof() { diff --git a/tests/test_coords.scad b/tests/test_coords.scad index 702a3515..4b89c123 100644 --- a/tests/test_coords.scad +++ b/tests/test_coords.scad @@ -1,4 +1,4 @@ -include +include <../std.scad> module test_point2d() { diff --git a/tests/test_cubetruss.scad b/tests/test_cubetruss.scad index 9b71f0ec..3c6965d3 100644 --- a/tests/test_cubetruss.scad +++ b/tests/test_cubetruss.scad @@ -1,5 +1,5 @@ -include -include +include <../std.scad> +include <../cubetruss.scad> module test_cubetruss_dist() { diff --git a/tests/test_debug.scad b/tests/test_debug.scad index 8c305645..b6351c5e 100644 --- a/tests/test_debug.scad +++ b/tests/test_debug.scad @@ -1,4 +1,4 @@ -include +include <../std.scad> module test_standard_anchors() { diff --git a/tests/test_edges.scad b/tests/test_edges.scad index 047c1f9a..0b755d6a 100644 --- a/tests/test_edges.scad +++ b/tests/test_edges.scad @@ -1,4 +1,4 @@ -include +include <../std.scad> module test_is_edge_array() { diff --git a/tests/test_errors.scad b/tests/test_errors.scad index ba82755b..fffe9ebb 100644 --- a/tests/test_errors.scad +++ b/tests/test_errors.scad @@ -1,4 +1,4 @@ -include +include <../std.scad> // Can't test echo output as yet. Include these for coverage calculations. diff --git a/tests/test_geometry.scad b/tests/test_geometry.scad index 4b6dd90f..78c3e575 100644 --- a/tests/test_geometry.scad +++ b/tests/test_geometry.scad @@ -1,4 +1,4 @@ -include +include <../std.scad> module test_point_on_segment2d() { @@ -87,7 +87,7 @@ module test_line_normal() { assert(line_normal([[0,0],[-10,0]]) == [0,-1]); assert(line_normal([[0,0],[0,-10]]) == [1,0]); assert(approx(line_normal([[0,0],[10,10]]), [-sqrt(2)/2,sqrt(2)/2])); - pts = [for (i=list_range(1000)) rands(-100,100,2,seed_value=4312)]; + pts = [for (p=pair(rands(-100,100,1000,seed_value=4312))) p]; for (p = pair_wrap(pts)) { p1 = p.x; p2 = p.y; @@ -351,18 +351,18 @@ module test_tri_functions() { opp = p.y; hyp = norm([opp,adj]); ang = atan2(opp,adj); - assert(approx(hyp_opp_to_adj(hyp,opp),adj)); - assert(approx(hyp_ang_to_adj(hyp,ang),adj)); - assert(approx(opp_ang_to_adj(opp,ang),adj)); - assert(approx(hyp_adj_to_opp(hyp,adj),opp)); - assert(approx(hyp_ang_to_opp(hyp,ang),opp)); - assert(approx(adj_ang_to_opp(adj,ang),opp)); - assert(approx(adj_opp_to_hyp(adj,opp),hyp)); - assert(approx(adj_ang_to_hyp(adj,ang),hyp)); - assert(approx(opp_ang_to_hyp(opp,ang),hyp)); - assert(approx(hyp_adj_to_ang(hyp,adj),ang)); - assert(approx(hyp_opp_to_ang(hyp,opp),ang)); - assert(approx(adj_opp_to_ang(adj,opp),ang)); + assert_approx(hyp_opp_to_adj(hyp,opp), adj); + assert_approx(hyp_ang_to_adj(hyp,ang), adj); + assert_approx(opp_ang_to_adj(opp,ang), adj); + assert_approx(hyp_adj_to_opp(hyp,adj), opp); + assert_approx(hyp_ang_to_opp(hyp,ang), opp); + assert_approx(adj_ang_to_opp(adj,ang), opp); + assert_approx(adj_opp_to_hyp(adj,opp), hyp); + assert_approx(adj_ang_to_hyp(adj,ang), hyp); + assert_approx(opp_ang_to_hyp(opp,ang), hyp); + assert_approx(hyp_adj_to_ang(hyp,adj), ang); + assert_approx(hyp_opp_to_ang(hyp,opp), ang); + assert_approx(adj_opp_to_ang(adj,opp), ang); } } test_tri_functions(); @@ -613,6 +613,23 @@ module test_pointlist_bounds() { [23,57,-42] ]; assert(pointlist_bounds(pts) == [[-63,-32,-42], [84,97,42]]); + pts2d = [ + [-53,12], + [-63,36], + [84,-5], + [63,42], + [23,-42] + ]; + assert(pointlist_bounds(pts2d) == [[-63,-42],[84,42]]); + pts5d = [ + [-53,27,12,-53,12], + [-63,97,36,-63,36], + [84,-32,-5,84,-5], + [63,-24,42,63,42], + [23,57,-42,23,-42] + ]; + assert(pointlist_bounds(pts5d) == [[-63,-32,-42,-63,-42],[84,97,42,84,42]]); + assert(pointlist_bounds([[3,4,5,6]]), [[3,4,5,6],[3,4,5,6]]); } test_pointlist_bounds(); diff --git a/tests/test_linear_bearings.scad b/tests/test_linear_bearings.scad index 19bc2f50..fb57ba53 100644 --- a/tests/test_linear_bearings.scad +++ b/tests/test_linear_bearings.scad @@ -1,5 +1,5 @@ -include -include +include <../std.scad> +include <../linear_bearings.scad> module test_get_lmXuu_bearing_diam() { diff --git a/tests/test_math.scad b/tests/test_math.scad index 16f49882..3233e6ff 100644 --- a/tests/test_math.scad +++ b/tests/test_math.scad @@ -1,4 +1,4 @@ -include +include <../std.scad> // Simple Calculations @@ -76,6 +76,12 @@ module test_is_matrix() { assert(is_matrix([[2,3,4],[5,6,7],[8,9,10]],square=false)); assert(is_matrix([[2,3],[5,6],[8,9]],m=3,n=2)); assert(is_matrix([[2,3,4],[5,6,7]],m=2,n=3)); + assert(is_matrix([[2,3,4],[5,6,7]],2,3)); + assert(is_matrix([[2,3,4],[5,6,7]],m=2)); + assert(is_matrix([[2,3,4],[5,6,7]],2)); + assert(is_matrix([[2,3,4],[5,6,7]],n=3)); + assert(!is_matrix([[2,3,4],[5,6,7]],m=4)); + assert(!is_matrix([[2,3,4],[5,6,7]],n=5)); assert(!is_matrix([[2,3,4],[5,6,7]],m=2,n=3,square=true)); assert(is_matrix([[2,3,4],[5,6,7],[8,9,10]],square=false)); assert(!is_matrix([[2,3],[5,6],[8,9]],m=2,n=3)); @@ -181,6 +187,9 @@ module test_sqr() { assert_equal(sqr(2.5), 6.25); assert_equal(sqr(3), 9); assert_equal(sqr(16), 256); + assert_equal(sqr([2,3,4]), [4,9,16]); + assert_equal(sqr([[2,3,4],[3,5,7]]), [[4,9,16],[9,25,49]]); + assert_equal(sqr([]),[]); } test_sqr(); @@ -525,6 +534,9 @@ module test_any() { assert_equal(any([1,5,true]), true); assert_equal(any([[0,0], [0,0]]), false); assert_equal(any([[0,0], [1,0]]), true); + assert_equal(any([[false,false],[[false,[false],[[[true]]]],false],[false,false]]), true); + assert_equal(any([[false,false],[[false,[false],[[[false]]]],false],[false,false]]), false); + assert_equal(any([]), false); } test_any(); @@ -536,6 +548,9 @@ module test_all() { assert_equal(all([[0,0], [0,0]]), false); assert_equal(all([[0,0], [1,0]]), false); assert_equal(all([[1,1], [1,1]]), true); + assert_equal(all([[true,true],[[true,[true],[[[true]]]],true],[true,true]]), true); + assert_equal(all([[true,true],[[true,[true],[[[false]]]],true],[true,true]]), false); + assert_equal(all([]), true); } test_all(); @@ -554,6 +569,7 @@ test_count_true(); module test_factorial() { + assert_equal(factorial(0), 1); assert_equal(factorial(1), 1); assert_equal(factorial(2), 2); assert_equal(factorial(3), 6); @@ -562,6 +578,8 @@ module test_factorial() { assert_equal(factorial(6), 720); assert_equal(factorial(7), 5040); assert_equal(factorial(8), 40320); + assert_equal(factorial(25,21), 303600); + assert_equal(factorial(25,25), 1); } test_factorial(); @@ -570,6 +588,11 @@ module test_gcd() { assert_equal(gcd(15,25), 5); assert_equal(gcd(15,27), 3); assert_equal(gcd(270,405), 135); + assert_equal(gcd(39, 101),1); + assert_equal(gcd(15,-25), 5); + assert_equal(gcd(-15,25), 5); + assert_equal(gcd(5,0),5); + assert_equal(gcd(0,5),5); } test_gcd(); @@ -578,9 +601,306 @@ module test_lcm() { assert_equal(lcm(15,25), 75); assert_equal(lcm(15,27), 135); assert_equal(lcm(270,405), 810); + assert_equal(lcm([3,5,15,25,35]),525); } 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]); +} +test_C_times(); + + +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(); + + +module test_back_substitute(){ + R = [[12,4,3,2], + [0,2,-4,2], + [0,0,4,5], + [0,0,0,15]]; + assert_approx(back_substitute(R, [1,2,3,3]), [-0.675, 1.8, 0.5, 0.2]); + assert_approx(back_substitute(R, [6, 3, 3.5, 37], transpose=true), [0.5, 0.5, 1, 2]); + assert_approx(back_substitute(R, [[38,101],[-6,-16], [31, 71], [45, 105]]), [[1, 4],[2,3],[4,9],[3,7]]); + assert_approx(back_substitute(R, [[12,48],[8,22],[11,36],[71,164]],transpose=true), [[1, 4],[2,3],[4,9],[3,7]]); + assert_approx(back_substitute([[2]], [4]), [2]); + sing1 =[[0,4,3,2], + [0,3,-4,2], + [0,0,4,5], + [0,0,0,15]]; + sing2 =[[12,4,3,2], + [0,0,-4,2], + [0,0,4,5], + [0,0,0,15]]; + sing3 = [[12,4,3,2], + [0,2,-4,2], + [0,0,4,5], + [0,0,0,0]]; + assert_approx(back_substitute(sing1, [1,2,3,4]), []); + assert_approx(back_substitute(sing2, [1,2,3,4]), []); + assert_approx(back_substitute(sing3, [1,2,3,4]), []); +} +test_back_substitute(); + + + +module test_linear_solve(){ + M = [[-2,-5,-1,3], + [3,7,6,2], + [6,5,-1,-6], + [-7,1,2,3]]; + assert_approx(linear_solve(M, [-3,43,-11,13]), [1,2,3,4]); + assert_approx(linear_solve(M, [[-5,8],[18,-61],[4,7],[-1,-12]]), [[1,-2],[1,-3],[1,-4],[1,-5]]); + assert_approx(linear_solve([[2]],[4]), [2]); + assert_approx(linear_solve([[2]],[[4,8]]), [[2, 4]]); + assert_approx(linear_solve(select(M,0,2), [2,4,4]), [ 2.254871220604705e+00, + -8.378819388897780e-01, + 2.330507118860985e-01, + 8.511278195488737e-01]); + assert_approx(linear_solve(subindex(M,[0:2]), [2,4,4,4]), + [-2.457142857142859e-01, + 5.200000000000000e-01, + 7.428571428571396e-02]); + assert_approx(linear_solve([[1,2,3,4]], [2]), [0.066666666666666, 0.13333333333, 0.2, 0.266666666666]); + assert_approx(linear_solve([[1],[2],[3],[4]], [4,3,2,1]), [2/3]); + rd = [[-2,-5,-1,3], + [3,7,6,2], + [3,7,6,2], + [-7,1,2,3]]; + assert_equal(linear_solve(rd,[1,2,3,4]),[]); + assert_equal(linear_solve(select(rd,0,2), [2,4,4]), []); + assert_equal(linear_solve(transpose(select(rd,0,2)), [2,4,3,4]), []); +} +test_linear_solve(); + + +module test_outer_product(){ + assert_equal(outer_product([1,2,3],[4,5,6]), [[4,5,6],[8,10,12],[12,15,18]]); + assert_equal(outer_product([9],[7]), [[63]]); +} +test_outer_product(); + + +module test_deriv(){ + pent = [for(x=[0:70:359]) [cos(x), sin(x)]]; + assert_approx(deriv(pent,closed=true), + [[-0.321393804843,0.556670399226], + [-0.883022221559,0.321393804843], + [-0.604022773555,-0.719846310393], + [0.469846310393,-0.813797681349], + [0.925416578398,0.163175911167], + [0.413175911167,0.492403876506]]); + assert_approx(deriv(pent,closed=true,h=2), + 0.5*[[-0.321393804843,0.556670399226], + [-0.883022221559,0.321393804843], + [-0.604022773555,-0.719846310393], + [0.469846310393,-0.813797681349], + [0.925416578398,0.163175911167], + [0.413175911167,0.492403876506]]); + assert_approx(deriv(pent,closed=false), + [[-0.432937491789,1.55799143673], + [-0.883022221559,0.321393804843], + [-0.604022773555,-0.719846310393], + [0.469846310393,-0.813797681349], + [0.925416578398,0.163175911167], + [0.696902572292,1.45914323952]]); + spent = yscale(8,pent); + lens = path_segment_lengths(spent,closed=true); + assert_approx(deriv(spent, closed=true, h=lens), + [[-0.0381285841663,0.998065839726], + [-0.254979378104,0.0449763331253], + [-0.216850793938,-0.953089506601], + [0.123993253223,-0.982919228715], + [0.191478335034,0.0131898128456], + [0.0674850818111,0.996109041561]]); + assert_approx(deriv(spent, closed=false, h=select(lens,0,-2)), + [[-0.0871925973657,0.996191473044], + [-0.254979378104,0.0449763331253], + [-0.216850793938,-0.953089506601], + [0.123993253223,-0.982919228715], + [0.191478335034,0.0131898128456], + [0.124034734589,0.992277876714]]); +} +test_deriv(); + + +module test_deriv2(){ + oct = [for(x=[0:45:359]) [cos(x), sin(x)]]; + assert_approx(deriv2(oct), + [[-0.828427124746,0.0719095841794],[-0.414213562373,-0.414213562373],[0,-0.585786437627], + [0.414213562373,-0.414213562373],[0.585786437627,0],[0.414213562373,0.414213562373], + [0,0.585786437627],[-0.636634192232,0.534938683021]]); + assert_approx(deriv2(oct,closed=false), + [[-0.828427124746,0.0719095841794],[-0.414213562373,-0.414213562373],[0,-0.585786437627], + [0.414213562373,-0.414213562373],[0.585786437627,0],[0.414213562373,0.414213562373], + [0,0.585786437627],[-0.636634192232,0.534938683021]]); + assert_approx(deriv2(oct,closed=true), + [[-0.585786437627,0],[-0.414213562373,-0.414213562373],[0,-0.585786437627], + [0.414213562373,-0.414213562373],[0.585786437627,0],[0.414213562373,0.414213562373], + [0,0.585786437627],[-0.414213562373,0.414213562373]]); + assert_approx(deriv2(oct,closed=false,h=2), + 0.25*[[-0.828427124746,0.0719095841794],[-0.414213562373,-0.414213562373],[0,-0.585786437627], + [0.414213562373,-0.414213562373],[0.585786437627,0],[0.414213562373,0.414213562373], + [0,0.585786437627],[-0.636634192232,0.534938683021]]); + assert_approx(deriv2(oct,closed=true,h=2), + 0.25* [[-0.585786437627,0],[-0.414213562373,-0.414213562373],[0,-0.585786437627], + [0.414213562373,-0.414213562373],[0.585786437627,0],[0.414213562373,0.414213562373], + [0,0.585786437627],[-0.414213562373,0.414213562373]]); +} +test_deriv2(); + + +module test_deriv3(){ + oct = [for(x=[0:45:359]) [cos(x), sin(x)]]; + assert_approx(deriv3(oct), + [[0.414213562373,-0.686291501015],[0.414213562373,-0.343145750508],[0.414213562373,0], + [0.292893218813,0.292893218813],[0,0.414213562373],[-0.292893218813,0.292893218813], + [-0.535533905933,0.0502525316942],[-0.778174593052,-0.192388155425]]); + assert_approx(deriv3(oct,closed=false), + [[0.414213562373,-0.686291501015],[0.414213562373,-0.343145750508],[0.414213562373,0], + [0.292893218813,0.292893218813],[0,0.414213562373],[-0.292893218813,0.292893218813], + [-0.535533905933,0.0502525316942],[-0.778174593052,-0.192388155425]]); + assert_approx(deriv3(oct,closed=false,h=2), + [[0.414213562373,-0.686291501015],[0.414213562373,-0.343145750508],[0.414213562373,0], + [0.292893218813,0.292893218813],[0,0.414213562373],[-0.292893218813,0.292893218813], + [-0.535533905933,0.0502525316942],[-0.778174593052,-0.192388155425]]/8); + assert_approx(deriv3(oct,closed=true), + [[0,-0.414213562373],[0.292893218813,-0.292893218813],[0.414213562373,0],[0.292893218813,0.292893218813], + [0,0.414213562373],[-0.292893218813,0.292893218813],[-0.414213562373,0],[-0.292893218813,-0.292893218813]]); + assert_approx(deriv3(oct,closed=true,h=2), + [[0,-0.414213562373],[0.292893218813,-0.292893218813],[0.414213562373,0],[0.292893218813,0.292893218813], + [0,0.414213562373],[-0.292893218813,0.292893218813],[-0.414213562373,0],[-0.292893218813,-0.292893218813]]/8); +} +test_deriv3(); + + + +module test_polynomial(){ + assert_equal(polynomial([],12),0); + assert_equal(polynomial([],[12,4]),[0,0]); + assert_equal(polynomial([1,2,3,4],3),58); + assert_equal(polynomial([1,2,3,4],[3,-1]),[47,-41]); + assert_equal(polynomial([0,0,2],4),2); +} +test_polynomial(); + + +module test_poly_roots(){ + // Fifth roots of unity + assert_approx( + poly_roots([1,0,0,0,0,-1]), + [[1,0],[0.309016994375,0.951056516295],[-0.809016994375,0.587785252292], + [-0.809016994375,-0.587785252292],[0.309016994375,-0.951056516295]]); + assert_approx(poly_roots(poly_mult([[1,-2,5],[12,-24,24],[-2, -12, -20],[1,-10,50]])), + [[1, 1], [5, 5], [1, 2], [-3, 1], [-3, -1], [1, -1], [1, -2], [5, -5]]); + assert_approx(poly_roots([.124,.231,.942, -.334]), + [[0.3242874219074053,0],[-1.093595323856930,2.666477428660098], [-1.093595323856930,-2.666477428660098]]); +} +test_poly_roots(); + +module test_real_roots(){ + // Wilkinson polynomial is a nasty test: + assert_approx( + sort(real_roots(poly_mult([[1,-1],[1,-2],[1,-3],[1,-4],[1,-5],[1,-6],[1,-7],[1,-8],[1,-9],[1,-10]]))), + list_range(n=10,s=1)); + assert_equal(real_roots([3]), []); + assert_equal(real_roots(poly_mult([[1,-2,5],[12,-24,24],[-2, -12, -20],[1,-10,50]])),[]); + assert_equal(real_roots(poly_mult([[1,-2,5],[12,-24,24],[-2, -12, -20],[1,-10,50],[1,0,0]])),[0,0]); + assert_approx(real_roots(poly_mult([[1,-2,5],[12,-24,24],[-2, -12, -20],[1,-10,50],[1,4]])),[-4]); + assert(approx(real_roots([1,-10,25]),[5,5],eps=5e-6)); + assert_approx(real_roots([4,-3]), [0.75]); + assert_approx(real_roots([0,0,0,4,-3]), [0.75]); +} +test_real_roots(); + +// Need decision about behavior for out of bounds ranges, empty ranges +module test_submatrix(){ + M = [[1,2,3,4,5], + [6,7,8,9,10], + [11,12,13,14,15], + [16,17,18,19,20], + [21,22,23,24,25]]; + assert_equal(submatrix(M,[1:2], [3:4]), [[9,10],[14,15]]); + assert_equal(submatrix(M,[1], [3,4]), [[9,10]]); + assert_equal(submatrix(M,1, [3,4]), [[9,10]]); + assert_equal(submatrix(M, [3,4],1), [[17],[22]]); + assert_equal(submatrix(M, [1,3],[2,4]), [[8,10],[18,20]]); +} +test_submatrix(); + + + +module test_qr_factor() { + // Check that R is upper triangular + function is_ut(R) = + let(bad = [for(i=[1:1:len(R)-1], j=[0:min(i-1, len(R[0])-1)]) if (!approx(R[i][j],0)) 1]) + bad == []; + + // Test the R is upper trianglar, Q is orthogonal and qr=M + function qrok(qr,M) = + is_ut(qr[1]) && approx(qr[0]*transpose(qr[0]), ident(len(qr[0]))) && approx(qr[0]*qr[1],M); + + M = [[1,2,9,4,5], + [6,7,8,19,10], + [11,12,13,14,15], + [1,17,18,19,20], + [21,22,10,24,25]]; + + assert(qrok(qr_factor(M),M)); + assert(qrok(qr_factor(select(M,0,3)),select(M,0,3))); + assert(qrok(qr_factor(transpose(select(M,0,3))),transpose(select(M,0,3)))); + + A = [[1,2,9,4,5], + [6,7,8,19,10], + [0,0,0,0,0], + [1,17,18,19,20], + [21,22,10,24,25]]; + assert(qrok(qr_factor(A),A)); + + B = [[1,2,0,4,5], + [6,7,0,19,10], + [0,0,0,0,0], + [1,17,0,19,20], + [21,22,0,24,25]]; + + assert(qrok(qr_factor(B),B)); + assert(qrok(qr_factor([[7]]), [[7]])); + assert(qrok(qr_factor([[1,2,3]]), [[1,2,3]])); + assert(qrok(qr_factor([[1],[2],[3]]), [[1],[2],[3]])); +} +test_qr_factor(); + + +module test_poly_mult(){ + assert_equal(poly_mult([3,2,1],[4,5,6,7]),[12,23,32,38,20,7]); + assert_equal(poly_mult([3,2,1],[]),[]); + assert_equal(poly_mult([[1,2],[3,4],[5,6]]), [15,68,100,48]); + assert_equal(poly_mult([[1,2],[],[5,6]]), []); + assert_equal(poly_mult([[3,4,5],[0,0,0]]),[]); +} +test_poly_mult(); + + +module test_poly_div(){ + assert_equal(poly_div(poly_mult([4,3,3,2],[2,1,3]), [2,1,3]),[[4,3,3,2],[]]); + assert_equal(poly_div([1,2,3,4],[1,2,3,4,5]), [[], [1,2,3,4]]); + assert_equal(poly_div(poly_add(poly_mult([1,2,3,4],[2,0,2]), [1,1,2]), [1,2,3,4]), [[2,0,2],[1,1,2]]); + assert_equal(poly_div([1,2,3,4], [1,-3]), [[1,5,18],[58]]); +} +test_poly_div(); + + +module test_poly_add(){ + assert_equal(poly_add([2,3,4],[3,4,5,6]),[3,6,8,10]); + assert_equal(poly_add([1,2,3,4],[-1,-2,3,4]), [6,8]); + assert_equal(poly_add([1,2,3],-[1,2,3]),[]); +} +test_poly_add(); // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap diff --git a/tests/test_mutators.scad b/tests/test_mutators.scad index 5854b529..73614344 100644 --- a/tests/test_mutators.scad +++ b/tests/test_mutators.scad @@ -1,4 +1,4 @@ -include +include <../std.scad> module test_HSL() { diff --git a/tests/test_primitives.scad b/tests/test_primitives.scad index 8780b5be..e903617a 100644 --- a/tests/test_primitives.scad +++ b/tests/test_primitives.scad @@ -1,4 +1,4 @@ -include +include <../std.scad> module test_square() { diff --git a/tests/test_quaternions.scad b/tests/test_quaternions.scad index b9b232fa..b1e51309 100644 --- a/tests/test_quaternions.scad +++ b/tests/test_quaternions.scad @@ -1,5 +1,5 @@ -include -include +include <../std.scad> +include <../strings.scad> function rec_cmp(a,b,eps=1e-9) = @@ -8,6 +8,7 @@ function rec_cmp(a,b,eps=1e-9) = is_list(a)? len(a)==len(b) && all([for (i=idx(a)) rec_cmp(a[i],b[i],eps=eps)]) : a == b; +function Qstandard(q) = sign([for(qi=q) if( ! approx(qi,0)) qi,0 ][0])*q; module verify_f(actual,expected) { if (!rec_cmp(actual,expected)) { @@ -22,6 +23,15 @@ module verify_f(actual,expected) { } +module test_is_quat() { + verify_f(Q_is_quat([0]),false); + verify_f(Q_is_quat([0,0,0,0]),false); + verify_f(Q_is_quat([1,0,2,0]),true); + verify_f(Q_is_quat([1,0,2,0,0]),false); +} +test_is_quat(); + + module test_Quat() { verify_f(Quat(UP,0),[0,0,0,1]); verify_f(Quat(FWD,0),[0,0,0,1]); @@ -92,6 +102,15 @@ module test_QuatXYZ() { test_QuatXYZ(); +module test_Q_From_to() { + verify_f(Q_Mul(Q_From_to([1,2,3], [4,5,2]),Q_From_to([4,5,2], [1,2,3])), Q_Ident()); + verify_f(Q_Matrix4(Q_From_to([1,2,3], [4,5,2])), rot(from=[1,2,3],to=[4,5,2])); + verify_f(Qrot(Q_From_to([1,2,3], -[1,2,3]),[1,2,3]), -[1,2,3]); + verify_f(unit(Qrot(Q_From_to([1,2,3], [4,5,2]),[1,2,3])), unit([4,5,2])); +} +test_Q_From_to(); + + module test_Q_Ident() { verify_f(Q_Ident(), [0,0,0,1]); } @@ -207,6 +226,16 @@ module test_Q_Conj() { test_Q_Conj(); +module test_Q_Inverse() { + + verify_f(Q_Inverse([1,0,0,1]),[-1,0,0,1]/sqrt(2)); + verify_f(Q_Inverse([0,1,1,0]),[0,-1,-1,0]/sqrt(2)); + verify_f(Q_Inverse(QuatXYZ([23,45,67])),Q_Conj(QuatXYZ([23,45,67]))); + verify_f(Q_Mul(Q_Inverse(QuatXYZ([23,45,67])),QuatXYZ([23,45,67])),Q_Ident()); +} +test_Q_Inverse(); + + module test_Q_Norm() { verify_f(Q_Norm([1,0,0,1]),1.414213562); verify_f(Q_Norm([0,1,1,0]),1.414213562); @@ -276,6 +305,10 @@ module test_Q_Angle() { verify_f(Q_Angle(QuatY(-37)),37); verify_f(Q_Angle(QuatZ(37)),37); verify_f(Q_Angle(QuatZ(-37)),37); + + verify_f(Q_Angle(QuatZ(-37),QuatZ(-37)), 0); + verify_f(Q_Angle(QuatZ( 37.123),QuatZ(-37.123)), 74.246); + verify_f(Q_Angle(QuatX( 37),QuatY(-37)), 51.86293283); } test_Q_Angle(); @@ -288,4 +321,87 @@ module test_Qrot() { test_Qrot(); +module test_Q_Rotation() { + verify_f(Qstandard(Q_Rotation(Q_Matrix3(Quat([12,34,56],33)))),Qstandard(Quat([12,34,56],33))); + verify_f(Q_Matrix3(Q_Rotation(Q_Matrix3(QuatXYZ([12,34,56])))), + Q_Matrix3(QuatXYZ([12,34,56]))); +} +test_Q_Rotation(); + + +module test_Q_Rotation_path() { + + verify_f(Q_Rotation_path(QuatX(135), 5, QuatY(13.5))[0] , Q_Matrix4(QuatX(135))); + verify_f(Q_Rotation_path(QuatX(135), 11, QuatY(13.5))[11] , yrot(13.5)); + verify_f(Q_Rotation_path(QuatX(135), 16, QuatY(13.5))[8] , Q_Rotation_path(QuatX(135), 8, QuatY(13.5))[4]); + verify_f(Q_Rotation_path(QuatX(135), 16, QuatY(13.5))[7] , + Q_Rotation_path(QuatY(13.5),16, QuatX(135))[9]); + + verify_f(Q_Rotation_path(QuatX(11), 5)[0] , xrot(11)); + verify_f(Q_Rotation_path(QuatX(11), 5)[4] , xrot(55)); + +} +test_Q_Rotation_path(); + + +module test_Q_Nlerp() { + verify_f(Q_Nlerp(QuatX(45),QuatY(30),0.0),QuatX(45)); + verify_f(Q_Nlerp(QuatX(45),QuatY(30),0.5),[0.1967063121, 0.1330377423, 0, 0.9713946602]); + verify_f(Q_Rotation_path(QuatX(135), 16, QuatY(13.5))[8] , Q_Matrix4(Q_Nlerp(QuatX(135), QuatY(13.5),0.5))); + verify_f(Q_Nlerp(QuatX(45),QuatY(30),1.0),QuatY(30)); +} +test_Q_Nlerp(); + + +module test_Q_Squad() { + verify_f(Q_Squad(QuatX(45),QuatZ(30),QuatX(90),QuatY(30),0.0),QuatX(45)); + verify_f(Q_Squad(QuatX(45),QuatZ(30),QuatX(90),QuatY(30),1.0),QuatY(30)); + verify_f(Q_Squad(QuatX(0),QuatX(30),QuatX(90),QuatX(120),0.5), + Q_Slerp(QuatX(0),QuatX(120),0.5)); + verify_f(Q_Squad(QuatY(0),QuatY(0),QuatX(120),QuatX(120),0.3), + Q_Slerp(QuatY(0),QuatX(120),0.3)); +} +test_Q_Squad(); + + +module test_Q_exp() { + verify_f(Q_exp(Q_Ident()), exp(1)*Q_Ident()); + verify_f(Q_exp([0,0,0,33.7]), exp(33.7)*Q_Ident()); + verify_f(Q_exp(Q_ln(Q_Ident())), Q_Ident()); + verify_f(Q_exp(Q_ln([1,2,3,0])), [1,2,3,0]); + verify_f(Q_exp(Q_ln(QuatXYZ([31,27,34]))), QuatXYZ([31,27,34])); + let(q=QuatXYZ([12,23,34])) + verify_f(Q_exp(q+Q_Inverse(q)),Q_Mul(Q_exp(q),Q_exp(Q_Inverse(q)))); + +} +test_Q_exp(); + + +module test_Q_ln() { + verify_f(Q_ln([1,2,3,0]), [24.0535117721, 48.1070235442, 72.1605353164, 1.31952866481]); + verify_f(Q_ln(Q_Ident()), [0,0,0,0]); + verify_f(Q_ln(5.5*Q_Ident()), [0,0,0,ln(5.5)]); + verify_f(Q_ln(Q_exp(QuatXYZ([13,37,43]))), QuatXYZ([13,37,43])); + verify_f(Q_ln(QuatXYZ([12,23,34]))+Q_ln(Q_Inverse(QuatXYZ([12,23,34]))), [0,0,0,0]); +} +test_Q_ln(); + + +module test_Q_pow() { + q = Quat([1,2,3],77); + verify_f(Q_pow(q,1), q); + verify_f(Q_pow(q,0), Q_Ident()); + verify_f(Q_pow(q,-1), Q_Inverse(q)); + verify_f(Q_pow(q,2), Q_Mul(q,q)); + verify_f(Q_pow(q,3), Q_Mul(q,Q_pow(q,2))); + verify_f(Q_Mul(Q_pow(q,0.456),Q_pow(q,0.544)), q); + verify_f(Q_Mul(Q_pow(q,0.335),Q_Mul(Q_pow(q,.552),Q_pow(q,.113))), q); +} +test_Q_pow(); + + + + + + // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap diff --git a/tests/test_queues.scad b/tests/test_queues.scad index cf45eefd..4ce3c8ce 100644 --- a/tests/test_queues.scad +++ b/tests/test_queues.scad @@ -1,5 +1,5 @@ -include -include +include <../std.scad> +include <../queues.scad> module test_queue_init() { diff --git a/tests/test_shapes.scad b/tests/test_shapes.scad index 0b66d60d..ec1d58be 100644 --- a/tests/test_shapes.scad +++ b/tests/test_shapes.scad @@ -1,5 +1,5 @@ -include -include +include <../std.scad> +include <../hull.scad> module test_prismoid() { diff --git a/tests/test_shapes2d.scad b/tests/test_shapes2d.scad index dc8e874b..1ee5ffa6 100644 --- a/tests/test_shapes2d.scad +++ b/tests/test_shapes2d.scad @@ -1,4 +1,4 @@ -include +include <../std.scad> module test_turtle() { diff --git a/tests/test_skin.scad b/tests/test_skin.scad index 347ccd23..c4a231d2 100644 --- a/tests/test_skin.scad +++ b/tests/test_skin.scad @@ -1,5 +1,5 @@ -include -include +include <../std.scad> +include <../skin.scad> module test_skin() { diff --git a/tests/test_stacks.scad b/tests/test_stacks.scad index 3f3f0d21..72412035 100644 --- a/tests/test_stacks.scad +++ b/tests/test_stacks.scad @@ -1,5 +1,5 @@ -include -include +include <../std.scad> +include <../stacks.scad> module test_stack_init() { diff --git a/tests/test_strings.scad b/tests/test_strings.scad index 673e8f08..a1d05d6a 100644 --- a/tests/test_strings.scad +++ b/tests/test_strings.scad @@ -1,5 +1,5 @@ -include -include +include <../std.scad> +include <../strings.scad> module test_upcase() { diff --git a/tests/test_structs.scad b/tests/test_structs.scad index 55db3b45..a4a4d924 100644 --- a/tests/test_structs.scad +++ b/tests/test_structs.scad @@ -1,5 +1,5 @@ -include -include +include <../std.scad> +include <../structs.scad> module test_struct_set() { diff --git a/tests/test_transforms.scad b/tests/test_transforms.scad index 722b65e7..a44e47cc 100644 --- a/tests/test_transforms.scad +++ b/tests/test_transforms.scad @@ -1,4 +1,4 @@ -include +include <../std.scad> module test_translate() { diff --git a/tests/test_vectors.scad b/tests/test_vectors.scad index 399be016..42f0a04c 100644 --- a/tests/test_vectors.scad +++ b/tests/test_vectors.scad @@ -1,4 +1,4 @@ -include +include <../std.scad> module test_is_vector() { @@ -9,6 +9,10 @@ module test_is_vector() { assert(is_vector(1) == false); assert(is_vector("foo") == false); assert(is_vector(true) == false); + assert(is_vector([0,0,0],zero=true) == true); + assert(is_vector([0,0,0],zero=false) == false); + assert(is_vector([0,1,0],zero=true) == false); + assert(is_vector([0,0,1],zero=false) == true); } test_is_vector(); @@ -56,7 +60,7 @@ module test_vabs() { } test_vabs(); -include +include <../strings.scad> module test_vang() { assert(vang([1,0])==0); assert(vang([0,1])==90); diff --git a/tests/test_version.scad b/tests/test_version.scad index f877b5ce..e82f13f2 100644 --- a/tests/test_version.scad +++ b/tests/test_version.scad @@ -1,4 +1,4 @@ -include +include <../std.scad> module test_bosl_version() { diff --git a/tests/test_vnf.scad b/tests/test_vnf.scad index 65d53503..753803cb 100644 --- a/tests/test_vnf.scad +++ b/tests/test_vnf.scad @@ -1,5 +1,5 @@ -include -include +include <../std.scad> +include <../vnf.scad> module test_is_vnf() { diff --git a/vectors.scad b/vectors.scad index 1ee1f5d1..bb28ead9 100644 --- a/vectors.scad +++ b/vectors.scad @@ -19,6 +19,8 @@ // Arguments: // v = The value to test to see if it is a vector. // length = If given, make sure the vector is `length` items long. +// zero = If false, require that the length of the vector is not approximately zero. If true, require the length of the vector to be approx zero-length. Default: `undef` (don't check vector length.) +// eps = The minimum vector length that is considered non-zero. Default: `EPSILON` (`1e-9`) // Example: // is_vector(4); // Returns false // is_vector([4,true,false]); // Returns false @@ -28,8 +30,14 @@ // is_vector([3,4,5],3); // Returns true // is_vector([3,4,5],4); // Returns true // is_vector([]); // Returns false -function is_vector(v,length) = - is_list(v) && is_num(0*(v*v)) && (is_undef(length)||len(v)==length); +// is_vector([0,0,0],zero=true); // Returns true +// is_vector([0,0,0],zero=false); // Returns false +// is_vector([0,1,0],zero=true); // Returns false +// is_vector([0,0,1],zero=false); // Returns true +function is_vector(v,length,zero,eps=EPSILON) = + is_list(v) && is_num(0*(v*v)) + && (is_undef(length) || len(v)==length) + && (is_undef(zero) || ((norm(v) >= eps) == !zero)); // Function: add_scalar() @@ -105,18 +113,25 @@ function vceil(v) = [for (x=v) ceil(x)]; // Function: unit() +// Usage: +// unit(v, [error]); // Description: -// Returns unit length normalized version of vector v. -// If passed a zero-length vector, returns the unchanged vector. +// Returns the unit length normalized version of vector v. If passed a zero-length vector, +// asserts an error unless `error` is given, in which case the value of `error` is returned. // Arguments: // v = The vector to normalize. +// error = If given, and input is a zero-length vector, this value is returned. Default: Assert error on zero-length vector. // Examples: // unit([10,0,0]); // Returns: [1,0,0] // unit([0,10,0]); // Returns: [0,1,0] // unit([0,0,10]); // Returns: [0,0,1] // unit([0,-10,0]); // Returns: [0,-1,0] -// unit([0,0,0]); // Returns: [0,0,0] -function unit(v) = assert(is_vector(v),str(v)) norm(v)<=EPSILON? v : v/norm(v); +// unit([0,0,0],[1,2,3]); // Returns: [1,2,3] +// unit([0,0,0]); // Asserts an error. +function unit(v, error=[[["ASSERT"]]]) = + assert(is_vector(v), str("Expected a vector. Got: ",v)) + norm(v)=EPSILON) : error) : + v/norm(v); // Function: vector_angle() @@ -181,22 +196,31 @@ function vector_angle(v1,v2,v3) = // vector_axis([10,0,10], [0,0,0], [-10,10,0]); // Returns: [-0.57735, -0.57735, 0.57735] // vector_axis([[10,0,10], [0,0,0], [-10,10,0]]); // Returns: [-0.57735, -0.57735, 0.57735] function vector_axis(v1,v2=undef,v3=undef) = - (is_list(v1) && is_list(v1[0]) && is_undef(v2) && is_undef(v3))? ( - assert(is_vector(v1.x)) - assert(is_vector(v1.y)) - len(v1)==3? assert(is_vector(v1.z)) vector_axis(v1.x, v1.y, v1.z) : - len(v1)==2? vector_axis(v1.x, v1.y) : - assert(false, "Bad arguments.") - ) : - (is_vector(v1) && is_vector(v2) && is_vector(v3))? vector_axis(v1-v2, v3-v2) : - (is_vector(v1) && is_vector(v2) && is_undef(v3))? let( + is_vector(v3) + ? assert(is_consistent([v3,v2,v1]), "Bad arguments.") + vector_axis(v1-v2, v3-v2) + : + assert( is_undef(v3), "Bad arguments.") + is_undef(v2) + ? assert( is_list(v1), "Bad arguments.") + len(v1) == 2 + ? vector_axis(v1[0],v1[1]) + : vector_axis(v1[0],v1[1],v1[2]) + : + assert( + is_vector(v1,zero=false) && + is_vector(v2,zero=false) && + is_consistent([v1,v2]), + "Bad arguments." + ) + let( eps = 1e-6, - v1 = point3d(v1/norm(v1)), - v2 = point3d(v2/norm(v2)), - v3 = (norm(v1-v2) > eps && norm(v1+v2) > eps)? v2 : - (norm(vabs(v2)-UP) > eps)? UP : - RIGHT - ) unit(cross(v1,v3)) : assert(false, "Bad arguments."); + w1 = point3d(v1/norm(v1)), + w2 = point3d(v2/norm(v2)), + w3 = (norm(w1-w2) > eps && norm(w1+w2) > eps) ? w2 + : (norm(vabs(w2)-UP) > eps) ? UP + : RIGHT + ) unit(cross(w1,w3)); // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap diff --git a/version.scad b/version.scad index 4402f4c7..406f9e4b 100644 --- a/version.scad +++ b/version.scad @@ -8,7 +8,7 @@ ////////////////////////////////////////////////////////////////////// -BOSL_VERSION = [2,0,389]; +BOSL_VERSION = [2,0,397]; // Section: BOSL Library Version Functions