Why doesn't my union Boolean modifier work using Python script?

Question

I am new to Blender. I tried to do a union operation on all beams using the python script. As shown in the following image, the model looks like:

However, after the union operation, one of the beams disappears:

I attached my code here:

import bpy;
import numpy as np
import bmesh;
import math;

# Connectivity
Point1List = [[0.0,0.0,0.0],[3.0,0.0,0.0],[2.0,2.0,0.0]]
Point2List = [[3.0,0.0,0.0],[2.0,2.0,0.0],[0.0,0.0,0.0]]

mylayers = [False] * 20;
mylayers[0] = True;

bpy.ops.mesh.primitive_uv_sphere_add(
            location = (0.0, 0.0, 0.0),
            size = 0.4,
            layers=mylayers);

bpy.context.active_object.name = "Point1ListSphere1";

bpy.ops.mesh.primitive_uv_sphere_add(
            location = (3.0, 0.0, 0.0),
            size = 0.4,
            layers=mylayers);

bpy.context.active_object.name = "Point2ListSphere1";

bpy.ops.mesh.primitive_uv_sphere_add(
            location = (2.0, 2.0, 0.0),
            size = 0.4,
            layers=mylayers);

bpy.context.active_object.name = "Point2ListSphere2";

def poly(z,L):
    return   (1.2/(L*L)) * z * z - (1.2/L)*z + 0.4

def rotation_matrix(axis, theta):
    axis = np.asarray(axis)
    axis = axis/math.sqrt(np.dot(axis, axis))
    a = math.cos(theta/2.0)
    b, c, d = -axis*math.sin(theta/2.0)
    aa, bb, cc, dd = a*a, b*b, c*c, d*d
    bc, ad, ac, ab, bd, cd = b*c, a*d, a*c, a*b, b*d, c*d
    return np.array([[aa+bb-cc-dd, 2*(bc+ad), 2*(bd-ac)],
                     [2*(bc-ad), aa+cc-bb-dd, 2*(cd+ab)],
                     [2*(bd+ac), 2*(cd-ab), aa+dd-bb-cc]])  

for i in range(len(Point1List)):
    if Point1List[i] != Point2List[i]:          
        Dist = 0.0
        for j in range(3):
            Dist+=(Point1List[i][j]-Point2List[i][j])**2

        Dist=math.sqrt(Dist)

        # Calculate ratational axis
        u = [0.0,0.0,Dist]
        v = [Point2List[i][0]-Point1List[i][0],Point2List[i][1]- 
            Point1List[i][1],Point2List[i][2]-Point1List[i][2]]
        R_axis = np.cross(u,v)
        R_axis = R_axis/np.linalg.norm(R_axis)

        # Calculate theta radians
        theta = math.acos(np.dot(u,v)/(np.linalg.norm(u)*np.linalg.norm(v)));

        # Create a new object with name
        mesh = bpy.data.meshes.new("Branch"+str(i+1))
        object = bpy.data.objects.new("Branch"+str(i+1), mesh)

        # Show the object in the scene
        scene = bpy.context.scene
        scene.objects.link(object)  # put the object into the scene (link)
        scene.objects.active = object  # set as the active object in the scene

        # Set the layer
        object.layers[0] = True
        for j in range(1,20):
            object.layers[j] = False

        object.select = True  # select object
        mesh = bpy.context.object.data
        bm = bmesh.new()

        # Translate and rotate the created column to the correct position

        v_coord=[poly(0,Dist), 0, 0]
        R_matrix = rotation_matrix(R_axis,theta)
        v_coord = np.dot(R_matrix,v_coord)
        # Point translation
        v_coord = v_coord + Point1List[i]

        v0 = bm.verts.new((v_coord[0], v_coord[1], v_coord[2]))

        for dz in range(1, 1000):
            z = dz / (1000/Dist)
            if z>=Dist:
                break
            v_coord=[poly(z,Dist), 0, z]

            v_coord = np.dot(R_matrix,v_coord)
            v_coord = v_coord + Point1List[i]
            v_new = bm.verts.new((v_coord[0], v_coord[1], v_coord[2]))


            edge = bm.edges.new((v0, v_new))
            v0 = v_new

        bmesh.ops.spin(bm, geom=bm.edges, axis=(v[0], v[1], v[2]), angle=math.radians(360), steps=20, cent=(Point1List[i][0],Point1List[i][1],Point1List[i][2])) 
        # Remove the duplicate vertices
        bmesh.ops.remove_doubles(bm, verts=bm.verts, dist=0.001)

        bm.to_mesh(mesh)
        mesh.update()

# Boolean operation
def UnionBoolean(target, opObj, i):

    bpy.ops.object.select_all(action='DESELECT')
    target.select = True

    bpy.context.scene.objects.active = target
    bpy.ops.object.modifier_add(type='BOOLEAN')

    mod = target.modifiers
    mod[0].name = "Boon"+str(i)
    mod[0].object = opObj
    mod[0].operation = 'UNION'

    bpy.ops.object.modifier_apply(apply_as='DATA', modifier=mod[0].name)

for i in range (1,len(Point1List)):
    target = bpy.data.objects['Branch1']
    opObj = bpy.data.objects['Branch'+str(i+1)]

    UnionBoolean(target,opObj, i)

Is there any problem in my code that makes the boolean not working? I used remove double. Normals not pointing outward could be the issue but I don't know how to change normals using python script. Thank you.