Today I received this letter.

Dear Roine,
I’m trying to do ellipses in OpenGL. The ellipses turn out fine, but I can’t seem to position them correctly. All my tries have turned into giant hairballs. Please advice.
Confused

Dear Confused. I know just what you are going through. Just the other day (well, today actually) I needed to do the same thing. I started out with a bunch of glRotatef(..) to position the ellipses with angles computed from the vectors with arcsine, but the cyclic issues and rounding issues with e.g. sin(asin(x)) made the ellipses go bonkers. Google didn’t turn up more than hints, so I had to write this beautiful piece of code myself:

// Draw an ellipse with focus points r1 and r2, with radius r
// 'quad' is a previously allocated GLUquadric object
// The code avoids using inverse trigonometric functions and the cyclic issues with them.
void
ellipse(float r1x, float r1y, float r1z, float r2x, float r2y, float r2z, float r, GLUquadric* quad )
{
  glPushMatrix();

  // As with all OpenGL code, this is done 'backwards', the first transformation being applied last:

  // Finally, create a coordinate system transformation for the ellipse,
  // where the X axis transforms to the desired vector (i.e. r2-r1).
  // 'by' and 'bz' are created by cross products, which is perfectly
  // valid since an ellipse is rotationally symmetric about (in this case)
  // the X-axis and so any orthonormal vectors to 'r2-r1' are valid choices.

  const float d = sqrt((r2x-r1x)*(r2x-r1x)+(r2y-r1y)*(r2y-r1y)+(r2z-r1z)*(r2z-r1z));  // d = |r2-r1|

  const float bx[] = { (r2x-r1x)/d, (r2y-r1y)/d, (r2z-r1z)/d };  // normalized r2-r1

  float by[3];
  const float dx2 = bx[2]*bx[2]+bx[1]*bx[1];          // |[ex x bx]|^2
  // Unfotunately, there is no vector guaranteed to be orthogonal to r2-r1;
  // we have to test for it.
  // If bx happens to be (nearly) parallel to ex then use ey instead to avoid degenerate case.
  if (dx2 > 0.01) {
    // normalized [ex x bx]
    const float dx = sqrt(dx2);                       // |[ex x bx]|
    by[0] = 0;
    by[1] = -bx[2]/dx;
    by[2] = bx[1]/dx;
  } else {
    // normalized [ey x bx]
    const float dy = sqrt(bx[0]*bx[0]+bx[2]*bx[2]);   // |[ey x bx]|
    by[0] = bx[2]/dy;
    by[1] = 0;
    by[2] = -bx[0]/dy;
  }

  // Now make one last orthonormal vector to both bx and by crossing them
  const float bz[] = { bx[1]*by[2]-bx[2]*by[1],
                       bx[2]*by[0]-bx[0]*by[2],
                       bx[0]*by[1]-bx[1]*by[0] };  // [bx x by]

  // This matrix both rotates the ellipse to the desired vector r2-r1
  // and translates it to r1 (last column)
  // 'm' is column major as required by OpenGL.
  const float m[16] = { bx[0], bx[1], bx[2], 0.0,
                        by[0], by[1], by[2], 0.0,
                        bz[0], bz[1], bz[2], 0.0,
                        r1x, r1y, r1z, 1.0};
  glMultMatrixf(m);

  // Reposition ellipse with left focus in origo, and right focus at 'd' on X axis
  glTranslatef(d/2.0,0.0,0.0);

  // Stretch sphere to an ellipse with semimajor radius 'r' and distance 'd' between foci
  glScalef((d+r)/2.0,r,r);    

  // Rotate sphere so the rotationally symmetric axis is in X direction (looks better)
  glRotatef(90,0.0f,1.0f,0.0f);  

  // Create unit sphere
  gluSphere(quad,1.0,12,8);    

  glPopMatrix();
}