solving set of linear equations using R for plane 3D equation

I have some trouble in order to solve my set of linear equations.

I have three 3D points (A, B, C) in my example and I want to automate the solving of my system. I want to create a plane with these 3 points. It's very simple manually (mathematically) but I don't see why I don't solve my problem when I code...

I have a system of cartesian equation which is the equation of a plane : ax+by+cz+d=0

xAx + yAy + zA*z +d = 0 #point A

xBx + yBy + zB*z +d = 0 #point B

etc

I use a matrix, for example A=(0,0,1) ; B=(4,2,3) and C=(-3,1,0). With manual solving, I have for this example this solution : x+3y-5z+5=0.

For resolving it in R : I wanted to use solve().

A <- c(0,0,1)
B <- c(4,2,3)
C <- c(-3,1,0)
res0 <- c(-d,-d,-d) #I don't know how having it so I tried c(0,0,0) cause each equation = 0. But I really don't know for that !

#' @param A vector 3x1 with the 3d coordinates of the point A
carteq <- function(A, B, C, res0) {
   matrixtest0 <- matrix(c(A[1], A[2], A[3], B[1], B[2], B[3],C[1], C[2], C[3]), ncol=3) #I tried to add the 4th column for solving "d" but that doesn't work.
   #checking the invertibility of my matrix
   out <- tryCatch(determinant(matrixtest0)$modulus<threshold, error = function(e) e)#or out <- tryCatch(solve(X) %*% X, error = function(e) e)
   abcd <- solve(matrixtest0, res0) #returns just 3 values
   abcd <- qr.solve(matrixtest0, res0) #returns just 3 values
}

That's not the good method... But I don't know how I can add the "d" in my problem.

The return that I need is : return(a, b, c, d)

I thing that my problem is classical and easy, but I don't find a function like solve() or qr.solve() which can solve my problem...

Solution

Your solution is actually wrong:

A <- c(0,0,1)
B <- c(4,2,3)
C <- c(-3,1,0)

CrossProduct3D <- function(x, y, i=1:3) {
  #http://stackoverflow.com/a/21736807/1412059

  To3D <- function(x) head(c(x, rep(0, 3)), 3)
  x <- To3D(x)
  y <- To3D(y)

  Index3D <- function(i) (i - 1) %% 3 + 1

  return (x[Index3D(i + 1)] * y[Index3D(i + 2)] -
            x[Index3D(i + 2)] * y[Index3D(i + 1)])
}



N <- CrossProduct3D(A - B, C - B)
#[1]   4   2 -10
d <- -sum(N * B)
#[1] 10

#test it:

crossprod(A, N) + d
#     [,1]
#[1,]    0

crossprod(B, N) + d
#     [,1]
#[1,]    0

crossprod(C, N) + d
#     [,1]
#[1,]    0