setwd("D:/PhD-cursus/RN")


results <- data.frame(corr = as.numeric(), st = as.numeric(), run = as.character())
results.h2 <- data.frame(corr = as.numeric(), st = as.numeric(), run = as.character())

run <- c("RN") ## name of your .as file
for (x in run){
varx1 <- 0 ## environmental factor 0
varx11 <- varx1^2 ## environmental factor squared for later calculation of variance and covariance
varx12 <- varx1*2 ## environmental factor multiplied by 2 for later calculation of variance and covariance
varx2 <- seq(-3,3,1.0) ## the range of environemtal values on the x-axis
varx21 <- varx2^2
varx22 <- varx2*2
varx3 <- varx1+varx2
varx4 <- varx1*varx2

pin1  <- paste("F var1 2 + 3 *", varx12, " + 4 *",varx11, sep = "") ## additive genetic variance given the environmental factor zero
pin2  <- paste("F covar 2 + 3 *", varx3, " + 4 *",varx4, sep = "") ## additive genetic covariance between the environment zero and any other environment
pin3  <- paste("F var2 2 + 3 *", varx22, " + 4 *",varx21, sep = "") ## additive genetic variance for each point of environmental factor

pin4  <- paste("R cor 5 ",6:12," ",13:19," ", sep = "") ## genetic correlation between the trait at environment zero and other environments
    pin <- c(pin1,pin2,pin3,pin4)


 write.table(pin, file = paste(x,".pin", sep = ""),quote=F,col.names=F,row.names=F)
  shell(paste(x,".pin", sep = ""))
  readLines(paste(x,".pvc", sep = ""))
}

# run first the pin-file in asreml and the part below to make the graph

resa <- (as.data.frame(as.numeric(substr(readLines(paste(x,".pvc", sep = ""))[25:31],56,63 )))  )
resb <- (as.data.frame(as.numeric(substr(readLines(paste(x,".pvc", sep = ""))[25:31],67,73 )))  )


res <- data.frame(resa,resb)
colnames(res) <- c("corr", "st")
res$run = x
res$low <- resa-resb
res$hi <- resa + resb

res$x <- varx2

require(lattice)
tiff("correlation.tif", width = 17.8, height = 11.4 , units = "cm",res = 300,family = "serif")
xyplot(corr ~ x  , data = res[res$run == "RN",],ylim = c(0,1.1), scales = list(cex = 1.5),
       col = "black",pch = 16,
       panel = function(x,y,...) {
         lims <- current.panel.limits()
         panel.abline(h= 0,v=0, lwd = 2 )
         panel.segments(x0 = res$x, y0 = res$cor - res$st, x1 = res$x, y1 = res$cor + res$st, col = "darkgrey")
         panel.xyplot(x,y,...)}
       ,par.settings = list( axis.line = list(col = "white")),
       xlab = list(label = "x",cex=1.5), ylab = list(label = "Correlation",cex=1.5)
)
dev.off()