## Do this after oak1.s and oak102.s.
##
######## 1) Pocket plot
## e = unit vector pointing north
zz = x1 ## we will do this on the raw totals, x1, 
## but could do this on the smoothed totals, a1, instead.
n1 = nrow(zz) ## we will show p[i,1] up to p[i,5]
y = matrix(0,nrow=n1,ncol=n1)
for(i in c(1:(n1-1))){
for(j in c((i+1):n1)){
y[i,j] = mean(sqrt(abs(zz[i,] - zz[j,])))
y[j,i] = y[i,j]
}
}

z = rep(0,(n1-1))  ### z[k] will represent y-double-bar[k]
for(k in c(1:(n1-2))){
cat(k,"\n")
temp = y[1,(1+k)]
for(i in c(2:(n1-k))){
temp = c(temp,y[i,(i+k)])
}
z[k] = mean(temp)
}
z[n1-1] = y[1,n1]

p = matrix(rep(0,(n1*n1)),ncol=n1)
for(i in c(1:(n1-1))){
for(j in c((i+1):n1)){
k = j-i
p[i,j] = y[i,j] - z[k]
p[j,i] = p[i,j]
}
}

par(mfrow=c(1,1))
plot(c(1,n1),range(p),xlab="Row number",ylab="P_jk",type="n",xaxt="n")
boxplot(p[,1],p[,2],p[,3],p[,4],p[,5],p[,6],p[,7],p[,8],p[,9],p[,10],p[,11],
    p[,12],p[,13],p[,14],
## I don't know how to do this other than by writing all of these out.
## If n1 > 10, then you have to keep going....
## For instance if x1 had 30 rows, you'd say
## boxplot(p[,1],p[,2],p[,3],p[,4],p[,5],p[,6],p[,7],p[,8],p[,9],p[,10],
## p[,11],p[,12],p[,13],p[,14],p[,15],p[,16],p[,17],p[,18],p[,19],p[,20],
## p[,21],p[,22],p[,23],p[,24],p[,25],p[,26],p[,27],p[,28],p[,29],p[,30],
names=as.character(seq(1,n1)),add=T)
mtext(s=3,"Pocket Plot. e = unit vector pointing North")

####### 2) VARIOGRAMS AND COVARIOGRAMS:
sill1 = var(c(zz))
data2 = data.frame(x=rep(-119+c(1:14)*.1-.05,times=rep(10,14)),
    y=rep(33.7+c(1:10)*.1-.05,14),z=c(t(zz)))
library(spatial) ## you might have to do install.packages("spatial") first
my.ls0 = surf.ls(0,data2)
c1 = correlogram(my.ls0,100,plot=F)  ## 100 lags -- you can change this
plot(c1$x, c1$y * sill1, xlab="distance", ylab="covariance",type="l") ## covariogram
m1 = max((1:length(c1$x))[c1$x<.5]) ## to only go up to distance 0.5
plot(c1$x[1:m1], c1$y[1:m1] * sill1, xlab="distance", ylab="covariance",type="l") 
plot(c1$x[1:m1], c1$y[1:m1], xlab="distance", ylab="correlation",type="l") ## correlogram
c2 = variogram(my.ls0,100,plot=F)  ## again, 100 lags -- you may change this.
plot(c2$x[1:m1], c2$y[1:m1], xlab="distance, h", ylab="gamma(h)",type="l") ## semi-variogram
my.ls1 = surf.ls(np=1,data2)
c3 = variogram(my.ls1,100,plot=F)  ## fits semi-variogram to LS residuals
plot(c3$x[1:m1], c3$y[1:m1], xlab="distance, h", ylab="gamma(h)",type="l") ## semi-variogram
my.matrix.ls = trmat(my.ls1,min(data2$x),max(data2$x),min(data2$y),max(data2$y),40)
image(zz,col=grey((64:20)/64))
mtext(s=3,"data")
image(my.matrix.ls,col=grey((64:20)/64))
mtext(s=3,"linear fit to data")

## Non-parametric and variogram estimates, using correlogram:
par(mfrow=c(2,2))
c1 = correlogram(my.ls0,nint=100),plot=F) 
m1 = max((1:length(c1$x))[c1$x<.5]) ## to only go up to distance 0.5 in what follows
a1 = seq(0.01,0.5,length=50)
plot(c1$x[1:m1], 2*sill1 - 2*c1$y[1:m1] * sill1, 
    xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=2,main="exponential, d=5") 
lines(a1, 2*sill1 - 2*expcov(a1,5)*sill1)
plot(c1$x[1:m1], 2*sill1 - 2*c1$y[1:m1] * sill1, xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=2,main="exponential, d=0.5") 
lines(a1, 2*sill1 - 2*expcov(a1,0.5)*sill1)
plot(c1$x[1:m1], 2*sill1 - 2*c1$y[1:m1] * sill1, xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=2,main="exponential, d=0.1") 
lines(a1, 2*sill1 - 2*expcov(a1,0.1)*sill1)
plot(c1$x[1:m1], 2*sill1 - 2*c1$y[1:m1] * sill1, xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=2,main="exponential, d=0.01") 
lines(a1, 2*sill1 - 2*expcov(a1,0.01)*sill1)

par(mfrow=c(2,2))
plot(c1$x[1:m1], 2*sill1 - 2*c1$y[1:m1] * sill1, xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=2,main="Gaussian, d=5") 
lines(a1, 2*sill1 - 2*gaucov(a1,5)*sill1)
plot(c1$x[1:m1], 2*sill1 - 2*c1$y[1:m1] * sill1, xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=2,main="Gaussian, d=0.5") 
lines(a1, 2*sill1 - 2*gaucov(a1,0.5)*sill1)
plot(c1$x[1:m1], 2*sill1 - 2*c1$y[1:m1] * sill1, xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=2,main="Gaussian, d=0.1") 
lines(a1, 2*sill1 - 2*gaucov(a1,0.1)*sill1)
plot(c1$x[1:m1], 2*sill1 - 2*c1$y[1:m1] * sill1, xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=2,main="Gaussian, d=0.01") 
lines(a1, 2*sill1 - 2*gaucov(a1,0.01)*sill1)

par(mfrow=c(2,2))
plot(c1$x[1:m1], 2*sill1 - 2*c1$y[1:m1] * sill1, xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=2,main="spherical, d=5") 
lines(a1, 2*sill1 - 2*sphercov(a1,5)*sill1)
plot(c1$x[1:m1], 2*sill1 - 2*c1$y[1:m1] * sill1, xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=2,main="spherical, d=0.5") 
lines(a1, 2*sill1 - 2*sphercov(a1,0.5)*sill1)
plot(c1$x[1:m1], 2*sill1 - 2*c1$y[1:m1] * sill1, xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=2,main="spherical, d=0.1") 
lines(a1, 2*sill1 - 2*sphercov(a1,0.1)*sill1)
plot(c1$x[1:m1], 2*sill1 - 2*c1$y[1:m1] * sill1, xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=2,main="spherical, d=0.01") 
lines(a1, 2*sill1 - 2*sphercov(a1,0.01)*sill1)

## Compare all 3 with d = 0.1:
par(mfrow=c(1,1))
plot(c1$x[1:m1], 2*sill1 - 2*c1$y[1:m1] * sill1, xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=2,main="variograms, d=0.1") 
lines(a1, 2*sill1 - 2*expcov(a1,0.1)*sill1,col=3)
lines(a1, 2*sill1 - 2*gaucov(a1,0.1)*sill1,col=4)
lines(a1, 2*sill1 - 2*sphercov(a1,0.1)*sill1,col=5)
legend(0.3,.9,c("nonparametric","exponential", "Gaussian", "spherical"),
    col=c(1,3,4,5),lty=c(2,1,1,1))

mp1 = med.polish(x1,1)
res = mp1[[1]]
par(mfrow=c(1,2))
plot(c1$x[1:m1], 2*sill1 - 2*c1$y[1:m1] * sill1, xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=2,main="variogram of data") 
data3 = data.frame(x=rep(1:14,times=rep(10,14)),y=rep(1:10,14),z=c(t(mp1[[1]])))
my.ls0 = surf.ls(0,data3)
c5 = variogram(my.ls0,100,plot=F) 
plot(c1$x[1:m1], 2*c5$y[1:m1], type="l",xlab="distance, h", ylab="2gamma(h)",lty=2) 
title(main = "variogram of \n median polish residuals")