This question assumes an ellipsoidal model of the earth. Its reference surface is obtained by rotating an ellipse around its minor axis (plotted vertically by convention). Such an ellipse is just a circle that has been stretched horizontally by a factor of a and vertically by a factor of b. Using the standard parameterization of the unit circle,
t --> (cos(t), sin(t))
(which defines cosine and sine), we obtain a parameterization
t --> (a cos(t), b sin(t)).
(The two components of this parameterization describe a trip around the curve: they specify, in Cartesian coordinates, our location at "time" t.)
The geodetic latitude, f, of any point is the angle that "up" makes to the equatorial plane. When a differs from b, the value of f differs from that of t (except along the equator and at the poles).
In this picture, the blue curve is one quadrant of such an ellipse (greatly exaggerated compared to the earth's eccentricity). The red dot at the lower left corner is its center. The dashed line designates the radius to one point on the surface. Its "up" direction there is shown with a black segment: it is, by definition, perpendicular to the ellipse at that point. Due to the exaggerated eccentricity, it is easy to see that "up" is not parallel to the radius.
In our terminology, t is related to the angle made by the radius to the horizontal and f is the angle made by that black segment. (Note that any point on the surface can be viewed from this perspective. This allows us to limit both t and f to lie between 0 and 90 degrees; their cosines and sines will be positive, so we don't have to worry about negative square roots in the formulas.)
The trick is to convert from the t-parameterization to one in terms of f, because in terms of t the radius R is easy to compute (via the Pythagorean theorem). Its square is the sum of squares of the components of the point,
R(t)^2 = a^2 cos(t)^2 + b^2 sin(t)^2.
To make this conversion we need to relate the "up" direction f to the parameter t. This direction is perpendicular to the tangent of the ellipse. By definition, a tangent to a curve (expressed as a vector) is obtained by differentiating its parameterization:
Tangent(t) = d/dt (a cos(t), b sin(t)) = (-a sin(t), b cos(t)).
(Differentiation computes the rate of change. The rate of change of our position as we travel around the curve is, of course, our velocity, and that always points along the curve.)
Rotate this clockwise by 90 degrees to obtain the perpendicular, called the "normal" vector:
Normal(t) = (b cos(t), a sin(t)).
The slope of this normal vector, equal to (a sin(t)) / (b cos(t)) ("rise over run"), is also the tangent of the angle it makes to the horizontal, whence
tan(f) = (a sin(t)) / (b cos(t)).
Equivalently,
(b/a) tan(f) = sin(t) / cos(t) = tan(t).
(If you have good insight into Euclidean geometry, you could obtain this relationship directly from the definition of an ellipse without going through any trig or calculus, simply by recognizing that the combined horizontal and vertical expansions by a and b respectively have the effect of changing all slopes by this factor b / a.)
Look again at the formula for R(t)^2: we know a and b -- they determine the shape and size of the ellipse -- so we only need to find cos(t)^2 and sin(t)^2 in terms of f, which the preceding equation lets us do easily:
cos(t)^2 = 1/(1 + tan(t)^2)
= 1 / (1 + (b/a)^2 tan(f)^2)
= a^2 / (a^2 + b^2 tan(f)^2);
sin(t)^2 = 1 - cos(t)^2
= b^2 tan(f)^2 / (a^2 + b^2 tan(f)^2).
(When tan(f) is infinite, we're at the pole, so just set f = t in that case.)
This is the connection we need. Substitute these values for cos(t)^2 and sin(t)^2 into the expression for R(t)^2 and simplify to get
R(f)^2 = ( a^4 cos(f)^2 + b^4 sin(f)^2 ) / ( a^2 cos(f)^2 + b^2 sin(f)^2 ).
A simple transformation shows that this equation is the same as the one found on Wikipedia. Because a^2 b^2 = (ab)^2 and (a^2)^2 = a^4,
R(f)^2 = ( (a^2 cos(f))^2 + (b^2 sin(f))^2 ) / ( (a cos(f))^2 + (b sin(f))^2 )
Best Answer
One could use either kind of latitude to locate points on the WGS 84 ellipsoid (used by the NED) or any other ellipsoid, but "everybody knows" that the values will always be given as geodetic latitudes. However, it is surprisingly hard to find an authoritative statement to that effect!
Before we go on, it helps to understand that although a datum like the WGS 84 ellipsoid describes a reference surface to locate everything on, above, and within the Earth, it is compatible with many different ways to associate coordinates with the geometric points on that ellipsoid. Two of those ways are the geocentric and geodetic systems.
The geocentric latitude of a point is the angle it makes with the Equatorial plane.
The geodetic latitude of a point is the angle made between the Equatorial plane and a perpendicular from that point to the surface of the datum beneath it.
The two angles are the same only when the point is located over the Equator (when they equal zero) or over one of the poles (when they equal -90 or +90 degrees). Otherwise the points denoted by the same numerical values of latitude and longitude can differ by as much as several tens of kilometers in the two systems (as I recall--I haven't rechecked that figure recently).
Both types of latitude can in principle be used, even though WGS 84 is referred to as a "geodetic" datum. (Indeed, both of them, as well as a "reduced latitude," are used in the governing NIMA technical report TR8350.2, available at http://earth-info.nga.mil/GandG/publications/tr8350.2/wgs84fin.pdf. See p. 4-4 and section 5.2 (on Gravity Potential) et seq.)
In the recent article Using NHDPlus as the Land Base for the Noah-distributed Model (David, Maidment, et al., Transactions in GIS 2009; 13(4): 363-377) a group of GIS luminaries writes
(Emphasis added.)
In the 2007 book Unmanned Aerial Vehicle Real-time Guidance System Via State-space Heuristic Search, Manuel Soto makes a similarly clear statement about which latitude is used by the NED (and employs similar figures to illustrate the distinction). The relevant pages (85-88) are available on Google books.