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
Best Answer
Let the helix be given by $(\cos t, \sin t, ht)$ (after scaling). If $P$ is your point $(a,b,c)$, and $Q = (\cos t, \sin t, ht)$ is the nearest point on the helix, then $PQ$ is perpendicular to the tangent at $Q$, which is just $(-\sin t, \cos t, h)$:
$-(\cos t - a)\sin t + (\sin t - b)\cos t + (ht - c)h = 0 $
This simplifies to $A \sin(t+B) + Ct + D = 0$ for some constants $A,B,C,D,$ as Moron said. But then you have to solve this numerically. There will be more than one solution in general, but (as Jonas Kibelbek pointed out in the comments) you only need to check the solutions with $z$-coordinate in the interval $[c-\pi h, c+\pi h)$.