Q:

Script method to trim mesh objects

I went away and learned some OO in python, and now I will pick up where I left off - creating mesh objects from mathematical expressions. In the simple test below, I would like to cleanly trim the tops and bottoms of the surfaces where they would intersect planes at Z = +1 and Z = -1. So far I just append faces which have at least one vertex inside. (there are probably unused vertices which I should clean up at some point).
I am hoping there is a simple way to trim this (or future) object at Z = +/-1 within the script. There are probably a number of possible ways to do this - what would be most helpful is both a few lines of script and a pointer to the relevant documentation on mesh object methods so I can dig in and learn. Any other pointers/suggestions would be helpful. I hope to plot some pretty complicated stuff, so helping me starting off on the right track would be greatly appreciated!
Note: An alternative approach would be to calculate the exact intersections myself, but for math more complicated than just np.tan() I'd have to use something like 
from scipy import fsolve

to find the intersection with the planes, and that's another issue.
import bpy
import numpy as np

pi = np.pi

nx, ny = 40, 21
x = pi * np.linspace(-0.475, 1.475, nx)
y = pi * np.linspace( 0.0,   2.0,   ny)
X, Y = np.meshgrid(x, y)
Z =  np.tan(X + Y)

verts = []    
for iy in range(ny):
    for ix in range(nx):
        vert = (X[iy,ix], Y[iy,ix], Z[iy,ix])
        verts.append(vert)

lo, hi = -1.0, 1.0

faces = []
for iy in range(ny-1):
    for ix in range(nx-1):
        i = ix + nx*iy
        face = (int(i), int(i+1), int(i+nx+1), int(i+nx))
        highpoint = max([verts[iv][2] for iv in face])   
        lowpoint  = min([verts[iv][2] for iv in face])
        # if face has at least one vertex inside hi-lo range, 
        # but doesn't cross Z=0 (by more than a little)
        if (highpoint > lo and lowpoint < hi) and highpoint*lowpoint >= -0.1: 
            faces.append(face)

# alternatively could I remove them later using some things such as
# faces.pop() and face.co[:][2]???

meshthing = bpy.data.meshes.new("Tan")
object = bpy.data.objects.new("Tan", meshthing)

object.location = (0,0,0)  #bpy.context.scene.cursor_location
bpy.context.scene.objects.link(object)

meshthing.from_pydata(verts,[],faces)
#meshthing.update(calc_edges=True)

mat = bpy.data.materials.new("PKHG")
mat.diffuse_color = (float(1.5),0.0,1.0)
mat.specular_color = (200.0,1.0,100.0)
object.active_material = mat

mypolys = meshthing.polygons
for p in mypolys:
    p.use_smooth = True

# this helps ensure it's visible and selected
bpy.data.objects['Tan'].select = False
bpy.data.objects['Tan'].select = True

print_info = False
for item in bpy.data.objects:
    print(item.name)
    if item.type == 'MESH' and print_info: 
        for vertex in item.data.vertices:
            print(vertex.co)
        # This doesn't work yet
        for face in item.data.tessfaces:
            print(face.co)

A:

The simple way to solve your problem is to use a Boolean modifier and intersect with a really big box spanning z=-1..1 , but eventually you will want to do something to a mesh that isn't so simple.
You're going to want the bmesh module:
http://www.blender.org/documentation/blender_python_api_2_67_release/bmesh.html
It has access to the vertices, edges, and faces of the mesh with the ability to add, modify, and remove any of them.
For an example I'll include the code from http://web.purplefrog.com/~thoth/blender/python-cookbook/squish-mesh-axis.html :
import bmesh
import bpy

def stretchAxis(mesh, idx, tgt_min, tgt_max):
    """ Squish/stretch the vertex coordinates of a mesh so that their .co[idx] ranges from tgt_min to tgt_max. """
    bm = bmesh.new()
    bm.from_mesh(mesh)

    min = bm.verts[0].co[idx]
    max = min

    for v in bm.verts:
        val = v.co[idx]
        if (val<min):
            min = val
        if (val > max):
            max = val

    # remap co[idx] from min..max to tgt_min..tgt_max
    for v in bm.verts:
        old = v.co[idx]
        v.co[idx] = (old-min)*(tgt_max-tgt_min) / (max-min) +tgt_min
    bm.to_mesh(mesh)
    bm.free()

#
# squish or stretch all the selected meshes to a specific thickness along a
# specific axis.
#

idx = 2 # the z coordinate
for obj in bpy.context.scene.objects:
    if (obj.select):
        stretchAxis(obj.data, idx, 0, 1)
        obj.location[idx] = 1
        obj.scale[idx] = 1

Although this example does not alter the TOPOLOGY (faces and edges).  It only adjusts the position of the vertices.

