You should use gluUnProject
:
First, compute the "unprojection" to the near plane:
GLdouble modelMatrix[16];
GLdouble projMatrix[16];
GLint viewport[4];
glGetIntegerv(GL_VIEWPORT, viewport);
glGetDoublev(GL_MODELVIEW_MATRIX, modelMatrix);
glGetDoublev(GL_PROJECTION_MATRIX, projMatrix);
GLdouble x, y, z;
gluUnProject(sx, viewport[1] + viewport[3] - sy, 0, modelMatrix, projMatrix, viewport, &x, &y, &z);
and then to the far plane:
// replace the above gluUnProject call with
gluUnProject(sx, viewport[1] + viewport[3] - sy, 1, modelMatrix, projMatrix, viewport, &x, &y, &z);
Now you've got a line in world coordinates that traces out all possible points you could have been clicking on. So now you just need to interpolate: suppose you're given the z-coordinate:
GLfloat nearv[3], farv[3]; // already computed as above
if(nearv[2] == farv[2]) // this means we have no solutions
return;
GLfloat t = (nearv[2] - z) / (nearv[2] - farv[2]);
// so here are the desired (x, y) coordinates
GLfloat x = nearv[0] + (farv[0] - nearv[0]) * t,
y = nearv[1] + (farv[1] - nearv[1]) * t;