LinePlaneIntersection LinePlaneIntersectionFunc(Point line_p1, Point line_p2,
			    Point plane_p1, Point plane_p2, Point plane_p3)
{
  Matrix m = { { line_p1.x-line_p2.x, plane_p2.x-plane_p1.x, plane_p3.x-plane_p1.x, 0.0,
	       line_p1.y-line_p2.y, plane_p2.y-plane_p1.y, plane_p3.y-plane_p1.y, 0.0,
	       line_p1.z-line_p2.z, plane_p2.z-plane_p1.z, plane_p3.z-plane_p1.z, 0.0,
	       0.0,                 0.0,                   0.0,                   1.0 } };
  
  //std::cout << "LinePlaneIntersectionFunc" << m << std::endl;
  Matrix mi = Matrix::Inverse(m);
  Vector v;
  v.dx = line_p1.x - plane_p1.x;
  v.dy = line_p1.y - plane_p1.y;
  v.dz = line_p1.z - plane_p1.z;
  Vector tuv = v*mi;
  //tuv.Normalize();
  LinePlaneIntersection sec;
  //tuv.dx = 0.5; 
  sec.tuv = tuv;
  //std::cout << tuv.dx << std::endl;
  return sec;

}


  virtual Point Pos(int i) const
  {
    Point p = pts->Pos(i);
    Vector v = light_vec;
    LinePlaneIntersection sect = LinePlaneIntersectionFunc(p, p+v,
							   pos, pos+u_x, pos+u_y);
    return p+v*sect.tuv.dx;
  }