I am working on writing a code to find the polar plane of an ellipse, and subsequently calculate the cap formed by subtending the point $p$ (red) onto a random ellipse.

I think I have managed to compute the polar plane/ellipse of intersection using the algorithm referenced here (https://hal.science/hal-01402302v3/document), though there are some plotting issues potentially, and now I would like to compute the surface area of the ellipsoidal cap on the side closest to the red point (blue point shows the translation of the ellipsoid into coordinate system with the ellipsoid at the origin)

Here is my code so far:

ClearAll[crossProduct, frisvadOrthoBasis, EllipsoidToEllipse]

crossProduct[a_, b_] := {
   a[[2]]*b[[3]] - a[[3]]*b[[2]],
   a[[3]]*b[[1]] - a[[1]]*b[[3]],
   a[[1]]*b[[2]] - a[[2]]*b[[1]]
   };

frisvadOrthoBasis[v0_] := Module[{v = v0},
  (*Normalize input vector*)
  v = Normalize[v];
  (*First orthonormal basis vector*)
  u1 = v;
  (*Second orthonormal basis vector*)
  If[
   Abs[u1[[1]]] > Abs[u1[[2]]],
   u2 = {-u1[[3]], 0, u1[[1]]};,
   u2 = {0, u1[[3]], -u1[[2]]};
   ];
  u2 = Normalize[u2];
  (*Third orthonormal basis vector*)
  u3 = Cross[u1, u2];
  (*Return the orthonormal basis vectors*)
  Return[{u2, u3}];
  ]

EllipsoidToEllipse[A10_, A20_, A30_, a10_, a20_, a30_, Pa0_] := 
 Module[{
   A1 = A10, A2 = A20, A3 = A30, a1 = a10, a2 = a20, a3 = a30, Pa = Pa0
   },
  
  
  (*This function assumes that the ellipsoid lies with its centroid \
at the origin*)
  
  (*3.1 ellipsoid to sphere*)
  M = Transpose[Join[{A1, A2, A3}]] . DiagonalMatrix[{a1, a2, a3}] . 
    Join[{A1, A2, A3}];
  (*Print["M",M];*)
  Minv = Inverse[M];
  Pb = Minv . Pa;
  (*3.2 sphere to disk*)
  theta = ArcSin[1/Norm[Pb]];
  (*Print["theta ",theta];*)
  Pc = (Cos[theta]^2)*Pb;
  (*Print["Pc ",Pc];*)
  radius = Tan[theta]*Norm[Pc];
  (*Print["radius ",radius];*)
  {C1, C2} = frisvadOrthoBasis[Pc/Norm[Pc]];
  (*3.3 disk to ellipse*)
  (*Print["C1 ", C1];*)
  (*Print["C2 ", C2];*)
  Pd = M . Pc;
  D1 = M . (radius*C1);
  (*Print["D1 first",D1];*)
  D2 = M . (radius*C2);
  (*Print["D2 first",D2];*)
  (*3.4 ellipse principal axes*)
  Q = {{Dot[D1, D1], Dot[D1, D2]}, {Dot[D1, D2], Dot[D2, D2]}};
  {eigenvalues, eigenvectors} = Eigensystem[Q];
  (*Print[eigenvectors];*)
  D1 = eigenvectors[[1, 1]]*D1 + eigenvectors[[2, 1]]*D2;
  (*Print["D1 before",D1];*)
  D1 /= Norm[D1];
  (*Print["D1 after",D1];*)
  D2 = eigenvectors[[1, 2]]*D1 + eigenvectors[[2, 2]]*D2;
  (*Print[D2];*)
  D2 /= Norm[D2];
  (*Print["Eigenvalues",eigenvalues];*)
  {d1, d2} = Sqrt[eigenvalues];
  Return[{Pd, D1, D2, d1, d2}];]

Manipulate[
 Show[
  
  Graphics3D[{
    {Red, PointSize[0.05], Point[{px, py, pz}]},
    {Opacity[0.5], Ellipsoid[{0, 0, 0}, {a1, a2, a3}]
     },
    
    
    
    (*Lets plot the point of the plane*)
    {Blue, PointSize[0.05], 
     Point[EllipsoidToEllipse[{1, 0, 0}, {0, 1, 0}, {0, 0, 1}, a1, a2,
          a3, {-px, -py, -pz}][[1]] + {px, py, pz}]},
    
    
    }, Boxed -> False],
  
  
  
  
  

  {
   ContourPlot3D[
    
    {
     Cross[
        Evaluate[(*EllipsoidToEllipse[{1,0,0},{0,1,0},{0,0,1},a1,a2,
         a3,{-px,-py,-pz}][[1]] +{px,py,
         pz}+*)(EllipsoidToEllipse[{1, 0, 0}, {0, 1, 0}, {0, 0, 1}, 
              a1, a2, a3, {-px, -py, -pz}][[3]]*
            
            EllipsoidToEllipse[{1, 0, 0}, {0, 1, 0}, {0, 0, 1}, a1, 
              a2, a3, {-px, -py, -pz}][[
             5]]) - (EllipsoidToEllipse[{1, 0, 0}, {0, 1, 0}, {0, 0, 
               1}, a1, a2, a3, {-px, -py, -pz}][[1]] + {px, py, pz}) ],
        
        Evaluate[(*EllipsoidToEllipse[{1,0,0},{0,1,0},{0,0,1},a1,a2,
         a3,{-px,-py,-pz}][[1]] +{px,py,
         pz}+*)(EllipsoidToEllipse[{1, 0, 0}, {0, 1, 0}, {0, 0, 1}, 
              a1, a2, a3, {-px, -py, -pz}][[2]]*
            
            EllipsoidToEllipse[{1, 0, 0}, {0, 1, 0}, {0, 0, 1}, a1, 
              a2, a3, {-px, -py, -pz}][[4]])
          
          - (EllipsoidToEllipse[{1, 0, 0}, {0, 1, 0}, {0, 0, 1}, a1, 
              a2, a3, {-px, -py, -pz}][[1]] + {px, py, pz})]
        ] . ({x, y, 
          z} - (EllipsoidToEllipse[{1, 0, 0}, {0, 1, 0}, {0, 0, 1}, 
             a1, a2, a3, {-px, -py, -pz}][[1]] + {px, py, pz})) == 0,
     
     
     (x/a1)^2 + (y/a2)^2 + (z/a3)^2 - 1 == 0
     
     }, {x, -0.5, 0.5}, {y, -0.5, 0.5}, {z, -0.6, 0.6}, 
    ContourStyle -> Opacity[0.5]
    
    , MeshFunctions -> Function[
      (*Start of Mesh Intersection contour*)
      {x, y, z},
      (*Start of Mesh Intersection function*)
      Cross[
         
         Evaluate[(*EllipsoidToEllipse[{1,0,0},{0,1,0},{0,0,1},a1,a2,
          a3,{-px,-py,-pz}][[1]] +{px,py,
          pz}+*)(EllipsoidToEllipse[{1, 0, 0}, {0, 1, 0}, {0, 0, 1}, 
               a1, a2, a3, {-px, -py, -pz}][[3]]*
             
             EllipsoidToEllipse[{1, 0, 0}, {0, 1, 0}, {0, 0, 1}, a1, 
               a2, a3, {-px, -py, -pz}][[
              5]]) - (EllipsoidToEllipse[{1, 0, 0}, {0, 1, 0}, {0, 0, 
                1}, a1, a2, a3, {-px, -py, -pz}][[1]] + {px, py, 
              pz}) ],
         
         
         Evaluate[(*EllipsoidToEllipse[{1,0,0},{0,1,0},{0,0,1},a1,a2,
          a3,{-px,-py,-pz}][[1]] +{px,py,
          pz}+*)(EllipsoidToEllipse[{1, 0, 0}, {0, 1, 0}, {0, 0, 1}, 
               a1, a2, a3, {-px, -py, -pz}][[2]]*
             
             EllipsoidToEllipse[{1, 0, 0}, {0, 1, 0}, {0, 0, 1}, a1, 
               a2, a3, {-px, -py, -pz}][[4]])
           
           - (EllipsoidToEllipse[{1, 0, 0}, {0, 1, 0}, {0, 0, 1}, a1, 
               a2, a3, {-px, -py, -pz}][[1]] + {px, py, pz})]
         ] . ({x, y, 
           z} - (EllipsoidToEllipse[{1, 0, 0}, {0, 1, 0}, {0, 0, 1}, 
              a1, a2, a3, {-px, -py, -pz}][[1]] + {px, py, pz}))
       
       - ((x/a1)^2 + (y/a2)^2 + (z/a3)^2 - 1)
      
      (*End of Mesh Intersection function*)
      ]
    , MeshStyle -> {Thick, Red}, Mesh -> {{0}}
    ]}
  
  
  
  
  ]
 
 
 
 ,
 {px, -1, 1},
 {py, -1, 1},
 {pz, -1, 1},
 {{a1, 0.15}, 0.1, 1},
 {{a2, 0.15}, 0.1, 1},
 {{a3, 0.3}, 0.1, 1}
 ]

Thanks for your help and assistance!