#
###################################################
# SCRIPT FOR TEMPLATE OF SCATTER PLOT   
# Created by Gottfried Pestal (Solv Consulting Ltd.)
# Version 1 - August 23, 2012
###################################################

# This script illustrates how to create a pdf file with 10 pages and 3 figures on each page.


# create an array and fill it with parameters for plotting (2 random normal variables, mean and sd changing over time)
sample.pars <- array(NA,dim=c(10,2,2), dimnames=list(1980:1989,c("mean","sd"),c("Var_1","Var_2")))
sample.pars[,"mean","Var_1"] <- seq(-2,2,length.out=10)
sample.pars[,"sd","Var_1"] <- seq(1,3,length.out=10)
sample.pars[1:5,"mean","Var_2"] <- seq(-1,2,length.out=5)
sample.pars[6:10,"mean","Var_2"] <- seq(2,-2,length.out=5)
sample.pars[,"sd","Var_2"] <- 1


# open the pdf file as a plotting device
# note: onefile = TRUE creates a new page in the same pdf file whenever the current device is full
pdf("FigRs_Templates_ScatterPlot.pdf",width=11, height=8.5,onefile=TRUE)


# specify margins 
par(omi=c(0.5,0.5,0.5,0.5) , mai=c(0.9,0.7,0.2,0.2))

# divide the plotting device 3 into panels
layout(matrix(c(1,2,1,3), byrow=TRUE, nrow=2,ncol=2), heights=c(1,1),widths=c(1.7,1))



# loop through the 10 years
for(year in 1:10){

# sample random data for this year
var1.tmp <- rnorm(500,mean=sample.pars[year,"mean","Var_1"],sd=sample.pars[year,"sd","Var_1"])
var2.tmp <- rnorm(500,mean=sample.pars[year,"mean","Var_2"],sd=sample.pars[year,"sd","Var_2"])


# create plot in panel 1
plot(var1.tmp,var2.tmp,xlim=c(-6,6), ylim=c(-6,6),pch=19,col="lightblue", xlab="Variable 1",ylab="Variable 2",bty="n",cex=1.2,cex.lab=1.6,cex.axis=1.6)

# add lines with value of 3 as a benchmark
abline(h=3,col="red")
abline(v=3,col="red")

#redraw points above benchmark in red
bm.flag<-var1.tmp>3 | var2.tmp>3
points(var1.tmp[bm.flag],var2.tmp[bm.flag],pch=19,col="red", cex=1.2)


# create plot in panel 2
hist(var1.tmp,xlim=c(-6,6),freq=FALSE, ylim=c(0,0.5),breaks=c(-100,seq(-6,6,by=0.5),100),col="red",border="red", axes=FALSE,xlab="",ylab="",main="Variable 1")
hist(var1.tmp[var1.tmp<3],xlim=c(-6,3),freq=FALSE, ylim=c(0,0.5),breaks=c(-100,seq(-6,6,by=0.5),100),add=TRUE,col="lightblue",border="lightblue")
axis(1,at=seq(-6,6,by=2),cex.axis=1.3)

# create plot in panel 3
hist(var2.tmp,xlim=c(-6,6),freq=FALSE, ylim=c(0,0.5),breaks=c(-100,seq(-6,6,by=0.5),100),col="red",border="red", axes=FALSE,xlab="",ylab="",main="Variable 2")
hist(var2.tmp[var2.tmp<3],xlim=c(-6,3),freq=FALSE, ylim=c(0,0.5),breaks=c(-100,seq(-6,6,by=0.5),100),add=TRUE,col="lightblue",border="lightblue")
axis(1,at=seq(-6,6,by=2),cex.axis=1.3)

# add vertical line to second plot, expanding them across the first plot
abline(v=3,xpd=NA, col="red")



# label the page
mtext(year+1979,side=3,line=0.5, outer=TRUE, cex=2.5)
mtext("www.solv.ca/FigRs",side=1,line=0.5, outer=TRUE, cex=1)



} # end looping across year


# turn off plotting device and save pdf file
dev.off()


#