reg3 <- function(d.lm, model, pcol = 2, pch=16, xlim=NULL, ylim=NULL, N = 20, ...)
{
z = d.lm$model[,1]
x = d.lm$model[,2]
y = d.lm$model[,3]
k = coef(d.lm)
if(is.null(xlim))
xlim = range(x)
if(is.null(ylim))
ylim = range(y)
X = seq(xlim[1], xlim[2], length=N)
Y = seq(ylim[1], ylim[2], length=N)
Z = outer(X, Y, model, k)
persp(X, Y, Z, ...) -> res
points(trans3d(x, y, z, pmat=res),  col = pcol, pch =pch)
}