Details such as mountains are not considered, we assume that the planet is locally flat. For calculating gravity, we assume that the planet's mass may be modelled by a point-mass at its centre. This model is exact for spherical objects, and practial for planets.

Consider a ball sitting somewhere on the surface of the planet. Two accelerations act upon the ball at that point: gravity , and normal reaction . The sum of these accelerations keeps the ball moving with the planet, in circular motion about the axis.

The normal reaction at a point on the surface is perpendicular to the surface.

We need only consider a cross-section of the planet, passing through the axis:

We know the formulae for the acceleration due to gravity, and the acceleration needed to achieve circular motion:

is the gravitational constant, is the mass of the planet, is the vector from the center of the planet (not used directly), is the magnitude of , is the unit vector in the direction of , is the rate at which the planet is spinning (radians per second), and is the perpendicular vector from the axis of the planet. The gravity vector is directed toward the centre of the planet. The circular motion vector is directed toward the nearest point on the axis of the planet.

We can express the normal reaction in these terms:

Consider a point on the surface , using cartesian coordinates with origin at the center of the planet, -axis corresponding to the planet's equator and -axis to the planet's axis, where and . We split the vectors into their and components, and consider the magnitude .

Now, considering the curve of the surface, it is perpendicular to the normal:

Let :

This equation determines the curve of the surface, and can be integrated as follows:

Substitute , , :

At the pole, and , so we can find :

So we have as a function of :

We can express in terms of , giving a parametric equation for the curve of the planet's surface, which we can plot using a computer:

At the equator, , so the equatorial radius satisfies:

This is a cubic in which can be solved exactly, but I have not done so here.

I wrote a computer program to draw planet shapes based on these formulae, and to calculate the amount of bulge at the equator. Here is the picture generated of Saturn, which has a visible bulge compared to the reference circle:

The values calculated for the bulges (equatorial diameter minus polar diameter) of Earth and Saturn are a bit different from what has been observed in real life:

Earth: computed bulge: 21.89 km, actual bulge: 42.72 km.

Saturn: computed bulge: 8239.7 km, actual bulge: 11808 km.

The calculated values are not close enough to the actual values, so I think I have made some mistake, or did not consider some significant factor. Nevertheless, I think my method is the right approach.

Here is a picture of the shape Saturn would have, if it were spinning 56% faster than it really does, at the critical point.

Here is a picture of the shape Saturn would have, if it were spinning 59% faster than it really does. The shape is no longer closed. Saturn could not spin that fast and remain in one piece.

Here is the program I used, it is written in brace, my dialect of C:

#!/usr/bin/env bx use b cstr usage[] = "m r0 p [steps=10000]", "mass polar-radius spin-period", "5.9736e24 6356800 86164.09 # Earth", "5.6869e26 54364000 37050.56 # Saturn", NULL num G = 6.67428e-11 # m^3/kg/s Main() num M, r0, p, s, a, b, R2, r2, r, x, y, yp, xp, R, dR, r1, steps steps = 10000 Getargs(num, M, r0, p) getargs(num, steps) warn("M=%g r0=%f p=%f", M, r0, p) s = 2*pi/p paper(600, 400) int max_wh = imax(w, h) zoom((max_wh * 0.6) / (r0*2)) blue() circle(0, 0, r0) black() a = 2*G*M/(s*s) b = 1/r0 x = 0, y = r0, yp = y, xp = 0 dR = r0 / steps for R=r0; R<r0*2; R+=dR: again R2 = R*R r2 = a*(b - 1/R) if r2 < -tinynum || r2 > R2 + tinynum: break r2 = clamp(r2, 0, R2) r = sqrt(r2) x = r y = sqrt(R2 - r2) if sd(x) > w_2: break point(x, y) point(-x, y) point(x, -y) point(-x, -y) Paint() if xp && hypot(yp-y, xp-x) > r0 / h / 2: R-=dR ; dR /= 2 ; R+=dR ; again eif xp && hypot(yp-y, xp-x) < r0 / h / 8: R-=dR ; dR *= 2 ; R+=dR ; again yp = y ; xp = x r1 = hypot(x, y) warn("r1 ~= %f", r1) warn("equatorial bulge = 2*(r1-r0) = %f km", 2*(r1-r0)/1000)