rm(list=ls())
par(mfrow=c(1,1))
n.points <- 4
#
# das Kovarianz- und Variogrammmodell (exponentielles Modell)
#
cov.mod <- function(x){        # Kovarianzmodell
  cov.mod <- exp(-abs(x)/5)
}
gamma.mod <- function(x) {     # Variogrammmodell
  gamma.mod <- cov.mod(0) - cov.mod(x)
}
#
plot(1,1,type="n",xlim=c(0,10),ylim=c(0,10),xlab="x",ylab="y")
# Koordinaten des zu schätzenden Punktes:
y.coords <- c(4,5)
points(t(y.coords),cex=2.5,pch=20,col="red")
#
# Anklicken der Punkte des "Messgitters"
#
p.x <- 1:n.points
p.y <- 1:n.points
for (i in 1:n.points){
  p <- locator(1)
  p.x[i] <- p$x
  p.y[i] <- p$y
  points(p,cex=2.5,pch=20,col="darkgreen")
}
p.coords <- matrix(c(p.x,p.y),nrow=2,byrow=T)
#
# Aufbau der Krigematrix und der rechten Seite des LGS
#
D <- matrix(1:(n.points+1)^2,ncol=n.points+1)
rs <- 1:(n.points+1)
for(i in 1:n.points) {
  for(j in 1:n.points) {
    dist <- p.coords[,i]-p.coords[,j]
    D[i,j] <- gamma.mod(sqrt(dist%*%dist))
  }
  dist <- y.coords-p.coords[,i]
  rs[i] <- gamma.mod(sqrt(dist%*%dist))
}
D[n.points+1,] <- 1; D[,n.points+1] <- 1; D[n.points+1,n.points+1] <- 0
rs[n.points+1] <- 1

#
# Berechnung der Wichtungsfaktoren
#
D.inv <- solve(D)
lambda <- D.inv%*%rs
cat("\n\n  Ordinary Kriging: Lösung des LGS:   \n")
cat("--------------------------------------\n")
cat("\n Wichtungsfaktoren: \n")
print(as.matrix(lambda[1:n.points]))
cat("\n Lagrange multiplier: \n")
print(lambda[n.points+1])
cat("\n Summe der Wichtungsfaktoren: \n")
print(round(sum(lambda[1:n.points]),3)) 
#
# Berechnung des Schätzwertes und der Schätzvarianz:
#
T.long <- c(18,20,17,13,18,23,19,22,16,25) 
# wenn n.points vergrößert wird, müssen hier mehr Werte stehen
T <- T.long[1:n.points]
T.val <- as.real(T%*%lambda[1:n.points])
cat("\n Schätzwert: \n")
print(T.val)
#
T.var <- as.real(rs%*%lambda)
cat("\n Schätzvarianz: \n")
print(T.var)
#
cat("\n++++ Ende Ordinary Kriging  ++++\n")
