How to draw a sphere. Suppose we want to draw something like a planet. This page describes a developmental story, with the following way-points
How can we draw a sphere - without glutSolidSphere? One way is like a geographer does it, with lines of longitude and latitude:

We can pick out a point on the sphere with an angle φ (Greek letter phi ) around from the North pole, and an angle θ ( Greek theta ) around the equator.
If we change phi and theta a bit, we get this:

That orange bit will be pretty flat if the changes in theta and phi are small - so we can draw it as an OpenGL quad. That is the basic idea - we just need some trigonometry.
Here is the start of our Planet class:
public class Planet {
/** location of the centre */
private double x, y, z;
/** radius */
private double radius;
TGABuffer imageBuffer;
int texture;
public Planet(double x, double y, double z, double radius) {
this.x = x;
this.y = y;
this.z = z;
this.radius = radius;
}
We set up the texture like this:
/**
* Set up textures for the planet
* This is called by the init() of the main displaying class. It
* loads and sets up required textures for all instances of this class
*
* @param gl The gl of the displaying window
*/
public void setupTexture(GL gl) {
final int[] tmp = new int[1];
gl.glGenTextures(1, tmp, 0);
String filename = "C:/..;
imageBuffer = TGABufferMaker.make(filename);
gl.glPixelStorei(GL.GL_UNPACK_ALIGNMENT, 1);
texture = tmp[0];
gl.glBindTexture(GL.GL_TEXTURE_2D, texture);
gl.glTexParameteri(GL.GL_TEXTURE_2D, GL.GL_TEXTURE_WRAP_S, GL.GL_REPEAT);
gl.glTexParameteri(GL.GL_TEXTURE_2D, GL.GL_TEXTURE_WRAP_T, GL.GL_REPEAT);
gl.glTexParameteri(GL.GL_TEXTURE_2D, GL.GL_TEXTURE_MAG_FILTER,
GL.GL_LINEAR);
gl.glTexParameteri(GL.GL_TEXTURE_2D, GL.GL_TEXTURE_MIN_FILTER,
GL.GL_LINEAR);
// connect to data
gl.glTexImage2D(GL.GL_TEXTURE_2D, 0, GL.GL_RGBA, imageBuffer.getWidth(),
imageBuffer.getHeight(), 0, GL.GL_RGBA, GL.GL_UNSIGNED_BYTE,
imageBuffer.getBuffer());
}
Then we draw it like this:
public void draw(GL gl) {
double theta, phi;
double twoPi = Math.PI * 2;
final int divs = 64;
double dTheta = twoPi / divs;
double dPhi = Math.PI / divs;
double x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4;
/** want texture */
gl.glEnable(GL.GL_TEXTURE_2D);
gl.glTexEnvf(GL.GL_TEXTURE_ENV, GL.GL_TEXTURE_ENV_MODE, GL.GL_MODULATE);
/** texturePos is the texture coord position of the start of the disc */
float texPosx = 0, texPosy = 0;
/** TextureStep is tex coord size to match angular steps */
float texStep = 1.0f / divs;
gl.glBindTexture(GL.GL_TEXTURE_2D, texture);
for (theta = 0; theta < twoPi; theta += dTheta) {
texPosy = 0;
for (phi = 0; phi < Math.PI; phi += dPhi) // phi is the latitude, as angle from North pole axis
{
y1 = radius * Math.cos(phi);
x1 = radius * Math.sin(phi) * Math.cos(theta);
z1 = radius * Math.sin(phi) * Math.sin(theta);
y2 = radius * Math.cos(phi);
x2 = radius * Math.sin(phi) * Math.cos(theta + dTheta);
z2 = radius * Math.sin(phi) * Math.sin(theta + dTheta);
y3 = radius * Math.cos(phi + dPhi);
x3 = radius * Math.sin(phi + dPhi) * Math.cos(theta + dTheta);
z3 = radius * Math.sin(phi + dPhi) * Math.sin(theta + dTheta);
y4 = radius * Math.cos(phi + dPhi);
x4 = radius * Math.sin(phi + dPhi) * Math.cos(theta);
z4 = radius * Math.sin(phi + dPhi) * Math.sin(theta);
gl.glBegin(GL.GL_QUADS);
gl.glNormal3dv(WMVector.normal(x2, y2, z2, x1, y1, z1, x4, y4, z4), 0);
gl.glTexCoord2f(texPosx, texPosy);
gl.glVertex3d(x1, y1, z1);
gl.glTexCoord2f(texPosx + texStep, texPosy);
gl.glNormal3dv(WMVector.normal(x3, y3, z3, x2, y2, z2, x1, y1, z1), 0);
gl.glVertex3d(x2, y2, z2);
gl.glNormal3dv(WMVector.normal(x4, y4, z4, x3, y3, z3, x2, y2, z2), 0);
gl.glTexCoord2f(texPosx + texStep, texPosy + texStep);
gl.glVertex3d(x3, y3, z3);
gl.glNormal3dv(WMVector.normal(x1, y1, z1, x4, y4, z4, x3, y3, z3), 0);
gl.glTexCoord2f(texPosx, texPosy + texStep);
gl.glVertex3d(x4, y4, z4);
gl.glEnd();
texPosy += texStep;
}
texPosx += texStep;
}
}
Which looks like 
We have a nested loop changing theta and phi, and in the middle we calculate and draw a quad with these vertices:

The texture image is a Mercator Projection photo of Jupiter. When flat it looks weird, but when you wrap it around a sphere like this it is perfect.
So we have stage one - we've drawn a sphere. But it does not look very good, and it is slow - about 25 frames per second.
One reason it is slow is because we calculate each point four times. For example, the 'x4' point, the next time around on the phi loop, becomes the 'x1' point - but we calculate it again even though we've done so already. There is only one new point for each quad, so we should be able to do it it with just one vertex calculation each time.
To do this, we must
1. First fill an array with calculated vertex locations down the first longitude.
2. Then down the loop, we use those values plus one new vertex each time. We draw those, and put the new vertex in the array for the next time around.
We use a Triple class for the array, like a C struct:
class Triple {
double x, y, z;
public Triple(double x, double y, double z) {
this.x = x;
this.y = y;
this.z = z;
}
}
/**
* Draw a sphere version 2.
* We store the vertices down a line of longitude, so that
* we only calculate each point once
*
*
*/
public void draw(GL gl) {
double theta, phi;
double twoPi = Math.PI * 2;
final int divs = 64;
double dTheta = twoPi / divs;
double dPhi = Math.PI / divs;
double x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4;
// the array..
Triple[] vertexData = new Triple[divs + 1];
for (int i = 0; i < divs; i++) {
vertexData[i] = new Triple(0, 0, 0);
}
/** want texture */
gl.glEnable(GL.GL_TEXTURE_2D);
gl.glTexEnvf(GL.GL_TEXTURE_ENV, GL.GL_TEXTURE_ENV_MODE, GL.GL_MODULATE);
/** texturePos is the texture coord position of the start of the disc */
float texPosx = 0, texPosy = 0;
/** TextureStep is tex coord size to match angular steps */
float texStep = 1.0f / divs;
gl.glBindTexture(GL.GL_TEXTURE_2D, texture);
// calculate the vertices down the Greenwich meridian, where theta=0;
int arrayPos;
phi = 0;
for (arrayPos = 0; arrayPos <
divs + 1; arrayPos++) {
y1 = radius * Math.cos(phi);
x1 = radius * Math.sin(phi);
z1 = 0;
phi += dPhi;
vertexData[arrayPos] = new Triple(x, y, z);
}
for (theta = 0; theta < twoPi; theta += dTheta) {
// first 2 - at N Pole
y2 = radius;
x2 = 0;
z2 = 0;
texPosy = 0;
phi = 0;
for (arrayPos = 0; arrayPos < divs; arrayPos++) // phi is the latitude, as angle from North pole axis
{
y1 = vertexData[arrayPos].y;
x1 = vertexData[arrayPos].x;
z1 = vertexData[arrayPos].z;
y3 = radius * Math.cos(phi + dPhi);
x3 = radius * Math.sin(phi + dPhi) * Math.cos(theta + dTheta);
z3 = radius * Math.sin(phi + dPhi) * Math.sin(theta + dTheta);
y4 = vertexData[arrayPos + 1].y;
x4 = vertexData[arrayPos + 1].x;
z4 = vertexData[arrayPos + 1].z;
gl.glBegin(GL.GL_QUADS);
gl.glNormal3dv(WMVector.normal(x2, y2, z2, x1, y1, z1, x4, y4, z4), 0);
gl.glTexCoord2f(texPosx, texPosy);
gl.glVertex3d(x1, y1, z1);
gl.glTexCoord2f(texPosx + texStep, texPosy);
gl.glNormal3dv(WMVector.normal(x3, y3, z3, x2, y2, z2, x1, y1, z1), 0);
gl.glVertex3d(x2, y2, z2);
gl.glNormal3dv(WMVector.normal(x4, y4, z4, x3, y3, z3, x2, y2, z2), 0);
gl.glTexCoord2f(texPosx + texStep, texPosy + texStep);
gl.glVertex3d(x3, y3, z3);
gl.glNormal3dv(WMVector.normal(x1, y1, z1, x4, y4, z4, x3, y3, z3), 0);
gl.glTexCoord2f(texPosx, texPosy + texStep);
gl.glVertex3d(x4, y4, z4);
gl.glEnd();
texPosy += texStep;
vertexData[arrayPos] = new Triple(x2, y2, z2);
x2 = x3;
y2 = y3;
z2 = z3;
phi += dPhi;
}
texPosx += texStep;
}
}
This looks just like the first version, since the drawing is identical. But at the cost of the storage of the array - and worse, the pain of actually coding this - we look for an improvement in speed.
In fact for this the fps is just 38. We've reduced the calculation by a factor of 4, but the speed improvement is only about 50%. This means the time is being taken by the drawing, not the calculations.
So the next improvement is to use a display list. We need an ID for it, and in the setUpTexture function, called from init(), we add:
// display list
displayList = gl.glGenLists(1);
gl.glNewList(displayList, GL.GL_COMPILE);
draw2(gl);
gl.glEndList();
where draw2 is just our last draw, renamed, and we have a new draw which is simply:
public void draw(GL gl) {
gl.glCallList(displayList);
}
Again this looks the same, because the drawing is unchanged, but we get over 100 fps.
So we've got a fast planet - but it still looks a bit rubbish. Suppose we
close into Jupiter and take a closer look at those quads: 
If you look at the edge of the planet, it is pretty curvy. Our 64 divisions looks fine in terms of the shape. The problem with it is the lighting on each quad. And that is because we (I'm including you in this) are working out the normals using just the four vertices of each quad. We actually need another 4 quads to do this properly:

Aargh!! How are we going to do this, without calculating everything maybe 8 times?
Again we can do it better if we use more storage. Instead of just holding the vertex data down one longitude line, we use a 2 D array, and calculate all vertices. Then we draw them all, like this:
public void draw3(GL gl) {
double theta, phi;
double twoPi = Math.PI * 2;
final int divs = 32;
double dTheta = twoPi / divs;
double dPhi = Math.PI / divs;
double x, y, z, x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4;
double x5, y5, z5, x6, y6, z6, x7, y7, z7, x8, y8, z8;
// the array..
Triple[][] vertexData = new Triple[divs + 1][divs + 1];
theta = 0;
for (int i = 0; i < divs+1; i++) {
phi = 0;
for (int j = 0; j < divs+1; j++) {
y = radius * Math.cos(phi);
x = radius * Math.sin(phi) * Math.cos(theta);
z = radius * Math.sin(phi) * Math.sin(theta);
vertexData[i][j] = new Triple(x, y, z);
phi += dPhi;
}
theta += dTheta;
}
/** want texture */
gl.glEnable(GL.GL_TEXTURE_2D);
gl.glTexEnvf(GL.GL_TEXTURE_ENV, GL.GL_TEXTURE_ENV_MODE, GL.GL_MODULATE);
/** texturePos is the texture coord position of the start of the disc */
float texPosx = 0, texPosy = 0;
/** TextureStep is tex coord size to match angular steps */
float texStep = 1.0f / divs;
gl.glBindTexture(GL.GL_TEXTURE_2D, texture);
for (int i = 0; i < divs; i++) {
texPosy = 0;
for (int j = 0; j < divs-1; j++) // phi is the latitude, as angle from North pole axis
{
y1 = vertexData[i][j].y;
x1 = vertexData[i][j].x;
z1 = vertexData[i][j].z;
y2 = vertexData[i+1][j].y;
x2 = vertexData[i+1][j].x;
z2 = vertexData[i+1][j].z;
y3 = vertexData[i+1][j+1].y;
x3 = vertexData[i+1][j+1].x;
z3 = vertexData[i+1][j+1].z;
y4 = vertexData[i][j+1].y;
x4 = vertexData[i][j+1].x;
z4 = vertexData[i][j+1].z;
y5 = vertexData[(i+2)%(divs+1)][j].y;
x5 = vertexData[(i+2)%(divs+1)][j].x;
z5 = vertexData[(i+2)%(divs+1)][j].z;
y6 = vertexData[(i+2)%(divs+1)][j+1].y;
x6 = vertexData[(i+2)%(divs+1)][j+1].x;
z6 = vertexData[(i+2)%(divs+1)][j+1].z;
y7 = vertexData[(i+2)%(divs+1)][j+2].y;
x7 = vertexData[(i+2)%(divs+1)][j+2].x;
z7 = vertexData[(i+2)%(divs+1)][j+2].z;
y8 = vertexData[i][j+2].y;
x8 = vertexData[i][j+2].x;
z8 = vertexData[i][j+2].z;
gl.glBegin(GL.GL_QUADS);
gl.glNormal3dv(WMVector.normal(x2, y2, z2, x1, y1, z1, x4, y4, z4), 0);
gl.glTexCoord2f(texPosx, texPosy);
gl.glVertex3d(x1, y1, z1);
gl.glTexCoord2f(texPosx + texStep, texPosy);
gl.glNormal3dv(WMVector.normal(x5, y5, z5, x2, y2, z2, x3, y3, z3), 0);
gl.glVertex3d(x2, y2, z2);
gl.glNormal3dv(WMVector.normal(x6, y6, z6, x3, y3, z3, x7, y7, z7), 0);
gl.glTexCoord2f(texPosx + texStep, texPosy + texStep);
gl.glVertex3d(x3, y3, z3);
gl.glNormal3dv(WMVector.normal(x3, y3, z3, x4, y4, z4, x8, y8, z8), 0);
gl.glTexCoord2f(texPosx, texPosy + texStep);
gl.glVertex3d(x4, y4, z4);
gl.glEnd();
texPosy += texStep;
}
texPosx += texStep;
}
}
We write this into a display list, and it looks like:
This looks much better, even though the divs have been reduced to 32, and the fps is around 150.