diff --git a/shapes.scad b/shapes.scad
index 31e35b1a..40a950b1 100644
--- a/shapes.scad
+++ b/shapes.scad
@@ -1689,112 +1689,118 @@ module arced_slot(
 }
 
 
-// Module: heightfield()
-// Usage:
-//   heightfield(heightfield, [size], [bottom]);
+// Function&Module: heightfield()
+// Usage: As Module
+//   heightfield(data, <size>, <bottom>, <maxz>, <xrange>, <yrange>, <style>, <convexity>);
+// Usage: As Function
+//   vnf = heightfield(data, <size>, <bottom>, <maxz>, <xrange>, <yrange>, <style>);
 // Description:
-//   Given a regular rectangular 2D grid of scalar values, generates a 3D surface where the height at
-//   any given point is the scalar value for that position.
+//   Given a regular rectangular 2D grid of scalar values, or a function literal, generates a 3D
+//   surface where the height at any given point is the scalar value for that position.
 // Arguments:
-//   heightfield = The 2D rectangular array of heights.
-//   size = The [X,Y] size of the surface to create.  If given as a scalar, use it for both X and Y sizes.
-//   bottom = The Z coordinate for the bottom of the heightfield object to create.  Must be less than the minimum heightfield value.  Default: 0
-//   convexity = Max number of times a line could intersect a wall of the surface being formed.
+//   data = This is either the 2D rectangular array of heights, or a function literal that takes X and Y arguments.
+//   size = The [X,Y] size of the surface to create.  If given as a scalar, use it for both X and Y sizes. Default: `[100,100]`
+//   bottom = The Z coordinate for the bottom of the heightfield object to create.  Any heights lower than this will be truncated to very slightly above this height.  Default: -20
+//   maxz = The maximum height to model.  Truncates anything taller to this height.  Default: 99
+//   xrange = A range of values to iterate X over when calculating a surface from a function literal.  Default: [-1 : 0.01 : 1]
+//   yrange = A range of values to iterate Y over when calculating a surface from a function literal.  Default: [-1 : 0.01 : 1]
+//   style = The style of subdividing the quads into faces.  Valid options are "default", "alt", and "quincunx".  Default: "default"
+//   convexity = Max number of times a line could intersect a wall of the surface being formed. Module only.  Default: 10
+//   anchor = Translate so anchor point is at origin (0,0,0).  See [anchor](attachments.scad#anchor).  Default: `CENTER`
+//   spin = Rotate this many degrees around the Z axis.  See [spin](attachments.scad#spin).  Default: `0`
+//   orient = Vector to rotate top towards.  See [orient](attachments.scad#orient).  Default: `UP`
 // Example:
-//   heightfield(size=[100,100], bottom=-20, heightfield=[
-//       for (x=[-180:4:180]) [for(y=[-180:4:180]) 10*cos(3*norm([x,y]))]
+//   heightfield(size=[100,100], bottom=-20, data=[
+//       for (y=[-180:4:180]) [for(x=[-180:4:180]) 10*cos(3*norm([x,y]))]
 //   ]);
 // Example:
 //   intersection() {
-//       heightfield(size=[100,100], heightfield=[
-//           for (x=[-180:5:180]) [for(y=[-180:5:180]) 10+5*cos(3*x)*sin(3*y)]
+//       heightfield(size=[100,100], data=[
+//           for (y=[-180:5:180]) [for(x=[-180:5:180]) 10+5*cos(3*x)*sin(3*y)]
 //       ]);
 //       cylinder(h=50,d=100);
 //   }
-module heightfield(heightfield, size=[100,100], bottom=0, convexity=10)
+// Example(NORENDER): Heightfield by Function
+//   fn = function (x,y) 10*sin(x*360)*cos(y*360);
+//   heightfield(size=[100,100], data=fn);
+// Example(NORENDER): Heightfield by Function, with Specific Ranges
+//   fn = function (x,y) 2*cos(5*norm([x,y]));
+//   heightfield(size=[100,100], bottom=-20, data=fn, xrange=[-180:2:180], yrange=[-180:2:180]);
+module heightfield(data, size=[100,100], xrange=[-1:0.04:1], yrange=[-1:0.04:1], bottom=-20, maxz=100, style="default", convexity=10, anchor=CENTER, spin=0, orient=UP)
 {
     size = is_num(size)? [size,size] : point2d(size);
-    dim = array_dim(heightfield);
-    assert(dim.x!=undef);
-    assert(dim.y!=undef);
-    assert(bottom<min(flatten(heightfield)), "bottom must be less than the minimum heightfield value.");
-    spacing = vdiv(size,dim-[1,1]);
-    vertices = concat(
-        [
-            for (i=[0:1:dim.x-1], j=[0:1:dim.y-1]) let(
-                pos = [i*spacing.x-size.x/2, j*spacing.y-size.y/2, heightfield[i][j]]
-            ) pos
-        ], [
-            for (i=[0:1:dim.x-1]) let(
-                pos = [i*spacing.x-size.x/2, -size.y/2, bottom]
-            ) pos
-        ], [
-            for (i=[0:1:dim.x-1]) let(
-                pos = [i*spacing.x-size.x/2, size.y/2, bottom]
-            ) pos
-        ], [
-            for (j=[0:1:dim.y-1]) let(
-                pos = [-size.x/2, j*spacing.y-size.y/2, bottom]
-            ) pos
-        ], [
-            for (j=[0:1:dim.y-1]) let(
-                pos = [size.x/2, j*spacing.y-size.y/2, bottom]
-            ) pos
-        ]
-    );
-    faces = concat(
-        [
-            for (i=[0:1:dim.x-2], j=[0:1:dim.y-2]) let(
-                idx1 = (i+0)*dim.y + j+0,
-                idx2 = (i+0)*dim.y + j+1,
-                idx3 = (i+1)*dim.y + j+0,
-                idx4 = (i+1)*dim.y + j+1
-            ) each [[idx1, idx2, idx4], [idx1, idx4, idx3]]
-        ], [
-            for (i=[0:1:dim.x-2]) let(
-                idx1 = dim.x*dim.y,
-                idx2 = dim.x*dim.y+dim.x+i,
-                idx3 = idx2+1
-            ) [idx1,idx3,idx2]
-        ], [
-            for (i=[0:1:dim.y-2]) let(
-                idx1 = dim.x*dim.y,
-                idx2 = dim.x*dim.y+dim.x*2+dim.y+i,
-                idx3 = idx2+1
-            ) [idx1,idx2,idx3]
-        ], [
-            for (i=[0:1:dim.x-2]) let(
-                idx1 = (i+0)*dim.y+0,
-                idx2 = (i+1)*dim.y+0,
-                idx3 = dim.x*dim.y+i,
-                idx4 = idx3+1
-            ) each [[idx1, idx2, idx4], [idx1, idx4, idx3]]
-        ], [
-            for (i=[0:1:dim.x-2]) let(
-                idx1 = (i+0)*dim.y+dim.y-1,
-                idx2 = (i+1)*dim.y+dim.y-1,
-                idx3 = dim.x*dim.y+dim.x+i,
-                idx4 = idx3+1
-            ) each [[idx1, idx4, idx2], [idx1, idx3, idx4]]
-        ], [
-            for (j=[0:1:dim.y-2]) let(
-                idx1 = j,
-                idx2 = j+1,
-                idx3 = dim.x*dim.y+dim.x*2+j,
-                idx4 = idx3+1
-            ) each [[idx1, idx4, idx2], [idx1, idx3, idx4]]
-        ], [
-            for (j=[0:1:dim.y-2]) let(
-                idx1 = (dim.x-1)*dim.y+j,
-                idx2 = idx1+1,
-                idx3 = dim.x*dim.y+dim.x*2+dim.y+j,
-                idx4 = idx3+1
-            ) each [[idx1, idx2, idx4], [idx1, idx4, idx3]]
-        ]
-    );
-    polyhedron(points=vertices, faces=faces, convexity=convexity);
+    vnf = heightfield(data=data, size=size, xrange=xrange, yrange=yrange, bottom=bottom, maxz=maxz, style=style);
+    attachable(anchor,spin,orient, vnf=vnf) {
+        vnf_polyhedron(vnf, convexity=convexity);
+        children();
+    }
 }
 
 
+function heightfield(data, size=[100,100], xrange=[-1:0.04:1], yrange=[-1:0.04:1], bottom=-20, maxz=100, style="default", anchor=CENTER, spin=0, orient=UP) =
+    assert(is_list(data) || is_function(data))
+    let(
+        size = is_num(size)? [size,size] : point2d(size),
+        xvals = is_list(data)
+          ? [for (i=idx(data[0])) i]
+          : assert(is_list(xrange)||is_range(xrange)) [for (x=xrange) x],
+        yvals = is_list(data)
+          ? [for (i=idx(data)) i]
+          : assert(is_list(yrange)||is_range(yrange)) [for (y=yrange) y],
+        xcnt = len(xvals),
+        minx = min(xvals),
+        maxx = max(xvals),
+        ycnt = len(yvals),
+        miny = min(yvals),
+        maxy = max(yvals),
+        verts = is_list(data) ? [
+                for (y = [0:1:ycnt-1]) [
+                    for (x = [0:1:xcnt-1]) [
+                        size.x * (x/(xcnt-1)-0.5),
+                        size.y * (y/(ycnt-1)-0.5),
+                        data[y][x]
+                    ]
+                ]
+            ] : [
+                for (y = yrange) [
+                    for (x = xrange) let(
+                        z = data(x,y)
+                    ) [
+                        size.x * ((x-minx)/(maxx-minx)-0.5),
+                        size.y * ((y-miny)/(maxy-miny)-0.5),
+                        min(maxz, max(bottom+0.1, default(z,0)))
+                    ]
+                ]
+            ],
+        vnf = vnf_merge([
+            vnf_vertex_array(verts, style=style, reverse=true),
+            vnf_vertex_array([
+                verts[0],
+                [for (v=verts[0]) [v.x, v.y, bottom]],
+            ]),
+            vnf_vertex_array([
+                [for (v=verts[ycnt-1]) [v.x, v.y, bottom]],
+                verts[ycnt-1],
+            ]),
+            vnf_vertex_array([
+                [for (r=verts) let(v=r[0]) [v.x, v.y, bottom]],
+                [for (r=verts) let(v=r[0]) v],
+            ]),
+            vnf_vertex_array([
+                [for (r=verts) let(v=r[xcnt-1]) v],
+                [for (r=verts) let(v=r[xcnt-1]) [v.x, v.y, bottom]],
+            ]),
+            vnf_vertex_array([
+                [
+                    for (v=verts[0]) [v.x, v.y, bottom],
+                    for (r=verts) let(v=r[xcnt-1]) [v.x, v.y, bottom],
+                ], [
+                    for (r=verts) let(v=r[0]) [v.x, v.y, bottom],
+                    for (v=verts[ycnt-1]) [v.x, v.y, bottom],
+                ]
+            ])
+        ])
+    ) reorient(anchor,spin,orient, vnf=vnf, p=vnf);
+
 
 // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap
diff --git a/version.scad b/version.scad
index 98921b7e..39cf19ed 100644
--- a/version.scad
+++ b/version.scad
@@ -8,7 +8,7 @@
 //////////////////////////////////////////////////////////////////////
 
 
-BOSL_VERSION = [2,0,500];
+BOSL_VERSION = [2,0,504];
 
 
 // Section: BOSL Library Version Functions