* For the ellipsoid equation 
$$\frac{(x-\text{p1})^2}{\text{a1}^2}+\frac{(y-\text{p2})^2}{\text{a2}^2}+\frac{(z-\text{p3})^2}{\text{a3}^2}=1$$
and a point `{o1,o2,o3}`, it's polar plane just has the form 

$$\frac{(\text{o1}-\text{p1}) (x-\text{p1})}{\text{a1}^2}+\frac{(\text{o2}-\text{p2})
   (y-\text{p2})}{\text{a2}^2}+\frac{(\text{o3}-\text{p3}) (z-\text{p3})}{\text{a3}^2}=1$$

The main reasion is that the normal vector of the ellipsoid (the normal vector of the tangent plane of such ellipsoid) is

```
Clear["Global`*"];
f[x_, y_,z_] := (x - p1)^2/a1^2 + (y - p2)^2/a2^2 + (z - p3)^2/a3^2 - 1;
Grad[f[x, y, z], {x, y, z}]
```
$$\left\{\frac{2 (x-\text{p1})}{\text{a1}^2},\frac{2 (y-\text{p2})}{\text{a2}^2},\frac{2
   (z-\text{p3})}{\text{a3}^2}\right\}$$

and the direction of the tangent lines is orthogonal to such normal vector.

```
Clear["Global`*"];
{p1, p2, p3} = {0.2, .4, .8};
{a1, a2, a3} = {.8, .4, .5};
{o1, o2, o3} = {0, 0, 0};
f[x_, y_, z_] := (x - p1)^2/a1^2 + (y - p2)^2/a2^2 + (z - p3)^2/
   a3^2 - 1;
plane[x_, y_, z_] := ((x - p1) (o1 - p1))/a1^2 + ((y - p2) (o2 - p2))/
   a2^2 + ((z - p3) (o3 - p3))/a3^2 - 1;
ellipsoid = 
  ContourPlot3D[
   f[x, y, z] == 0 // Evaluate, {x, p1 - a1, p1 + a1}, {y, p2 - a2, 
    p2 + a2}, {z, p3 - a3, p3 + a3}, 
   MeshFunctions -> Function[{x, y, z}, plane[x, y, z]], 
   Mesh -> {{0}}, 
   MeshShading -> {Directive@{Opacity[.5], Green}, 
     Directive@{Opacity[.5], Blue}}, 
   MeshStyle -> Directive@{Thick, Yellow}, BoxRatios -> Automatic, 
   Boxed -> False, Axes -> False];
polarplane = 
  ContourPlot3D[
   plane[x, y, z] == 0 // Evaluate, {x, p1 - a1, p1 + a1}, {y, 
    p2 - a2, p2 + a2}, {z, p3 - a3, p3 + a3}, Mesh -> None];
point = Graphics3D[{Red, AbsolutePointSize[5], Point@{o1, o2, o3}}];
Show[ellipsoid, polarplane, point, PlotRange -> All]
```

[![enter image description here][1]][1]

and the areas of the two parts (blue and green) are

```
area1 = ImplicitRegion[{f[x, y, z] == 0, plane[x, y, z] >= 0}, {x, y, 
     z}] // Area;
area2 = ImplicitRegion[{f[x, y, z] == 0, plane[x, y, z] <= 0}, {x, y, 
     z}] // Area;
{area1, area2}
```

> `{0.970791, 2.96791}`

* Such method also work for the the 2D cases. In the 2D cases just the poles of quadratic curves.


  [1]: https://i.sstatic.net/HibX5qOy.png