I am stuck trying to make a cylinder given its endpoint locations and a radius.  I need to compute the rotation for the following function call: 

bpy.ops.mesh.primitive_cylinder_add(radius = radius, depth = depth, location = center, rotation = rotation)

How do I do this?  For example, given endpoints (0,0,0) and (2,2,2), what value do I use for "rotation"?  Thanks for any pointers.  

Tough one, eh?  :-)  I've googled all over the place, and I've found some code, but it does not work. that is, the resulting quaternion, when used as obj.rotate_quaternion, does not result in properly rotated cylinder with ep1 = (0,0,0) and ep1 = (0,100,0)

def FindRotation(ep1, ep2):

First, create quaternion per

http://stackoverflow.com/questions/1171849/finding-quaternion-representing-the-rotation-from-one-vector-to-another

vec = numpy.array(ep2) - numpy.array(ep1)

vec = vec/numpy.linalg.norm(vec)

if not vec[0] and not vec[1]:
return (0,0,0)
print ("vec: %s"%str(vec))
Z = numpy.array((0,0,1))
kcostheta = numpy.dot(Z, vec)
k = Length(vec) # sqrt(Length(vec)^2 * Length(Z)^2) = Length(vec)
Q = Quaternion(numpy.cross(vec,Z), kcostheta+k)
Q.normalize()
print ("kcostheta: %f, k: %f, Q: %s"%(kcostheta, k, str(Q)))
return Q

Found a solution thanks to  http://mcngraphics.com/thelab/blender/connect/

def FindRotation(ep1, ep2):
ep1 = numpy.array(ep1)
ep2 = numpy.array(ep2)
center = ep1 + 0.5 * (ep2 - ep1)

Position of ep2 if an imaginary sphere encompassing

both objects was transformed to 0, 0, 0

ep2new = ep2 - center

--------------------------------------------

Spherical coordinates (radius r, inclination theta, azimuth phi) from

the Cartesian coordinates (x, y, z) of o1 on our imaginary sphere

r = numpy.linalg.norm(ep2new) # radius of sphere
theta = math.acos(ep2new[2]/r)
phi = math.atan2( ep2new[1], ep2new[0])
return (0,theta, phi)