I declared one point in each octant of the globe, transformed it to XYZ using the equations I had, and then tested martin f's answer. It didn't returned the same points. Then I delved deeper into Wikipedia's equations, and I finally understood how they worked. Then I adapted them.
latitude(r, x, y, z) = arcsin(z/r)(180/π)
longitude(r, x, y, z) =
if (x > 0) {
arctan(y/x)(180/π)
} else if (y > 0) {
arctan(y/x)(180/π) + 180
} else {
arctan(y/x)(180/π) - 180
}
(I wanted to write a piecewise function, but LaTeX isn't supported here) This is considering North and East as positive and right hand rule for XYZ, where Y is transverse to the plane that contains the Greenwhich meridian, positive to the East as well, and Z is aligned to North-South.
Thanks for helping me get there.
Edit: I just checked on Wikipedia, and it seems I learned the wrong right hand rule. Here, I mean X in the thumb, Y in the index finger and Z in the middle finger.
Latitude and Longitude are defined upon an ellipsoid.
Essentially you can have different lon/lat. NAD 27 lon/lat and WGS84 lon/lat will have subtle differences. As a rule, when people nowadays refer to lon/lat the mean WGS84 datum and WGS84 ellipsoid.
(The example is extracted from the eye opening book of PostGIS In Action from Regina O. Obe and Leo S. Hsu.)
So back to your question and numbering of thigs:
You start by modeling the earth using some variant of a reference ellipsoid, which should be the ellipsoid that deviates least from the geoid for the regions on earth you care about. The most common ellipsoids are GRS80 (specially used in North America) and WGS84 (used by GPS systems).
You use a datum to pin the ellipsoid to an actual place on earth and you assign a coordinate reference system to the ellipsoid so you can identify every point on the surface.
Best Answer
I have done this before using forward azimuths, here is a link that has descriptions and algorithms that may be helpful:
Inverse/Forward Utilities
A forward azimuth calculates a new point that is a specified distance and compass bearing from a starting point. The basic idea is that you have a point in Lat/Lon and you calculate a series of forward azimuths with a range of angles. The number of calculations depends on how many vertices you want in your circle, normally I would use 360 so I get one vertex for every degree of the circle which should look good for most applications.