I needed to do this recently. I think the following code should do the trick.

`p` is our point in 3D space.

    camera = bpy.data.cameras["Camera"]
    scene = bpy.context.scene
    
    f = camera.cycles.fisheye_lens
    w = camera.sensor_width
    h = (3/4) * w # Depends on aspect ratio. For some reason I don't get the correct value from camera.sensor_height
    
    camera_ob = bpy.data.objects["Camera"]
    
    p = camera_ob.matrix_world.inverted() * p
    p.normalize()

    # Calculate our angles
    phi = math.atan2(p.y, p.x)
    l = (p.x**2 + p.y**2)**(1/2)
    theta = math.asin(l)

    # Equisolid projection
    r = 2.0 * f * math.sin(theta / 2)

    u = r * math.cos(phi) / w + 0.5
    v = r * math.sin(phi) / h + 0.5

    x = u * scene.render.resolution_x
    y = v * scene.render.resolution_y