How to compute the solid angle underlying the intersection of a cone and a sphere?

Question

I want to compute the solid angle underlying the intersection between a cone and a sphere. To do that, I am using Blender to ease the calculations.

I parametrise a cone by specifying its apex location apex, half-angle opening halfOpeningDeg, and axis orientation with azimuth and elevation. I parametrise a sphere by specifying its centre location centreShpere and radius radiusSphere. As an example, here is a snippet creating the objects (the modules Cone and Sphere implement a class of the same name, with the Blender objects created upon instanciation)

from Cone import Cone
from Sphere import Sphere

apex: tuple[float, float, float] = (0.0, 0.0, 0.0)
height: float = 1.0e1
halfOpeningDeg: float = 15.0
azimuth: float = 45.0
elevation: float = 15.0
cone: Cone = Cone(
                apex = apex,
                height = height,
                halfOpeningDeg = halfOpeningDeg,
                azDeg = azimuth,
                elDeg = elevation,
            )

centreSphere: tuple[float, float, float] = (5.0, 3.0, 0.0)
radiusSphere: float = 3.5
sphere: Sphere = Sphere(
                        radius = radiusSphere,
                        centre = centreSphere,
                    )

Based on those objects, I want to compute the solid angle underlying their intersection. To do so, I have the following functions

import bpy
import mathutils
import bmesh
import numpy as np
from Cone import Cone
from Sphere import Sphere

def check_intersection(cone: Cone,
                        sphere: Sphere,
                    ) -> bool:
    """
    Function checking if a (Blender) cone intersects a (Blender) sphere

    Parameters
    ----------
    cone : bpy.types.Object
        Blender cone
    sphere : bpy.types.Object
        Blender sphere

    Return
    ------
    - : bool
        Condition checking if the cone and sphere intersect
    """
    coneBlend: bpy.types.Object = cone.cone
    sphereBlend: bpy.types.Object = sphere.sphere

    # Get the cone's axis orientation
    coneAxis: mathutils.Vector = mathutils.Vector((0.0, 0.0, 1.0))
    coneAxis.rotate(coneBlend.rotation_quaternion)

    # Project vector onto the cone axis
    apexToCentre: mathutils.Vector = sphereBlend.location - coneBlend.location
    projLength: float = apexToCentre.dot(coneAxis)
    closestPtAxis: mathutils.Vector = coneBlend.location + projLength * coneAxis

    # Calculate the distance from the sphere centre to the cone axis
    distToAxis: float = (sphereBlend.location - closestPtAxis).length

    # Get the radius of the sphere
    sphereRadius: float = (0.5 * sphereBlend.dimensions.x) * sphereBlend.scale.x

    # Calculate the cone radius at the height of the closest point on the axis
    coneDims: mathutils.Vector = coneBlend.dimensions
    coneScales: mathutils.Vector = coneBlend.scale
    coneRadiusHeight = (projLength / (coneDims[2] * coneScales[2])) * (0.5 * coneScales[0]) * coneScales[0]

    return distToAxis < (sphereRadius + coneRadiusHeight)

def calculate_intersection_volume(cone: Cone,
                                    sphere: Sphere,
                                ) -> float:
    """
    Function computing the volume of intersection between a (Blender) cone and a
    (Blender) sphere

    Parameters
    ----------
    cone : bpy.types.Object
        Blender cone
    sphere : bpy.types.Object
        Blender sphere

    Return
    ------
    volume : float
        Volume of intersection between the (Blender) cone and (Blender) sphere
        (if any)
    """
    bpy.ops.object.select_all(action = 'DESELECT')

    coneBlend: bpy.types.Object = cone.cone
    sphereBlend: bpy.types.Object = sphere.sphere

    # Duplicate the cone object
    cone_copy = coneBlend.copy()
    cone_copy.data = coneBlend.data.copy()
    bpy.context.collection.objects.link(cone_copy)

    # Select the cone copy and sphere
    cone_copy.select_set(True)
    # sphere.select_set(True)
    bpy.context.view_layer.objects.active = cone_copy
    bpy.ops.object.modifier_add(type = 'BOOLEAN')
    bpy.context.object.modifiers["Boolean"].operation = 'INTERSECT'
    bpy.context.object.modifiers["Boolean"].object = sphereBlend

    # Apply the boolean modifier
    bpy.ops.object.modifier_apply(modifier = "Boolean")

    # Get the intersection object
    intersection: bpy.types.Object = bpy.context.active_object
    
    # Use volume calculation if required
    bm: bmesh.types.BMesh = bmesh.new()
    bm.from_mesh(intersection.data)
    volume: float = bm.calc_volume()
    bm.free()

    # Clean up: remove the intersection object
    bpy.ops.object.select_all(action = 'DESELECT')
    intersection.select_set(True)
    bpy.ops.object.delete()
    
    return volume

def calculate_solid_angle(cone: Cone, sphere: Sphere) -> float:
    """
    Compute the solid angle underlying the intersection of a cone and a sphere

    Parameters
    ----------
    cone : Cone
        Cone
    sphere : bpy.types.Object
        Blender sphere

    Return
    ------
    solidAngIntersec : float
        Solid angle of the intersection, in steradians
    """
    # Calculate the cone's half-opening angle
    coneSolidAng: float = cone.get_solidAngle()

    coneBlend: bpy.types.Object = cone.cone
    sphereBlend: bpy.types.Object = sphere.sphere
    
    # Sphere and cone properties
    sphereRadius: float = (0.5 * sphereBlend.dimensions.x) * sphere.sphere.scale.x
    coneHeight: float = coneBlend.dimensions.z * coneBlend.scale.z

    # Distance from cone apex to the base
    halfOpeningRad: float = np.radians(cone.halfOpeningDeg)
    apexBaseDist: float = coneHeight / np.cos(halfOpeningRad)

    # # Distance from cone apex to the intersection along the axis
    distToIntersec: float = (sphereBlend.location - coneBlend.location).length

    if (distToIntersec + apexBaseDist) <= sphereRadius:
        # If the entire cone is within the sphere, use the solid angle of the cone
        return coneSolidAng
    else:
        # Fraction of the cone height inside the sphere
        distToIntersec -= sphereRadius  # Adjust distance to intersection
        fraction_inside_sphere: float = max(0.0, min(1.0, (sphereRadius - distToIntersec) / coneHeight))

        # Calculate the solid angle of the intersection
        solid_angle_intersection: float = coneSolidAng * fraction_inside_sphere
        return solid_angle_intersection

Running that script gives the output shown here below

As can be seen, the cone intersects "completely" the sphere, hence the solid angle underlying the intersection should be the same as that of the cone, but it is not. I do not get what is wrong in my calculations, can someone help me out?