public static int PointOfIntersection(Point3D La, Point3D Lb, Point3D P0, Point3D P1, Point3D P2, out Point3D intersection)
{
    //thanks to Paul Bourke and Bryan Hanson
    //la,Lb constitute the line. P0,P1,P2 constitute the plane
    intersection = new Point3D(0, 0, 0);
    Vector3D N = CrossProduct(P0, P1, P2);  //N is the normal of the plane        
    double n = DotProduct(N, new Vector3D(P2.X - La.X, P2.Y - La.Y, P2.Z - La.Z));
    Vector3D LbMinusLa = new Vector3D(Lb.X - La.X, Lb.Y - La.Y, Lb.Z - La.Z);
    double d = DotProduct(N, LbMinusLa);
    if (d == 0)
        return -1;//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
    intersection = new Point3D(La.X + u * LbMinusLa.X, La.Y + u * LbMinusLa.Y, La.Z + u * LbMinusLa.Z);
    return 1;
}
public static Vector3D CrossProduct(Point3D p0, Point3D p1, Point3D p2) //3 points constituting the plane
{//parallel planes have the same normal
    Vector3D A = new Vector3D(p1.X - p0.X, p1.Y - p0.Y, p1.Z - p0.Z);
    Vector3D B = new Vector3D(p2.X - p0.X, p2.Y - p0.Y, p2.Z - p0.Z);
    return CrossProduct(A, B);
}
public static Vector3D CrossProduct(Vector3D A, Vector3D B) //A and B are vectors ie magnitude and direction only
{
    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 Vector3D(X, Y, Z);
}
public static double DotProduct(Vector3D A, Vector3D B)
{
    double answer = A.X * B.X + A.Y * B.Y + A.Z * B.Z;
    return answer;
}