# Intersection of oblate spheroid $\frac{x^2+y^2}{R_e^2}+\frac{z^2}{R_p^2} =1$ and plane $n_xx+n_yy+n_zz=0$

To calculate the distance between two points on Earth, I used 3 different approaches.

For small distances, I used the Euclidean distance. For medium distances, I used the arc length on the circle obtained by intersecting the spherical earth with the plane passing through the origin and the two points.

For large distances, I assumed earth to be an oblate spheroid. So, I have to calculate an elliptic arc length. The angle (I assumed) stays the same as a spherical Earth, but I need the axes of the ellipse. The problem becomes as follows.

I have an oblate spheroid $$\frac{x^2+y^2}{R_e^2}+\frac{z^2}{R_p^2} =1$$ and a plane passing through the origin $$n_xx+n_yy+n_zz=0$$.

Their intersection is an ellipse. I want to calculate the semi-major and semi-minor axes of the ellipse.

I am only interested in the two axes.

I found here that the semi-major axis is the same. What about the semi-minor axis?

EDIT:

I replaced $$z = -\frac{n_xx+n_yy}{n_z}$$ in the oblate spheroid equation and I got:

$$\frac{x^2+y^2}{R_e^2} + \frac{1}{R_p^2 n_z^2}\left( n_x^2x^2 + n_y^2 y^2 + 2n_xn_yxy \right)= 1$$

$$\left( \frac{1}{R_e^2} + \frac{n_x^2}{R_p^2 n_z^2} \right)x^2 + \left( \frac{1}{R_e^2} + \frac{n_y^2}{R_p^2 n_z^2} \right)x^2 + \frac{2n_xn_y}{R_p^2 n_z^2}xy = 1$$

Now, I have to rotate the ellipse to lie in the $$xy$$ plane.

The expression will involve elliptic integrals, which cannot be expressed in terms of elementary functions. In a computer program, you will most likely have to evaluate it numerically.

Instead, I would recommend combining the approaches: approximate the curve with short line segments.

First, note that the expression for the surface of the oblate spheroid can be rewritten as $$(R_p x)^2 + (R_p y)^2 + (R_c z)^2 = (R_c R_p)^2$$ If you have some vector $$(x, y, z)$$, you can use the above to scale it (multiply each of its coordinate) by $$c$$ to stretch or shrink it to the surface of the oblate spheroid: $$c = \sqrt{\frac{R_c^2 R_p^2}{(R_p x)^2 + (R_p y)^2 + (R_c z)^2}}$$ This is useful, because you can use the latitude and longitude on an unit sphere to give you the direction, then scale that direction by $$c$$ to get a point on the oblate spheroid.

Second, if you have point $$(x_0, y_0, z_0)$$ on the unit sphere ($$x_0^2 + y_0^2 + z_0^2 = 1$$) corresponding to the latitude and longitude of the starting point on the oblate spheroid, and point $$(x_1, y_1, z_1)$$ on the unit sphere corresponding to the latitude and longitude of the finishing point on the oblate spheroid, you can calculate the axis of rotation $$(x_n, y_n, z_n)$$ between these two points using \left\lbrace ~ \begin{aligned} x_n &= y_0 z_1 - z_0 y_1 \\ y_n &= z_0 x_1 - x_0 z_1 \\ z_n &= x_0 y_1 - y_0 x_1 \\ \end{aligned} \right . which is just a cross product between the two directions. Note that you'll want to normalize this to unit length, \left\lbrace ~ \begin{aligned} X_n &= \displaystyle \frac{x_n}{x_n^2 + y_n^2 + z_n^2} \\ Y_n &= \displaystyle \frac{y_n}{x_n^2 + y_n^2 + z_n^2} \\ Z_n &= \displaystyle \frac{z_n}{x_n^2 + y_n^2 + z_n^2} \\ \end{aligned} \right. The angular distance between the two directions (assuming unit lengths, $$x_0^2 + y_0^2 + z_0^2 = 1$$ and $$x_1^2 + y_1^2 + z_1^2 = 1$$) is $$\theta$$, $$\cos\theta = x_0 x_1 + y_0 y_1 + z_0 z_1$$ i.e. $$\theta = \arccos(x_0 x_1 + y_0 y_1 + z_0 z_1)$$

If you want to split the path into $$N$$ line segments, you'll need $$N+1$$ points along the curve on the oblate spheroid between the starting and ending points. These line segments will vary slightly in length, as the oblateness of the spheroid will affect things, but it works very well for approximating the curve length. To get point $$k = 0, 1, \dots, N$$, use $$\varphi_k = \frac{k \theta}{N}$$ and apply Rodrigues' rotation formula to the starting point (on unit sphere) to find the $$k$$'th direction, then scale it to the oblate spheroid. In Cartesian coordinate form, the rotation is \left\lbrace ~ \begin{aligned} x_k &= x_0 \cos(\varphi_k) + (y_n z_0 - z_n y_0) \sin(\varphi_k) + x_n P \bigr( 1 - \cos(\varphi_k) \bigr) \\ y_k &= y_0 \cos(\varphi_k) + (z_n x_0 - x_n z_0) \sin(\varphi_k) + y_n P \bigr( 1 - \cos(\varphi_k) \bigr) \\ z_k &= z_0 \cos(\varphi_k) + (x_n y_0 - y_n x_0) \sin(\varphi_k) + z_n P \bigr( 1 - \cos(\varphi_k) \bigr) \\ \end{aligned} \right. where $$P$$ is a constant (dot product between the starting point and the rotation axis), $$P = x_0 x_n + y_0 y_n + z_0 z_n$$

Note that since $$\theta$$ is the angular distance, you can make a very sensible choice for $$N$$ using $$N = \left\lceil \frac{\theta}{\phi} \right\rceil$$ where $$\lceil~\rceil$$ denote rounding up (ceil()) and $$\phi$$ is the largest angular distance you wish to approximate with a straight line. For very small angular distances, you'll get $$N = 1$$, i.e. just one line segment, or two points, the starting and ending points. Its exact value depends on how precisely you want to calculate the curve length.

Finally, if you have the latitude ($$\theta$$) and longitude ($$\varphi$$) of some location on Earth, you can convert those to a point on the unit sphere using \left\lbrace ~ \begin{aligned} x &= \cos(\theta) \cos(\varphi) \\ y &= \cos(\theta) \sin(\varphi) \\ z &= \sin(\theta) \\ \end{aligned} \right. so that $$z$$ axis is North (with $$\theta \gt 0$$ on Northern hemisphere), $$x-y$$ plane is on the equator, $$x$$ axis towards zero longitude (zero meridian). If positive $$\varphi$$ go West and negative East from zero longitude (zero meridian), then positive $$y$$ axis is in the Americas, and negative $$y$$ axis in Asia.

It may feel tempting to use $$x = R_c \cos(\theta) \cos(\varphi)$$ and so on, but that will introduce angular error. Latitude and longitude on Earth are measured as if on an unit sphere. (The error is similar to drawing an ellipse with semimajor axis $$2$$ and semiminor axis $$1$$ using $$x = 2 \cos\phi$$, $$y = \sin\phi$$. The line through origin to point $$(2\cos 45°, \sin 45°)$$ does not make a $$45°$$ angle with the $$x$$ axis; the angle is only about $$26.565°$$. Of course, the error in Earth's case is smaller, as the difference in the radiuses is just over $$21\text{ km}$$, or about $$0.34\%$$.)

Hope this helps!

Make two rotations of the axes to align the the $$z$$-axis along the normal direction of the given plane $$(n_x,n_y,n_z)$$, which is to rotate $$\alpha$$-degrees around the $$z$$-axis, followed by $$\beta$$-degrees around the $$y$$-axis with $$\alpha$$ and $$\beta$$ given by

$$\cos\alpha = \frac{n_x}{\sqrt{n_x^2+n_y^2}},\>\>\>\sin\alpha = \frac{n_y}{\sqrt{n_x^2+n_y^2}}, \>\>\>\cos\beta= \frac{n_z}{\sqrt{n_x^2+n_y^2+n_z^2}},\>\>\>\sin\beta= \frac{\sqrt{n_x^2+n_y^2}}{\sqrt{n_x^2+n_y^2+n_z^2}}$$

The corresponding coordinate changes are

\begin{align} & x=(w\sin\beta+u\cos\beta)\cos\alpha - v\sin\alpha \\ & y=(w\sin\beta+u\cos\beta)\sin\alpha + v\cos\alpha \\ & z= w\cos\beta-u\sin\beta \end{align}

Then, substitute above into $$\frac{x^2+y^2}{R_e^2}+\frac{z^2}{R_p^2} =1$$ and set $$w=0$$ to obtain the equation of ellipse

$$\left(\frac{\cos^2\beta}{R_e^2}+\frac{\sin^2\beta}{R_p^2} \right)u^2+\frac{v^2}{R_e^2}=1$$

Thus, the axes of the ellipse are $$R_e$$ and $$\left(\frac{\cos^2\beta}{R_e^2}+\frac{\sin^2\beta}{R_p^2} \right)^{-\frac12}$$.