The Earth is assumed to be a WGS84 ellipsoid (oblate spheroid) $E$. $E: (x^2+y^2)/R_{eq}^2+z^2/R_{pol}^2=1$ With $R_{eq} > R_{pol}$, a point $M$ is outside $E$, and $P$ is its polar plane, that intersects $E$, forming an ellipse $S$.
How to get the coordinates of this $S$ in parametric form or in spherical form?
My goal is to iterate over the parameters of $S$ to get the coordinates of all the points forming $S$.
I tried many methods, adviced from others' help, but they all complicated
For a spherical Earth defined by $x^2+y^2+x^2=R^2$ and a point $M$ at $(D,0,0)$ the polar plane is defined by $x=d=\frac{R^2}D$, which intersects the sphere in a circle with radius $\rho=R\sqrt{1-(\frac RD)^2}$ that can be parametrized by $(x,y,z)=(d,\rho \cos t,\rho \sin t)$ with $0\le t \lt 2\pi$.
With point $M$ rotated upwards to $(D\cos\phi,0,D\sin\phi)$ the parametrization for the similarly rotated circle becomes $(x,y,z)=(d\cos\phi-\rho\sin\phi\sin t,\rho\cos t,d\sin \phi+\rho\cos\phi\sin t)$
A final rotation around the $z$-axis relocates point $M$ to $(D\cos\lambda\cos\phi,D\sin\lambda\cos\phi,D\sin\phi)$ and the circle to points $(d\cos\lambda\cos\phi-\rho\cos\lambda\sin\phi\sin t-\rho\sin\lambda\cos t,d\sin\lambda\cos\phi-\rho\sin\lambda\sin\phi\sin t+\rho\cos\lambda\cos t,d\sin \phi+\rho\cos\phi\sin t)$.
So for the ellipsoid case one could first transform $M=(x,y,z)$ to $(x,y,\frac{R_{eq}}{R_{pol}}z)$, find the circle parametrisation for a sphere with radius $R_{eq}$ and multiply its $z$-component with $\frac{R_{pol}}{R_{eq}}$.