I am trying to make a conical frustum follow an ellipsoid.
To do so, I defined the following data in Blender

import bpy
import mathutils
import numpy as np

# ELLIPSOID
# =========
ce = mathutils.Vector((2.5, 2.5, 2.5))
l = mathutils.Vector((1.5, 1.0, 0.75))
nhor: int = 60
nver: int = 60
bpy.ops.mesh.primitive_uv_sphere_add(segments = nhor, ring_count = nver, radius = 1.0, location = ce)
ell = bpy.context.object
ell.scale = (l.x, l.y, l.z)
ell.name = "Ell"


# CONICAL FRUSTUM
# ================
cf = mathutils.Vector((0.0, 0.0, 0.0))
fru_direction = ce - cf
fru_direction.normalize()
r1: float = 0.15
h: float = (ce - cf).length
div: float = np.pi / 60.0
r2: float = r1 + h * div
bpy.ops.mesh.primitive_cone_add(vertices = 32, radius1 = r1, radius2 = r2, depth = h)
fru = bpy.context.object
fru.name = "FRU"

up = mathutils.Vector((0.0, 0.0, 1.0))
rotation_quat = up.rotation_difference(fru_direction)
fru.rotation_mode = 'QUATERNION'
fru.rotation_quaternion = rotation_quat

fru.location = cf + fru_direction * (0.5 * fru.dimensions.z)
bpy.context.view_layer.objects.active = fru
bpy.context.scene.cursor.location = cf
bpy.ops.object.origin_set(type = 'ORIGIN_CURSOR', center = 'BOUNDS')
bpy.context.scene.cursor.location = cursor_origin

I then defined some key positions at certain frames for the ellipsoid centre position

# ANIMATION STUFF
# ===============
ell_positions = [
                ce.to_tuple(),
                (2.0, 2.0, 2.0),
                (4.0, 0.0, 4.0),
                (6.0, -2.0, 2.0),
                (8.0, 0.0, 0.0)
            ]
frames = [0, 20, 40, 60, 80]

The following snippet animates the bodies, and the frustum follows well the centre of the ellipsoid.
However, the scaling is not good, and I do not understand why?

for frame, ell_position in zip(frames, ell_positions):
    bpy.context.scene.frame_set(frame)
    ell.location = ell_position
    ell.keyframe_insert(data_path = 'location', frame = frame)
    fru_direction = mathutils.Vector(ell_position) - cf
    fru_direction.normalize()
    fru.rotation_quaternion = up.rotation_difference(fru_direction)
    h = (ell.location - cf).length
    r2 = r1 + h * div
    fru.scale = (r2, r2, h / fru.dimensions.z)
    fru.keyframe_insert(data_path = 'rotation_quaternion', frame = frame)
    fru.keyframe_insert(data_path = 'scale', frame = frame)

bpy.context.scene.frame_start = frames[0]
bpy.context.scene.frame_end = frames[-1]