[R] how to create model matrix of second order terms

[Stephen Bond]

I am sending to this forum as stackoverflow has devolved into sth
pretty bad.
Below code shows how to get what I want in a clumsy way.

cols <- letters[1:4]
a1 <- outer(cols,cols,paste0)
b1 <- a1[!lower.tri(a1)]

X <- matrix(rnorm(80),ncol=4)
colnames(X) <- cols
X <- as.data.frame(X)
XX <- matrix(0,nrow=nrow(X),ncol=length(b1))
colnames(XX) <- b1

for (k in 1:length(b1)){
    XX[,k] <- X[,substr(b1[k],1,1)]*X[,substr(b1[k],2,2)]
}



Is there a way to get that using a formula or some neat trick? The
above will not work for factors, so I will need to create the factor
crossings using formula a*b*c and then cross with the numerics, which
is even more clumsy.
Thanks everybody