Not exactly based on the example given, but from the description it sounds like you want to create a bezier curve with control points along normals:
import bpy
from bpy import context as C, data as D
strength = 5
num_points = 20
if C.mode != 'EDIT_MESH':
raise Exception("Select two objects, go to Edit Mode and select a face for each")
bpy.ops.object.mode_set(mode='OBJECT')
ob1, ob2 = C.selected_objects
f1, f2 = (f for ob in (ob1, ob2) for f in ob.data.polygons if f.select)
m1, m2 = (o.matrix_world for o in C.selected_objects)
curve_data = D.curves.new('connector', type='CURVE')
curve_data.dimensions = '3D'
curve_data.resolution_u = num_points
spline = curve_data.splines.new('BEZIER')
spline.bezier_points.add(1)
start, end = spline.bezier_points
start.handle_left = m1 @ f1.center
start.co = start.handle_left
start.handle_right = m1 @ (f1.center + f1.normal*strength)
end.handle_left = m2 @ (f2.center + f2.normal*strength)
end.co = m2 @ f2.center
end.handle_right = end.co
curve_ob = D.objects.new('connector', curve_data)
C.collection.objects.link(curve_ob)
# bpy.ops.object.mode_set(mode='EDIT')

You say you want to generate "a number of points", which you can do by converting the curve to the mesh, though I'm not sure if this is something you really need, as you say later you want to convert it to a curve anyway... You can add to the code:
# Convert to a mesh:
ob1.select_set(False)
ob2.select_set(False)
C.view_layer.objects.active = curve_ob
curve_ob.select_set(True)
bpy.ops.object.convert(target='MESH')
Mathematics
The maths are as simple as defining two lines and interpolating (lerping) between them, e.g. define a line going from the center of face1 towards its normal, and another line from the center of face2 towards its normal. Then just go away from face1 on each point, and towards face2 on each point on those lines, and interpolate between the results. This is pretty much what the quadratic bezier does, except in its formula both lines point to the same point which is defined by the control point, so it's just a matter of when you calculate the control point, but the result is the same. This quadratic bezier is similar to the example given by OP. For a cubic bezier, just interpolate between two quadratic beziers. For higher order beziers interpolate between lower order beziers...
import bpy
from bpy import context as C, data as D
strength = 5
num_points = 20
if C.mode != 'EDIT_MESH':
raise Exception("Select two objects, go to Edit Mode and select a face for each")
bpy.ops.object.mode_set(mode='OBJECT')
ob1, ob2 = C.selected_objects
f1, f2 = (f for ob in (ob1, ob2) for f in ob.data.polygons if f.select)
m1, m2 = (o.matrix_world for o in C.selected_objects)
co1, n1 = f1.center, f1.normal * strength
co2, n2 = f2.center, f2.normal * strength
verts = []
edges = []
last_i = num_points - 1
def lerp_between_straight_lines(i):
weight1 = i/last_i
weight2 = (last_i - i)/last_i
along_n1 = co1 + n1*weight2
along_n2 = co2 + n2*weight1
return along_n1 * weight1 + along_n2 * weight2
def quad_bezier(i):
t = i/last_i # same as weight1
# weight2 became (1-t)
co = (co1 + co2) / 2
n = (n1 + n2) / 2
control_point = co + n
# again two straight lines, but this time not along normals but towards the control point
from_co1_to_co = (1-t)*co1 + t*control_point
from_co2_to_co = t*co2 + (1-t)*control_point
return (1-t)*from_co1_to_co + t*from_co2_to_co
def cubic_bezier(i):
t = i/last_i
cp1 = co1 + n1
cp2 = co2 + n2
from_co1_to_cp1 = (1-t)*co1 + t*cp1
from_cp2_to_cp1 = t*cp2 + (1-t)*cp1
quad_bezier_1 = (1-t)*from_co1_to_cp1 + t*from_cp2_to_cp1
from_cp1_to_cp2 = (1-t)*cp1 + t*cp2
from_co2_to_cp2 = t*co2 + (1-t)*cp2
quad_bezier_2 = (1-t)*from_cp1_to_cp2 + t*from_co2_to_cp2
return (1-t)*quad_bezier_1 + t*quad_bezier_2
for i in range(num_points):
p = lerp_between_straight_lines(i)
# p = quad_bezier(i)
# p = cubic_bezier(i)
verts.append(p.to_tuple())
if i > 0:
edges.append((i-1, i))
me = D.meshes.new('connector')
me.from_pydata(verts, edges, [])
ob = D.objects.new('connector', me)
C.collection.objects.link(ob)