cov.th.e2 <- function(d,par,var.z) { # covariance exponentielle, 2 paramètres + pépite
	eps = 10 
	cov = par[1]*exp(-(d/par[2]))
	wd = which(d<eps)
	cov[wd] = var.z # pepite
	return(cov)
}


ll_gauss_mult2 = function(par,x,y,z,covariance=cov.th.e) {
	n = length(z)
	C = matrix(0,n,n)
	for (i in 1:n) {
		d = m_lldist_p(x,y,x[i],y[i])
		C[i,] = covariance(d,par,var.z=var(z)) # deux paramètres
	}
	#browser()
	z = matrix(z,n,1)
	ll =  -log(2*pi)*n/2-log(det(C))/2-t(z-mu)%*%solve(C)%*%(z-mu)/2
	return(-ll)
}

par0 = c(16,300)
par.opt = optim(par0,ll_gauss_mult2,gr=NULL,x=x,y=y,z=z,method="L-BFGS-B",lower=c(1,50),upper=c(var(z),500))