rm(list=ls())
#
# das Kovarianzmodell (exponentielles Modell)
#
cov.mod <- function(h){
  cov.mod <- exp(-abs(h)/5)
}
#
# Koordinaten des zu schätzenden Punktes:
#
y.coords <- c(3,4)
#
# Punkte des "Messgitters"
#
p.x <- c(2,4,0)
p.y <- c(6,4,2)
p.coords <- matrix(c(p.x,p.y),nrow=2,byrow=T)
#
# Aufbau der Krigematrix und der rechten Seite des LGS
#
D <- matrix(1:9,ncol=3)
dist <- p.coords[,1]-p.coords[,1]; D[1,1] <- cov.mod(sqrt(dist%*%dist))
dist <- p.coords[,1]-p.coords[,2]; D[1,2] <- cov.mod(sqrt(dist%*%dist))
dist <- p.coords[,1]-p.coords[,3]; D[1,3] <- cov.mod(sqrt(dist%*%dist))
dist <- p.coords[,2]-p.coords[,1]; D[2,1] <- cov.mod(sqrt(dist%*%dist))
dist <- p.coords[,2]-p.coords[,2]; D[2,2] <- cov.mod(sqrt(dist%*%dist))
dist <- p.coords[,2]-p.coords[,3]; D[2,3] <- cov.mod(sqrt(dist%*%dist))
dist <- p.coords[,3]-p.coords[,1]; D[3,1] <- cov.mod(sqrt(dist%*%dist))
dist <- p.coords[,3]-p.coords[,2]; D[3,2] <- cov.mod(sqrt(dist%*%dist))
dist <- p.coords[,3]-p.coords[,3]; D[3,3] <- cov.mod(sqrt(dist%*%dist))
#
# Die reche Seite des Gleichungssystems
#
rs <- 1:3
dist <- p.coords[,1]-y.coords; rs[1] <- cov.mod(sqrt(dist%*%dist))
dist <- p.coords[,2]-y.coords; rs[2] <- cov.mod(sqrt(dist%*%dist))
dist <- p.coords[,3]-y.coords; rs[3] <- cov.mod(sqrt(dist%*%dist))
#
# Berechnung der Wichtungsfaktoren
#
D.inv <- solve(D)
lambda <- D.inv%*%rs
cat("\n\n  Simple Kriging: Lösung des LGS:   \n")
cat("------------------------------------\n")
cat("\n Wichtungsfaktoren: \n")
print(lambda)
#
cat("\n Summe der Wichtungsfaktoren: \n")
print(sum(lambda)) 

# Einbeziehen von Messwerten:
T.mean = 19
T <- c(18,20,17) # wenn n.points vergrößert wird, müssen hier mehr Werte stehen
#
# Berechnung des Schätzwertes und der Schätzvarianz:
#
T.val <- as.real((T-T.mean)%*%lambda + T.mean)
cat("\n Schätzwert: \n")
print(T.val)
#
T.var <- cov.mod(0) - as.real(rs%*%lambda)
cat("\n Schätzvarianz: \n")
print(T.var)
cat("\n++++ finito ++++\n")

