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