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 )
