The simplest general-purpose mathematical model would be to use spherical coordinates with standard latitude $\phi$ (north of equator) and longitude $\theta$ (say east of Greenwich):
$$
X(\theta,\phi)
=\left[\matrix{x\\y\\z\\}\right]
=\left[
\matrix{
r\cos\phi\cos\theta\\
r\cos\phi\sin\theta\\
r\sin\phi
}\right]
=r\left[
\matrix{
\cos\phi\cos\theta\\
\cos\phi\sin\theta\\
\sin\phi
}\right]
$$
From an arbitrary point $X$ as above, you can compute the east and north directions on the tangent plane from the partial derivatives $X_\theta$
and $X_\phi$, respectively (useful for local comparisons with compatible map projections):
$$
X_\theta
=\frac{\partial}{\partial\theta}X(\theta,\phi)
=r\left[
\matrix{
-\cos\phi\sin\theta\\
\cos\phi\cos\theta\\
0
}\right]
\qquad\text{(east)}
$$
$$
X_\phi
=\frac{\partial}{\partial\phi}X(\theta,\phi)
=r\left[
\matrix{
-\sin\phi\cos\theta\\
-\sin\phi\sin\theta\\
\cos\phi
}\right]
\qquad\text{(north)}
$$
The zenith (outward radial vector direction) is of course just
$$
X_r
=\frac{\partial}{\partial r}X(\theta,\phi)
=\left[
\matrix{
-\sin\phi\cos\theta\\
-\sin\phi\sin\theta\\
\cos\phi
}\right]
\qquad\text{(zenith).}
$$
Now given two points $X_1(\theta_1,\phi_1)$ and $X_2(\theta_2,\phi_2)$,
the great arc between them is on a circle intersecting
the (assumed) spherical surface of the earth
(with assumed constant radius $r$, which is a simplification)
and the plane through the center of the earth
(the origin of our coordinate system) and the two points $X_1$ & $X_2$.
The normal to this plane can therefore be found by taking the cross product:
$$
N = X_1 \times X_2
$$
And the distance along this great arc is just the arc length from the angle $\alpha$ between the radial vectors of the two points, which is given by the arc cosine of their dot product (a caveat, however, is that it is ill-conditioned for angles near $0$ or $\pi\text{ rad}=180^\circ$):
$$s=r\alpha\qquad\text{for}\qquad\cos\alpha=X_1 \cdot X_2$$
If you want the bearing (starting direction on the map) from $X_1$ to travel on this great arc path, well that, too, can be easily found. Take the tangent vector $V$ at $X_1$ in the direction of $X_2$ by crossing the tangent plane's normal, $N$, with $X_1$ (possibly up to sign):
$$V = N \times X_1 = (X_1 \times X_2) \times X_1$$
Plugging in the above definitions or using some vector and trigonometric identities will give you explicit formulas for each of these in terms of $\theta_i$ and $\phi_i$.
So now suppose you have the point $X_0$ and the line (or great arc) through points $X_1$ and $X_2$ on the surface of the sphere. To find the distance from $X_0$ to this great circle, we essentially need to project (or better, rotate) $X_0$ onto the plane of the great arc (or onto the great arc itself). First we need to find the closest point, $F$ (for 'foot'), to $X_0$ on that arc. This will be the normalized version of
$$
M = N \times \left(X_0 \times N\right)
$$
scaled by the radius $r$,
$$
F = r\,\widehat{M} = r\,\frac{M}{||M||}
$$
and thus our distance $d=r\beta$, from $X_0$ to the "line" through $X_1,X_2$,
can be obtained from the great arc angle $\beta$ from $X_0$ to $F$,
by taking the arc (inverse) cosine of
$$
\cos\beta=X_0\cdot\widehat{M}=\frac{X_0\cdot M}{||M||}=X_0\cdot\frac{F}{r}\,.
$$
Another explanation can be found at this stackoverflow post.
You can find more background on this, with good diagrams, under wikipedia's spherical coordinates, geometry, trigonometry, law of cosines, half-side formula and haversine formula pages, along with some pages to get a sense for the physical realism of this model such as the earth's radius, among others.
Playing around in sage, I got a distance of $4841.93165384$ feet, or $1613.97721795$ yards, or $1475.72976066$ meters:
var('x,y,z,lat,lon')
rkm = 6371 # km
rmi = 3959 # mi
d2r = pi / 180
ll = Matrix(RDF, 3, 2, [\
28.503946, -99.453589,\
28.485044, -99.453709,\
28.498230, -99.468340])
xyzll = lambda lat,lon: [\
(cos(d2r*lat)*cos(d2r*lon)).n(),\
(cos(d2r*lat)*sin(d2r*lon)).n(),\
(sin(d2r*lat)).n()]
xyz = Matrix(RDF, 3, 3)
for i in range(3):
xyz[i] = xyzll(ll[i,0],ll[i,1])
print xyz
v0 = vector(xyz[0])
v1 = vector(xyz[1])
v2 = vector(xyz[2])
vn = v1.cross_product(v2); vn=vn/vn.norm(); print 'vn:', vn, vn.norm()
vm = vn.cross_product(v0.cross_product(vn)); vm=vm/vm.norm(); print 'vm:', vm, vm.norm()
a0 = arccos(vm.dot_product(v0)); print 'a0:', a0, (a0/d2r).n()
print 'distance:', a0*rmi, a0*rkm
# The OUTPUT:
[-0.144339114164 -0.86684945361 0.477219283871]
[-0.144366780752 -0.867004401598 0.476929345107]
[-0.144570113877 -0.866859220201 0.477131611326]
vn: (-0.760935788228, -0.21083794795, -0.613615584791) 1.0
vm: (-0.144515375391, -0.866898313757, 0.477077163445) 1.0
a0: 0.000231632359231 0.0132715565826116
distance: 0.917032510197 1.47572976066
Here is another way. Choose $t \in [0,1]$, then let $p(t) = (x_1+ t (x_2-x_1), y_1+t(y_2-y_1))$.
Then $p(0) = (x_1,y_1)$, $p(1) = (x_2,y_2)$ and for $t \in (0,1)$, $p(t)$ will be a point in between.
This scheme works even if the two $x$ coordinates are the same.
Best Answer
Let's assume that points $A$ and $B$ are $(x_1,y_1)$ and $(x_2,y_2)$ in Cartesian coordinates and $(r_1,\theta_1)$ and $(r_2,\theta_2)$ in polar coordinates.
Then the direction angle of the line segment is given by
$$\theta_0=\begin{cases} \tan^{-1}\left(\frac{y2-y1}{x2-x1}\right), & x_1\ne x_2 \\ \frac{\pi}2, & x_1=x_2 \end{cases} $$ And the signed distance from the origin to the line containing points $A$ and $B$ is
$$d_0=r_1\sin(\theta_1-\theta_0)=r_2\sin(\theta_2-\theta_0)$$
(Use either formula for $d_0$: you get the same result.) If $d_0\ne 0$ then the line containing points $A$ and $B$ does not pass through the origin, and the polar equation of that line is
and to get only the segment on the line just require that $\theta$ is between $\theta_1$ and $\theta_2$.
If $d_0=0$ then the line containing points $A$ and $B$ does pass through the origin, and thus $r$ is not a function of $\theta$ and another kind of polar equation will be required to specify the line, such as $\theta=\theta_1$ with $r$ limited to be between $r_1$ and $r_2$.
Note that the "distance" $d_0$ may be negative. If that is not desired, replace $d_0$ with $-d_0$ and $\theta_0$ with $\theta_0+\pi$ or $\theta_0-\pi$.
This diagram explains much of those formulas.