public static Point3D PointOfIntersection(Point3D La, Point3D Lb, Point3D P0, Point3D P1, Point3D P2)
        {
            //thanks to Paul Bourke and Bryan Hanson 
            //la,Lb constitute the line. P0,P1,P2 constitute the plane
            Point3D N = CrossProduct(P0, P1, P2);  //N is the normal of the plane        
            double n = DotProduct(N, new(P2.X - La.X, P2.Y - La.Y, P2.Z - La.Z));
            Point3D LbMinusLa = new Point3D(Lb.X - La.X, Lb.Y - La.Y, Lb.Z - La.Z);
            double d = DotProduct(N, LbMinusLa);
            if (d == 0) return null;//Line is parallel to or in the plane
            double u = n / d;
            if ((u >= 0.0) && (u <= 1.0))
            {//Plane is between the two points
            }//can be used for checking but does not influence the outcome
            Point3D IntersectingPoint = new Point3D(La.X + u * LbMinusLa.X, La.Y + u * LbMinusLa.Y, La.Z + u * LbMinusLa.Z);
            return IntersectingPoint;
        }
        public static Point3D CrossProduct(Point3D p0, Point3D p1, Point3D p2) //3 points constituting the plane
        {//parallel planes have the same normal
            Point3D A = new Point3D(p1.X - p0.X, p1.Y - p0.Y, p1.Z - p0.Z);
            Point3D B = new Point3D(p2.X - p0.X, p2.Y - p0.Y, p2.Z - p0.Z);
            double X = A.Y * B.Z - A.Z * B.Y;
            double Y = A.Z * B.X - A.X * B.Z;
            double Z = A.X * B.Y - A.Y * B.X;
            return new Point3D(X, Y, Z);
        }
        public static Point3D CrossProduct(Point3D A, Point3D B) //A and B are vectors ie magnitude and direction same as Point3D
        {
            double X = A.Y * B.Z - A.Z * B.Y;
            double Y = A.Z * B.X - A.X * B.Z;
            double Z = A.X * B.Y - A.Y * B.X;
            return new Point3D(X, Y, Z);
        }
        public static double DotProduct(Point3D A, Point3D B)
        {
            double answer = A.X * B.X + A.Y * B.Y + A.Z * B.Z;
            return answer;
        }