diff --git a/attachments.scad b/attachments.scad
index f36cce5..fc8be8e 100644
--- a/attachments.scad
+++ b/attachments.scad
@@ -326,7 +326,7 @@ function attach_transform(anchor=CENTER, spin=0, orient=UP, geom, p) =
 					ang = vector_angle(anch[2], DOWN),
 					axis = vector_axis(anch[2], DOWN),
 					ang2 = (anch[2]==UP || anch[2]==DOWN)? 0 : 180-anch[3],
-					axis2 = rotate_points3d([axis],[0,0,ang2])[0]
+					axis2 = rot(p=axis,[0,0,ang2])
 				)
 				affine3d_rot_by_axis(axis2,ang) *
 				affine3d_zrot(ang2+spin) *
diff --git a/beziers.scad b/beziers.scad
index 5ab58ad..b4bcc4a 100644
--- a/beziers.scad
+++ b/beziers.scad
@@ -884,12 +884,14 @@ function _bezier_triangle(tri, splinesteps=16, vnf=EMPTY_VNF) =
 //   trace_bezier_patches([patch], size=1, showcps=true);
 function bezier_patch_flat(size=[100,100], N=4, spin=0, orient=UP, trans=[0,0,0]) =
 	let(
-		patch = [for (x=[0:1:N]) [for (y=[0:1:N]) vmul(point3d(size),[x/N-0.5, 0.5-y/N, 0])]]
-	) [for (row=patch)
-		translate_points(v=trans,
-			rotate_points3d(a=spin, from=UP, to=orient, row)
-		)
-	];
+		patch = [
+			for (x=[0:1:N]) [
+				for (y=[0:1:N])
+				vmul(point3d(size), [x/N-0.5, 0.5-y/N, 0])
+			]
+		],
+		m = move(trans) * rot(a=spin, from=UP, to=orient)
+	) [for (row=patch) apply(m, row)];
 
 
 
diff --git a/coords.scad b/coords.scad
index 2d3edd0..f1d707f 100644
--- a/coords.scad
+++ b/coords.scad
@@ -100,121 +100,6 @@ function path4d(points, fill=0) =
     result + repeat(addition, len(result));
 
 
-// Function: translate_points()
-// Usage:
-//   translate_points(pts, v);
-// Description:
-//   Moves each point in an array by a given amount.
-// Arguments:
-//   pts = List of points to translate.
-//   v = Amount to translate points by.
-function translate_points(pts, v=[0,0,0]) =
-	pts==[]? [] : let(
-		v=point3d(v)
-	) [for (pt = pts) pt+v];
-
-
-// Function: scale_points()
-// Usage:
-//   scale_points(pts, v, [cp]);
-// Description:
-//   Scales each point in an array by a given amount, around a given centerpoint.
-// Arguments:
-//   pts = List of points to scale.
-//   v = A vector with a scaling factor for each axis.
-//   cp = Centerpoint to scale around.
-function scale_points(pts, v=[1,1,1], cp=[0,0,0]) =
-	pts==[]? [] : let(
-		cp = point3d(cp),
-		v = point3d(v,fill=1)
-	) [for (pt = pts) vmul(pt-cp,v)+cp];
-
-
-// Function: rotate_points2d()
-// Usage:
-//   rotate_points2d(pts, a, [cp]);
-// Description:
-//   Rotates each 2D point in an array by a given amount, around an optional centerpoint.
-// Arguments:
-//   pts = List of 3D points to rotate.
-//   a = Angle to rotate by.
-//   cp = 2D Centerpoint to rotate around.  Default: `[0,0]`
-function rotate_points2d(pts, a, cp=[0,0]) =
-	approx(a,0)? pts :
-	let(
-		cp = point2d(cp),
-		pts = path2d(pts),
-		m = affine2d_zrot(a)
-	) [for (pt = pts) point2d(m*concat(pt-cp, [1])+cp)];
-
-
-// Function: rotate_points3d()
-// Usage:
-//   rotate_points3d(pts, a, [cp], [reverse]);
-//   rotate_points3d(pts, a, v, [cp], [reverse]);
-//   rotate_points3d(pts, from, to, [a], [cp], [reverse]);
-// Description:
-//   Rotates each 3D point in an array by a given amount, around a given centerpoint.
-// Arguments:
-//   pts = List of points to rotate.
-//   a = Rotation angle(s) in degrees.
-//   v = If given, axis vector to rotate around.
-//   cp = Centerpoint to rotate around.
-//   from = If given, the vector to rotate something from.  Used with `to`.
-//   to = If given, the vector to rotate something to.  Used with `from`.
-//   reverse = If true, performs an exactly reversed rotation.
-function rotate_points3d(pts, a=0, v=undef, cp=[0,0,0], from=undef, to=undef, reverse=false) =
-	assert(is_undef(from)==is_undef(to), "`from` and `to` must be given together.")
-	(is_undef(from) && (a==0 || a==[0,0,0]))? pts :
-	let (
-		from = is_undef(from)? undef : (from / norm(from)),
-		to = is_undef(to)? undef : (to / norm(to)),
-		cp = point3d(cp),
-		pts2 = path3d(pts)
-	)
-	(!is_undef(from) && approx(from,to) && (a==0 || a == [0,0,0]))? pts2 :
-	let (
-		mrot = reverse? (
-			!is_undef(from)? (
-				assert(norm(from)>0, "The from argument cannot equal [0,0] or [0,0,0]")
-				assert(norm(to)>0, "The to argument cannot equal [0,0] or [0,0,0]")
-				let (
-					ang = vector_angle(from, to),
-					v = vector_axis(from, to)
-				)
-				affine3d_rot_by_axis(from, -a) * affine3d_rot_by_axis(v, -ang)
-			) : !is_undef(v)? (
-				affine3d_rot_by_axis(v, -a)
-			) : is_num(a)? (
-				affine3d_zrot(-a)
-			) : (
-				affine3d_xrot(-a.x) * affine3d_yrot(-a.y) * affine3d_zrot(-a.z)
-			)
-		) : (
-			!is_undef(from)? (
-				assert(norm(from)>0, "The from argument cannot equal [0,0] or [0,0,0]")
-				assert(norm(to)>0, "The to argument cannot equal [0,0] or [0,0,0]")
-				let (
-					from = from / norm(from),
-					to = to / norm(from),
-					ang = vector_angle(from, to),
-					v = vector_axis(from, to)
-				)
-				affine3d_rot_by_axis(v, ang) * affine3d_rot_by_axis(from, a)
-			) : !is_undef(v)? (
-				affine3d_rot_by_axis(v, a)
-			) : is_num(a)? (
-				affine3d_zrot(a)
-			) : (
-				affine3d_zrot(a.z) * affine3d_yrot(a.y) * affine3d_xrot(a.x)
-			)
-		),
-		m = affine3d_translate(cp) * mrot * affine3d_translate(-cp)
-	)
-	[for (pt = pts2) point3d(m*concat(pt, fill=1))];
-
-
-
 // Section: Coordinate Systems
 
 // Function: polar_to_xy()
@@ -289,7 +174,7 @@ function project_plane(point, a, b, c) =
 		v = unit(c-a),
 		n = unit(cross(u,v)),
 		w = unit(cross(n,u)),
-		relpoint = is_vector(point)? (point-a) : translate_points(point,-a)
+		relpoint = is_vector(point)? (point-a) : move(-a,p=point)
 	) relpoint * transpose([w,u]);
 
 
@@ -320,7 +205,7 @@ function lift_plane(point, a, b, c) =
 		n = unit(cross(u,v)),
 		w = unit(cross(n,u)),
 		remapped = point*[w,u]
-	) is_vector(remapped)? (a+remapped) : translate_points(remapped,a);
+	) is_vector(remapped)? (a+remapped) : move(a,p=remapped);
 
 
 // Function: cylindrical_to_xyz()
diff --git a/examples/orientations.scad b/examples/orientations.scad
index 9cabdfd..2a55d6d 100644
--- a/examples/orientations.scad
+++ b/examples/orientations.scad
@@ -65,7 +65,7 @@ module orient_cubes() {
 	}
 
 	for (ang = [0:90:270]) {
-		off = rotate_points3d([40*BACK],ang)[0];
+		off = rot(p=40*BACK,ang);
 		translate(off) {
 			orient_cube(ang);
 		}
diff --git a/examples/randomized_fractal_tree.scad b/examples/randomized_fractal_tree.scad
index 0cb34c2..c3d2592 100644
--- a/examples/randomized_fractal_tree.scad
+++ b/examples/randomized_fractal_tree.scad
@@ -10,7 +10,7 @@ module leaf(s) {
 	];
 	xrot(90)
 	linear_sweep_bezier(
-		scale_points(path, [s,s]/2),
+		path * s/2,
 		height=0.02
 	);
 }
diff --git a/polyhedra.scad b/polyhedra.scad
index 9e10b50..2afc770 100644
--- a/polyhedra.scad
+++ b/polyhedra.scad
@@ -320,14 +320,14 @@ module regular_polyhedron(
 	in_radius = entry[5];
 	if (draw){
 		if (rounding==0)
-			polyhedron(translate_points(scaled_points, translation), faces = face_triangles);
+			polyhedron(move(p=scaled_points, translation), faces = face_triangles);
 		else {
 			fn = segs(rounding);
 			rounding = rounding/cos(180/fn);
 			adjusted_scale = 1 - rounding / in_radius;
 			minkowski(){
 				sphere(r=rounding, $fn=fn);
-				polyhedron(translate_points(adjusted_scale*scaled_points,translation), faces = face_triangles);
+				polyhedron(move(p=adjusted_scale*scaled_points,translation), faces = face_triangles);
 			}
 		}
 	}
@@ -335,13 +335,13 @@ module regular_polyhedron(
 		maxrange = repeat ? len(faces)-1 : $children-1;
 		for(i=[0:1:maxrange]) {
 			// Would like to orient so an edge (longest edge?) is parallel to x axis
-			facepts = translate_points(select(scaled_points, faces[i]), translation);
+			facepts = move(p=select(scaled_points, faces[i]), translation);
 			center = mean(facepts);
-			rotatedface = rotate_points3d(translate_points(facepts,-center), from=face_normals[i], to=[0,0,1]);
+			rotatedface = rot(p=move(p=facepts,-center), from=face_normals[i], to=[0,0,1]);
 			clockwise = sortidx([for(pt=rotatedface) -atan2(pt.y,pt.x)]);
 			$face = rotate_children?
 						path2d(select(rotatedface,clockwise)) :
-						select(translate_points(facepts,-center), clockwise);
+						select(move(p=facepts,-center), clockwise);
 			$faceindex = i;
 			$center = -translation-center;
 			translate(center)
@@ -681,15 +681,15 @@ function regular_polyhedron_info(
 		facedown = facedown == true ? (stellate==false? entry[facevertices][0] : 3) : facedown,
 		down_direction = facedown == false?  [0,0,-1] :
 			faces_normals_vertices[1][search(facedown, faces_vertex_count)[0]],
-		scaled_points = scalefactor * rotate_points3d(faces_normals_vertices[2], from=down_direction, to=[0,0,-1]),
+		scaled_points = scalefactor * rot(p=faces_normals_vertices[2], from=down_direction, to=[0,0,-1]),
 		bounds = pointlist_bounds(scaled_points),
 		boundtable = [bounds[0], [0,0,0], bounds[1]],
 		translation = [for(i=[0:2]) -boundtable[1+anchor[i]][i]],
-		face_normals = rotate_points3d(faces_normals_vertices[1], from=down_direction, to=[0,0,-1]),
+		face_normals = rot(p=faces_normals_vertices[1], from=down_direction, to=[0,0,-1]),
 		side_length = scalefactor * entry[edgelen]
 	)
 	info == "fullentry" ? [scaled_points, translation,stellate ? faces : face_triangles, faces, face_normals, side_length*entry[in_radius]] :
-	info == "vertices" ? translate_points(scaled_points,translation) :
+	info == "vertices" ? move(p=scaled_points,translation) :
 	info == "faces" ? faces :
 	info == "face normals" ? face_normals :
 	info == "in_radius" ? side_length * entry[in_radius] :
diff --git a/rounding.scad b/rounding.scad
index 68c5708..08a95b4 100644
--- a/rounding.scad
+++ b/rounding.scad
@@ -763,7 +763,7 @@ module offset_sweep(
 	top_start_ind = len(vertices_faces_bot[0]);
 	initial_vertices_top = zip(path, repeat(middle,len(path)));
 	vertices_faces_top = make_polyhedron(
-		path, translate_points(offsets_top,[0,middle]),
+		path, move(p=offsets_top,[0,middle]),
 		struct_val(top,"offset"), !clockwise,
 		struct_val(top,"quality"),
 		struct_val(top,"check_valid"),
@@ -1278,7 +1278,7 @@ function _stroke_end(width,left, right, spec) =
 		angle = struct_val(spec,"absolute")?
 			angle_between_lines(left[0]-right[0],[cos(user_angle),sin(user_angle)]) :
 			user_angle,
-		endseg = [center, rotate_points2d([left[0]],angle, cp=center)[0]],
+		endseg = [center, rot(p=[left[0]], angle, cp=center)[0]],
 		intright = angle>0,
 		pathclip = _path_line_intersection(intright? right : left, endseg),
 		pathextend = line_intersection(endseg, select(intright? left:right,0,1))
diff --git a/tests/test_coords.scad b/tests/test_coords.scad
index 2dc5f9d..0b321d5 100644
--- a/tests/test_coords.scad
+++ b/tests/test_coords.scad
@@ -54,56 +54,6 @@ module test_path4d() {
 test_path4d();
 
 
-module test_translate_points() {
-	pts = [[0,0,1], [0,1,0], [1,0,0], [0,0,-1], [0,-1,0], [-1,0,0]];
-	assert(translate_points(pts, v=[1,2,3]) == [[1,2,4], [1,3,3], [2,2,3], [1,2,2], [1,1,3], [0,2,3]]);
-	assert(translate_points(pts, v=[-1,-2,-3]) == [[-1,-2,-2], [-1,-1,-3], [0,-2,-3], [-1,-2,-4], [-1,-3,-3], [-2,-2,-3]]);
-	assert(translate_points(pts, v=[1,2]) == [[1,2,1], [1,3,0], [2,2,0], [1,2,-1], [1,1,0], [0,2,0]]);
-	pts2 = [[0,1], [1,0], [0,-1], [-1,0]];
-	assert(translate_points(pts2, v=[1,2]) == [[1,3], [2,2], [1,1], [0,2]]);
-}
-test_translate_points();
-
-
-module test_scale_points() {
-	pts = [[0,0,1], [0,1,0], [1,0,0], [0,0,-1], [0,-1,0], [-1,0,0]];
-	assert(scale_points(pts, v=[2,3,4]) == [[0,0,4], [0,3,0], [2,0,0], [0,0,-4], [0,-3,0], [-2,0,0]]);
-	assert(scale_points(pts, v=[-2,-3,-4]) == [[0,0,-4], [0,-3,0], [-2,0,0], [0,0,4], [0,3,0], [2,0,0]]);
-	assert(scale_points(pts, v=[1,1,1]) == [[0,0,1], [0,1,0], [1,0,0], [0,0,-1], [0,-1,0], [-1,0,0]]);
-	assert(scale_points(pts, v=[-1,-1,-1]) == [[0,0,-1], [0,-1,0], [-1,0,0], [0,0,1], [0,1,0], [1,0,0]]);
-	pts2 = [[0,1], [1,0], [0,-1], [-1,0]];
-	assert(scale_points(pts2, v=[2,3]) == [[0,3], [2,0], [0,-3], [-2,0]]);
-}
-test_scale_points();
-
-
-module test_rotate_points2d() {
-	pts = [[0,1], [1,0], [0,-1], [-1,0]];
-	s = sin(45);
-	assert(rotate_points2d(pts,45) == [[-s,s],[s,s],[s,-s],[-s,-s]]);
-	assert(rotate_points2d(pts,90) == [[-1,0],[0,1],[1,0],[0,-1]]);
-	assert(rotate_points2d(pts,90,cp=[1,0]) == [[0,-1],[1,0],[2,-1],[1,-2]]);
-}
-test_rotate_points2d();
-
-
-module test_rotate_points3d() {
-	pts = [[0,0,1], [0,1,0], [1,0,0], [0,0,-1], [0,-1,0], [-1,0,0]];
-	assert(rotate_points3d(pts, [90,0,0]) == [[0,-1,0], [0,0,1], [1,0,0], [0,1,0], [0,0,-1], [-1,0,0]]);
-	assert(rotate_points3d(pts, [0,90,0]) == [[1,0,0], [0,1,0], [0,0,-1], [-1,0,0], [0,-1,0], [0,0,1]]);
-	assert(rotate_points3d(pts, [0,0,90]) == [[0,0,1], [-1,0,0], [0,1,0], [0,0,-1], [1,0,0], [0,-1,0]]);
-	assert(rotate_points3d(pts, [0,0,90],cp=[2,0,0]) == [[2,-2,1], [1,-2,0], [2,-1,0], [2,-2,-1], [3,-2,0], [2,-3,0]]);
-	assert(rotate_points3d(pts, 90, v=UP) == [[0,0,1], [-1,0,0], [0,1,0], [0,0,-1], [1,0,0], [0,-1,0]]);
-	assert(rotate_points3d(pts, 90, v=DOWN) == [[0,0,1], [1,0,0], [0,-1,0], [0,0,-1], [-1,0,0], [0,1,0]]);
-	assert(rotate_points3d(pts, 90, v=RIGHT)  == [[0,-1,0], [0,0,1], [1,0,0], [0,1,0], [0,0,-1], [-1,0,0]]);
-	assert(rotate_points3d(pts, from=UP, to=BACK) == [[0,1,0], [0,0,-1], [1,0,0], [0,-1,0], [0,0,1], [-1,0,0]]);
-	assert(rotate_points3d(pts, 90, from=UP, to=BACK), [[0,1,0], [-1,0,0], [0,0,-1], [0,-1,0], [1,0,0], [0,0,1]]);
-	assert(rotate_points3d(pts, from=UP, to=UP*2) == [[0,0,1], [0,1,0], [1,0,0], [0,0,-1], [0,-1,0], [-1,0,0]]);
-	assert(rotate_points3d(pts, from=UP, to=DOWN*2) == [[0,0,-1], [0,1,0], [-1,0,0], [0,0,1], [0,-1,0], [1,0,0]]);
-}
-test_rotate_points3d();
-
-
 module test_polar_to_xy() {
 	assert(approx(polar_to_xy(20,45), [20/sqrt(2), 20/sqrt(2)]));
 	assert(approx(polar_to_xy(20,135), [-20/sqrt(2), 20/sqrt(2)]));
diff --git a/tests/test_quaternions.scad b/tests/test_quaternions.scad
index 4331ff0..4aeb08c 100644
--- a/tests/test_quaternions.scad
+++ b/tests/test_quaternions.scad
@@ -11,11 +11,11 @@ function rec_cmp(a,b,eps=1e-9) =
 
 module verify_f(actual,expected) {
 	if (!rec_cmp(actual,expected)) {
-		echo(str("Expected: ",fmtf(expected,10)));
+		echo(str("Expected: ",fmt_float(expected,10)));
 		echo(str("        : ",expected));
-		echo(str("Actual  : ",fmtf(actual,10)));
+		echo(str("Actual  : ",fmt_float(actual,10)));
 		echo(str("        : ",actual));
-		echo(str("Delta   : ",fmtf(expected-actual,10)));
+		echo(str("Delta   : ",fmt_float(expected-actual,10)));
 		echo(str("        : ",expected-actual));
 		assert(approx(expected,actual));
 	}
diff --git a/tests/test_transforms.scad b/tests/test_transforms.scad
new file mode 100644
index 0000000..1db7ff5
--- /dev/null
+++ b/tests/test_transforms.scad
@@ -0,0 +1,159 @@
+include <BOSL2/std.scad>
+
+module test(got,expect,extra_info) {
+	if (
+		is_undef(expect) != is_undef(got) ||
+		expect*0 != got*0 ||
+		(is_vnf(expect) && !all([for (i=idx(expect[0])) approx(got[0][i],expect[0][i])]) && got[1]!=expect[1]) ||
+		(is_matrix(expect) && !all([for (i=idx(expect)) approx(got[i],expect[i])])) ||
+		(got!=expect && !approx(got, expect))
+	) {
+		fmt = is_int(expect)? "{:.14i}" :
+			is_num(expect)? "{:.14g}" :
+			is_vector(expect)? "{:.14g}" :
+			"{}";
+		echofmt(str("Expected: ",fmt),[expect]);
+		echofmt(str("But Got : ",fmt),[got]);
+		if (expect*0 == got*0) {
+			echofmt(str("Delta is: ",fmt),[expect-got]);
+		}
+		if (!is_undef(extra_info)) {
+			echo(str("Extra Info: ",extra_info));
+		}
+		assert(false, "TEST FAILED!");
+	}
+}
+
+
+module test_rot() {
+	pts2d = 50 * [for (x=[-1,0,1],y=[-1,0,1]) [x,y]];
+	pts3d = 50 * [for (x=[-1,0,1],y=[-1,0,1],z=[-1,0,1]) [x,y,z]];
+	vecs2d = [
+		for (x=[-1,0,1], y=[-1,0,1]) if(x!=0||y!=0) [x,y],
+		polar_to_xy(1, -75),
+		polar_to_xy(1,  75)
+	];
+	vecs3d = [
+		LEFT, RIGHT, FRONT, BACK, DOWN, UP,
+		spherical_to_xyz(1, -30,  45),
+		spherical_to_xyz(1,   0,  45),
+		spherical_to_xyz(1,  30,  45),
+		spherical_to_xyz(2,  30,  45),
+		spherical_to_xyz(1, -30, 135),
+		spherical_to_xyz(2, -30, 135),
+		spherical_to_xyz(1,   0, 135),
+		spherical_to_xyz(1,  30, 135),
+		spherical_to_xyz(1, -30,  75),
+		spherical_to_xyz(1,  45,  45),
+	];
+	angs = [-180, -90, -45, 0, 30, 45, 90];
+	for (a = [-360*3:360:360*3]) {
+		test(rot(a), affine3d_identity(), extra_info=str("rot(",a,") != identity"));
+		test(rot(a,p=pts2d), pts2d, extra_info=str("rot(",a,",p=...), 2D"));
+		test(rot(a,p=pts3d), pts3d, extra_info=str("rot(",a,",p=...), 3D"));
+	}
+	test(rot(90), [[0,-1,0,0],[1,0,0,0],[0,0,1,0],[0,0,0,1]])
+	for (a=angs) {
+		test(rot(a), affine3d_zrot(a), extra_info=str("Z angle (only) = ",a));
+		test(rot([a,0,0]), affine3d_xrot(a), extra_info=str("X angle = ",a));
+		test(rot([0,a,0]), affine3d_yrot(a), extra_info=str("Y angle = ",a));
+		test(rot([0,0,a]), affine3d_zrot(a), extra_info=str("Z angle = ",a));
+
+		test(rot(a,p=pts2d), apply(affine3d_zrot(a),pts2d), extra_info=str("Z angle (only) = ",a, ", p=..., 2D"));
+		test(rot([0,0,a],p=pts2d), apply(affine3d_zrot(a),pts2d), extra_info=str("Z angle = ",a, ", p=..., 2D"));
+
+		test(rot(a,p=pts3d), apply(affine3d_zrot(a),pts3d), extra_info=str("Z angle (only) = ",a, ", p=..., 3D"));
+		test(rot([a,0,0],p=pts3d), apply(affine3d_xrot(a),pts3d), extra_info=str("X angle = ",a, ", p=..., 3D"));
+		test(rot([0,a,0],p=pts3d), apply(affine3d_yrot(a),pts3d), extra_info=str("Y angle = ",a, ", p=..., 3D"));
+		test(rot([0,0,a],p=pts3d), apply(affine3d_zrot(a),pts3d), extra_info=str("Z angle = ",a, ", p=..., 3D"));
+	}
+	for (xa=angs, ya=angs, za=angs) {
+		test(
+			rot([xa,ya,za]),
+			affine3d_chain([
+				affine3d_xrot(xa),
+				affine3d_yrot(ya),
+				affine3d_zrot(za)
+			]),
+			extra_info=str("[X,Y,Z] = ",[xa,ya,za])
+		);
+		test(
+			rot([xa,ya,za],p=pts3d),
+			apply(
+				affine3d_chain([
+					affine3d_xrot(xa),
+					affine3d_yrot(ya),
+					affine3d_zrot(za)
+				]),
+				pts3d
+			),
+			extra_info=str("[X,Y,Z] = ",[xa,ya,za], ", p=...")
+		);
+	}
+	for (vec1 = vecs3d) {
+		for (ang = angs) {
+			test(
+				rot(a=ang, v=vec1),
+				affine3d_rot_by_axis(vec1,ang),
+				extra_info=str("a = ",ang,", v = ", vec1)
+			);
+			test(
+				rot(a=ang, v=vec1, p=pts3d),
+				apply(affine3d_rot_by_axis(vec1,ang), pts3d),
+				extra_info=str("a = ",ang,", v = ", vec1, ", p=...")
+			);
+		}
+	}
+	for (vec1 = vecs2d) {
+		for (vec2 = vecs2d) {
+			test(
+				rot(from=vec1, to=vec2, p=pts2d, planar=true),
+				apply(affine2d_zrot(vang(vec2)-vang(vec1)), pts2d),
+				extra_info=str(
+					"from = ", vec1, ", ",
+					"to = ", vec2, ", ",
+					"planar = ", true, ", ",
+					"p=..., 2D"
+				)
+			);
+		}
+	}
+	for (vec1 = vecs3d) {
+		for (vec2 = vecs3d) {
+			for (a = angs) {
+				test(
+					rot(from=vec1, to=vec2, a=a),
+					affine3d_chain([
+						affine3d_zrot(a),
+						affine3d_rot_from_to(vec1,vec2)
+					]),
+					extra_info=str(
+						"from = ", vec1, ", ",
+						"to = ", vec2, ", ",
+						"a = ", a
+					)
+				);
+				test(
+					rot(from=vec1, to=vec2, a=a, p=pts3d),
+					apply(
+						affine3d_chain([
+							affine3d_zrot(a),
+							affine3d_rot_from_to(vec1,vec2)
+						]),
+						pts3d
+					),
+					extra_info=str(
+						"from = ", vec1, ", ",
+						"to = ", vec2, ", ",
+						"a = ", a, ", ",
+						"p=..., 3D"
+					)
+				);
+			}
+		}
+	}
+}
+test_rot();
+
+
+// vim: noexpandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap
diff --git a/threading.scad b/threading.scad
index 1d5a1bf..0698a5b 100644
--- a/threading.scad
+++ b/threading.scad
@@ -60,7 +60,7 @@ module thread_helix(base_d, pitch, thread_depth=undef, thread_angle=15, twist=72
 			[0, -cap/2-dz],
 		]
 	);
-	pline = scale_points(profile, [1,1,1]*pitch);
+	pline = profile * pitch;
 	dir = left_handed? -1 : 1;
 	idir = internal? -1 : 1;
 	attachable(anchor,spin,orient, r=r, l=h) {
diff --git a/transforms.scad b/transforms.scad
index abe7269..93a0f77 100644
--- a/transforms.scad
+++ b/transforms.scad
@@ -337,102 +337,56 @@ function up(z=0,p=undef) = move([0,0,z],p=p);
 //   stroke(rot(30,p=path), closed=true);
 module rot(a=0, v=undef, cp=undef, from=undef, to=undef, reverse=false)
 {
-	if (!is_undef(cp)) {
-		translate(cp) rot(a=a, v=v, from=from, to=to, reverse=reverse) translate(-cp) children();
-	} else if (!is_undef(from)) {
-		assert(!is_undef(to), "`from` and `to` should be used together.");
-		from = point3d(from);
-		to = point3d(to);
-		axis = vector_axis(from, to);
-		ang = vector_angle(from, to);
-		if (ang < 0.0001 && a == 0) {
-			children();  // May be slightly faster?
-		} else if (reverse) {
-			rotate(a=-ang, v=axis) rotate(a=-a, v=from) children();
-		} else {
-			rotate(a=ang, v=axis) rotate(a=a, v=from) children();
-		}
-	} else if (a == 0) {
-		children();  // May be slightly faster?
-	} else if (reverse) {
-		if (!is_undef(v)) {
-			rotate(a=-a, v=v) children();
-		} else if (is_num(a)) {
-			rotate(-a) children();
-		} else {
-			rotate([-a[0],0,0]) rotate([0,-a[1],0]) rotate([0,0,-a[2]]) children();
-		}
-	} else {
-		rotate(a=a, v=v) children();
-	}
+	m = rot(a=a, v=v, cp=cp, from=from, to=to, reverse=reverse, planar=false);
+	multmatrix(m) children();
 }
 
-function rot(a=0, v=undef, cp=undef, from=undef, to=undef, reverse=false, p=undef, planar=false) =
+function rot(a=0, v, cp, from, to, reverse=false, planar=false, p, _m) =
 	assert(is_undef(from)==is_undef(to), "from and to must be specified together.")
-	let(
-		rev = reverse? -1 : 1,
-		from = is_undef(from)? undef : point3d(from),
-		to = is_undef(to)? undef : point3d(to)
-	)
 	is_undef(p)? (
-		is_undef(cp)? (
-			planar? (
-				is_undef(from)? affine2d_zrot(a*rev) :
-				affine2d_zrot(vector_angle(from,to)*sign(vector_axis(from,to)[2])*rev)
-			) : (
-				!is_undef(from)? affine3d_chain([
-					affine3d_zrot(a*rev),
-					affine3d_rot_by_axis(
-						vector_axis(from,to),
-						vector_angle(from,to)*rev
-					)
-				]) :
-				!is_undef(v)? affine3d_rot_by_axis(v,a*rev) :
-				is_num(a)? affine3d_zrot(a*rev) :
-				reverse? affine3d_chain([affine3d_zrot(-a.z),affine3d_yrot(-a.y),affine3d_xrot(-a.x)]) :
-				affine3d_chain([affine3d_xrot(a.x),affine3d_yrot(a.y),affine3d_zrot(a.z)])
-			)
-		) : (
-			planar? (
-				affine2d_chain([
-					move(-cp),
-					rot(a=a, v=v, from=from, to=to, reverse=reverse, planar=true),
-					move(cp)
-				])
-			) : (
-				affine3d_chain([
-					move(-cp),
-					rot(a=a, v=v, from=from, to=to, reverse=reverse),
-					move(cp)
-				])
-			)
-		)
+		planar? let(
+			cp = is_undef(cp)? cp : point2d(cp),
+			m1 = is_undef(from)? affine2d_zrot(a) :
+				assert(is_vector(from))
+				assert(!approx(norm(from),0))
+				assert(approx(point3d(from).z, 0))
+				assert(is_vector(to))
+				assert(!approx(norm(to),0))
+				assert(approx(point3d(to).z, 0))
+				affine2d_zrot(
+					vang(point2d(to)) -
+					vang(point2d(from))
+				),
+			m2 = is_undef(cp)? m1 : (move(cp) * m1 * move(-cp)),
+			m3 = reverse? matrix_inverse(m2) : m2
+		) m3 : let(
+			from = is_undef(from)? undef : point3d(from),
+			to = is_undef(to)? undef : point3d(to),
+			cp = is_undef(cp)? undef : point3d(cp),
+			m1 = !is_undef(from)? (
+					assert(is_vector(from))
+					assert(!approx(norm(from),0))
+					assert(is_vector(to))
+					assert(!approx(norm(to),0))
+					affine3d_rot_from_to(from,to) * affine3d_zrot(a)
+				) :
+				!is_undef(v)? affine3d_rot_by_axis(v,a) :
+				is_num(a)? affine3d_zrot(a) :
+				affine3d_zrot(a.z) * affine3d_yrot(a.y) * affine3d_xrot(a.x),
+			m2 = is_undef(cp)? m1 : (move(cp) * m1 * move(-cp)),
+			m3 = reverse? matrix_inverse(m2) : m2
+		) m3
 	) : (
 		assert(is_list(p))
-		is_num(p.x)? (
-			rot(a=a, v=v, cp=cp, from=from, to=to, reverse=reverse, p=[p], planar=planar)[0]
-		) : is_vnf(p)? (
-			[rot(a=a, v=v, cp=cp, from=from, to=to, reverse=reverse, p=p.x, planar=planar), p.y]
-		) : is_list(p.x) && is_list(p.x.x)? (
-			[for (l=p) rot(a=a, v=v, cp=cp, from=from, to=to, reverse=reverse, p=l, planar=planar)]
-		) : (
-			(
-				(planar || (p!=[] && len(p[0])==2)) && !(
-					(is_vector(a) && norm(point2d(a))>0) ||
-					(!is_undef(v) && norm(point2d(v))>0 && !approx(a,0)) ||
-					(!is_undef(from) && !approx(from,to) && !(abs(point3d(from).z)>0 || abs(point3d(to).z))) ||
-					(!is_undef(from) && approx(from,to) && norm(point2d(from))>0 && a!=0)
-				)
-			)? (
-				is_undef(from)? rotate_points2d(p, a=a*rev, cp=cp) : (
-					approx(from,to)&&approx(a,0)? p :
-					rotate_points2d(p, a=vector_angle(from,to)*sign(vector_axis(from,to)[2])*rev, cp=cp)
-				)
-			) : (
-				let( cp = is_undef(cp)? [0,0,0] : cp )
-				rotate_points3d(p, a=a, v=v, cp=cp, from=from, to=to, reverse=reverse)
-			)
-		)
+		let(
+			m = !is_undef(_m)? _m :
+				rot(a=a, v=v, cp=cp, from=from, to=to, reverse=reverse, planar=planar),
+			res = p==[]? [] :
+				is_vector(p)? apply(m, p) :
+				is_vnf(p)? [apply(m, p[0]), p[1]] :
+				is_list(p[0])? [for (pp=p) rot(p=pp, _m=m)] :
+				assert(false, "The p argument for rot() is not a point, path, patch, matrix, or VNF.")
+		) res
 	);
 
 
diff --git a/version.scad b/version.scad
index 8b47e5c..253fbfb 100644
--- a/version.scad
+++ b/version.scad
@@ -8,7 +8,7 @@
 //////////////////////////////////////////////////////////////////////
 
 
-BOSL_VERSION = [2,0,215];
+BOSL_VERSION = [2,0,216];
 
 
 // Section: BOSL Library Version Functions