This can be accomplished with Thaddeus Vincenty's inverse solution or with the haversine distance formula. If simplicity and speed of calculation is more important than accuracy, then use the haversine distance formula. Otherwise, go with Vincenty's inverse solution.
Vincenty's Inverse Solution
$\alpha$ length of the semi-major axis of the ellipsoid
$\beta$ length of the semi-minor axis of the ellipsoid
$\gamma=\frac{1}{\alpha}(\alpha-\beta)$ flattening of the ellipsoid
$x_1, x_2$ latitude of the points
in radians
$y_1, y_2$ longitude of the points
in radians
$\psi=y_2-y_1$ difference in longitude
$\lambda=\psi$ first and current approximation
$\lambda_0$ previous approximation
Below are some trigonometric optimizations, for $k=1,2$
\[ \tan\omega_k= (1-\gamma)\cdot\tan x_k \]
\[ \cos\omega_k= \frac{1}{\sqrt{1+\tan^2\omega_k}} \]
\[ \sin\omega_k=\tan\omega_k\cdot\cos\omega_k \]
Now we iterate the following calculations until $\lambda-\lambda_0 > 10^{-12}$mm
\[ \sin\phi=\sqrt{(\cos\omega_2\cdot\sin\lambda)^2+(\cos\omega_1\cdot\sin\omega_2-\sin\omega_1\cdot\cos\omega_2\cdot\cos\lambda)^2} \]
\[ \cos\phi=\sin\omega_1\cdot\sin\omega_2+\cos\omega_1\cdot\cos\omega_2\cdot\cos\lambda \]
\[ \phi=\arctan\left(\frac{\sin\phi}{\cos\phi}\right) \]
\[ \sin z = \frac{\cos\omega_1\cdot\cos\omega_2\cdot\sin\lambda}{\sin\phi} \]
\[ \cos^2 z = 1-\sin^2 z \]
\[ \cos 2\phi_m = \cos\phi-\frac{2\sin\omega_1\cdot\sin\omega_2}{\cos^2 z} \]
\[ \delta = \frac{\gamma}{16}\cos^2 z\cdot(4+\gamma\cdot(4-3\cos^2 z)) \]
\[ \lambda_0 = \lambda \]
\[ \zeta=\phi+\delta\cdot\sin\phi\cdot(\cos 2\phi_m+\delta\cdot\cos\phi\cdot(-1+2\cos^2 2\phi_m)) \]
\[ \lambda = \psi+(1-\delta)\cdot\gamma\cdot\zeta\cdot\sin z \]
Once $\lambda$ converges, calculate the following
\[ \mu^2=\frac{\cos^2 z\cdot(\alpha^2-\beta^2)}{\beta^2} \]
\[ A=1+\frac{\mu^2}{16384}\left(4096+\mu^2(-768+\mu^2(320-175\mu^2))\right) \]
\[ B=\frac{\mu^2}{1024}\left(256+\mu^2(-128+\mu^2(74-47\mu^2))\right) \]
\[ C=\frac{B}{6}\cos 2\phi_m(4\sin^2\phi-3)(4\cos^2 2\phi_m-3)\]
\[ D=\cos 2\phi_m+\frac{B}{4}\left(\cos\phi(2\cos^2 2\phi_m-1)-C\right) \]
\[ \Delta\phi=B\cdot D\cdot\sin\phi \]
\[ d = A\cdot \beta\cdot(\phi-\Delta\phi) \]
Finally we now have $d$, which is the ellipsoidal distance between $(x_1, y_1)$ and $(x_2, y_2)$ in meters. To convert distance $d$ to miles, just multiply $d$ by $0.000621371$. What this algorithm lacks in speed and simplicity, is made up by its accuracy of 0.5mm!
Haversine Distance Formula
$x_1, x_2$ latitude of the points
in radians
$y_1, y_2$ longitude of the points
in radians
$R$ radius of the earth in meters
\[ \alpha = \sin^2\left(\frac{x_2-x_1}{2}\right)+\cos x_1\cdot\cos x_2\cdot\sin^2\left(\frac{y_2-y_1}{2}\right) \]
\[ \beta = 2\cdot{\rm atan2}(\sqrt{\alpha}, \sqrt{1-\alpha}) \]
\[ d = R\cdot\beta \]
Now we easily have $d$, the ellipsoidal distance between $(x_1, x_2)$ and $(y_1, y_2)$ in meters. Again to convert distance $d$ to miles, just multiply $d$ by $0.000621371$.