function z = rayplane(r0,rd,p0,p1,origin)
    v1 = p0-origin;
    v2 = p1-origin;
    N = cross(p0,p1);
    N = N/sqrt(dot(N,N));
    D = -dot(N,origin);
    t = -(dot(N,r0)+D)/(dot(N,rd));
    z = t;