How do I reverse-project 2D points into 3D?

void YCamera :: CalculateWorldCoordinates(float x, float y, YVector3 *vec)
{
    //	START
    GLint viewport[4];
    GLdouble mvmatrix[16], projmatrix[16];

GLint real_y;
GLdouble mx, my, mz;

glGetIntegerv(GL_VIEWPORT, viewport);
glGetDoublev(GL_MODELVIEW_MATRIX, mvmatrix);
glGetDoublev(GL_PROJECTION_MATRIX, projmatrix);

real_y = viewport[3] - (GLint) y - 1;	// viewport[3] is height of window in pixels
gluUnProject((GLdouble) x, (GLdouble) real_y, 1.0, mvmatrix, projmatrix, viewport, &mx, &my, &mz);

/*	'mouse' is the point where mouse projection reaches FAR_PLANE.
	World coordinates is intersection of line(camera->mouse) with plane(z=0) (see LaMothe 306)

	Equation of line in 3D:
		(x-x0)/a = (y-y0)/b = (z-z0)/c

	Intersection of line with plane:
		z = 0
		x-x0 = a(z-z0)/c  <=> x = x0+a(0-z0)/c  <=> x = x0 -a*z0/c
		y = y0 - b*z0/c

*/
double lx = fPosition.x - mx;
double ly = fPosition.y - my;
double lz = fPosition.z - mz;
double sum = lx*lx + ly*ly + lz*lz;
double normal = sqrt(sum);
double z0_c = fPosition.z / (lz/normal);

vec->x = (float) (fPosition.x - (lx/normal)*z0_c);
vec->y = (float) (fPosition.y - (ly/normal)*z0_c);
vec->z = 0.0f;

GLint real_y;
GLdouble mx, my, mz;

glGetIntegerv(GL_VIEWPORT, viewport);
glGetDoublev(GL_MODELVIEW_MATRIX, mvmatrix);
glGetDoublev(GL_PROJECTION_MATRIX, projmatrix);

real_y = viewport[3] - (GLint) y - 1;	// viewport[3] is height of window in pixels
gluUnProject((GLdouble) x, (GLdouble) real_y, 1.0, mvmatrix, projmatrix, viewport, &mx, &my, &mz);

/*	'mouse' is the point where mouse projection reaches FAR_PLANE.
	World coordinates is intersection of line(camera->mouse) with plane(z=0) (see LaMothe 306)

	Equation of line in 3D:
		(x-x0)/a = (y-y0)/b = (z-z0)/c

	Intersection of line with plane:
		z = 0
		x-x0 = a(z-z0)/c  <=> x = x0+a(0-z0)/c  <=> x = x0 -a*z0/c
		y = y0 - b*z0/c

*/
double lx = fPosition.x - mx;
double ly = fPosition.y - my;
double lz = fPosition.z - mz;
double sum = lx*lx + ly*ly + lz*lz;
double normal = sqrt(sum);
double z0_c = fPosition.z / (lz/normal);

vec->x = (float) (fPosition.x - (lx/normal)*z0_c);
vec->y = (float) (fPosition.y - (ly/normal)*z0_c);
vec->z = 0.0f;