OpenGL: projecting mouse click onto geometry

glMatrixMode(GL_MODELVIEW); //Switch to the drawing perspective
glLoadIdentity(); //Reset the drawing perspective

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);

// replace the above gluUnProject call with
gluUnProject(sx, viewport[1] + viewport[3] - sy, 1, modelMatrix, projMatrix, viewport, &x, &y, &z);

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;