Calculate rotation of bones

Hi!
I have the following problem:
I have an armature and for each bone I have set a curve per script to show where the tail is located in every frame. (Every point on the curve represents the position per frame of a bone tail).
Now I want to make it the other way round. I have a curve with points and I want to calculate the quaternion rotation each frame for the bone so that the bone is pointing to the curve point. So I want that the normalized vector v = tail - head = the normalized vector between the point and the head of the bone.
What I have done is the following:


#In each Frame first the Rotation of the bone is set to (1,0,0,0) to calculate the vector between head and tail
bpy.context.scene.frame_set(start_frame + i - 1)
current_frame = start_frame + i - 1
old_rot = armature_object.pose.bones[bone.name].rotation_quaternion
armature_object.pose.bones[bone.name].rotation_quaternion = mathutils.Quaternion((1,0,0,0))
armature_object.pose.bones[bone.name].keyframe_insert(data_path = 'rotation_quaternion', frame = current_frame)
bpy.context.scene.frame_set(start_frame + i - 1)

invert_matrix = world_matrix.copy()
global_tail_position = invert_matrix * armature_object.pose.bones[bone.name].tail
global_head_position = invert_matrix * armature_object.pose.bones[bone.name].head

vec1 = global_tail_position - global_head_position
vec1.normalize()

#This is the point of the curve where the Point should point to
point_pos = point.co

#Der Bone soll in Richtung dieses Vektors zeigen
vec = mathutils.Vector((0,0,0))
vec.x = point_pos.x - global_head_position.x
vec.y = point_pos.y - global_head_position.y
vec.z = point_pos.z - global_head_position.z
vec.normalize()

#Here the rotation is calculated to rotate the bone from vec1 to vec
scalar = vec.dot(vec1)
angle = math.acos(scalar)
tmp = vec1.cross(vec)
tmp.normalize()
                            
k = 0
angle = vec1.angle(vec, k)
angle = angle * 0.5

angle_to_set_quat = mathutils.Quaternion((math.cos(angle), tmp.x * math.sin(angle), tmp.y * math.sin(angle), tmp.z * math.sin(angle)))
rot_matrix = angle_to_set_quat.to_matrix()

global_quat = mathutils.Quaternion((1,0,0,0))
rot_matrix_global = global_quat.to_matrix()

total_rotation = rot_matrix * rot_matrix_global
angle_to_set_quat = total_rotation.to_quaternion()

#The rotation is set
armature_object.pose.bones[bone.name].rotation_quaternion = angle_to_set_quat
armature_object.pose.bones[bone.name].keyframe_insert(data_path = 'rotation_quaternion', frame = current_frame)

But this isn’t working correctly. The result depends on the curve I have. In some cases the rotation is set correctly for half of the curve, but than it kind of rotates back. It seems that only determined angles are allowed in Blender.
Can anyone help me with this problem?