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