Connecting two points with a line/curve via Python script This problem is absolutely driving me crazy for over a week now. All I need is a bezier curve that starts at point with coordinates (x1, y1, z1,) and ends at a point with coordinates (x2, y2, z2). 

I attempted to write a function using basic coordinate geometry, and it works pretty well 90% of the times, but fails depending on domain/range of numpy.atan. (It creates a bezier curve between the two points, then scales and rotates it appropriately)

def connect_points(v1, v2):
    midpoint = [(1 / 2) * (v1[0] + v2[0]), (1 / 2) * (v1[1] + v2[1]), (1 / 2) *(v1[2] + v2[2])]
    distance = math.sqrt((v2[0] - v1[0]) ** 2 + (v2[1] - v1[1]) ** 2 + (v2[2] - v1[2]) ** 2)
    z_rotation = math.atan((v2[1] - v1[1]) / ((v2[0] - v1[0]) + 0.00001))
    y_rotation = -math.atan((v2[2] - v1[2]) / math.sqrt((v2[1] - v1[1]) ** 2 + (v2[0] - v1[0]) ** 2))
    bpy.ops.curve.primitive_bezier_curve_add(location=midpoint, radius=distance / 2, rotation=[0, y_rotation, z_rotation])
    bpy.context.object.data.bevel_depth = 0.010


I am not willing to spend another week trying to write different cases for arctan. Is there a simpler way?
Using API methods rather than operator.

Using methods outlined similar to those in Bezier-Spline with Python adds unwanted point .

Method to create a 2 point bezier curve object, with each end at defined global locations and the origin in the middle.  Notice the left and right handles types are set to auto resulting in a straight line.. (perhaps a poly would be simpler in this situation)

import bpy
from mathutils import Vector

def add_bezier(v0 , v1):
    v0, v1 = Vector(v0), Vector(v1)  
    o = (v1 + v0) / 2  

    curve = bpy.data.curves.new('Curve', 'CURVE')
    spline = curve.splines.new('BEZIER')
    bp0 = spline.bezier_points[0]
    bp0.co = v0 - o
    bp0.handle_left_type = bp0.handle_right_type = 'AUTO'

    spline.bezier_points.add(count=1)
    bp1 = spline.bezier_points[1]
    bp1.co = v1 - o
    bp1.handle_left_type = bp1.handle_right_type = 'AUTO'
    ob = bpy.data.objects.new('Curve', curve)
    ob.matrix_world.translation = o
    return ob


As an example of usage, select a mesh object and run  script below to add curves matching each edge.

context = bpy.context
scene = context.scene
ob = context.object
mw = ob.matrix_world
me = ob.data

for e in me.edges:
    pts = [mw * me.vertices[v].co for v in e.vertices]
    o = add_bezier(*pts)
    curve = o.data    
    curve.dimensions = '3D'
    curve.bevel_depth = 0.010
    curve.bevel_resolution = 3
    scene.objects.link(o)


Note should you wish to not use the global coordinates, could make a single axis aligned curve, and align it using the methods outlined in Align Object to Vector using python (consider using an end as origin)