#dev.new(width=9,height=5) # only needed once
pdf("f08_res_vs_co3_dep.pdf", width = 9, height = 5 ,encoding="WinAnsi")
par(mfrow=c(1,2))
par(mai=c(0.8,0.9,0.2,0.0))

infile <- "data_dx1000_dz5_1.1.dat"
source("read_mapped_data.R")
deltanatAll <- c13sednatb - c13wcnatb
lm.deltanat.co3All <- lm(deltanatAll ~ co3b)
deltanatLH <- c13sednatLHb - c13wcnatLHb
lm.deltanat.co3LH <- lm(deltanatLH ~ co3LHb)

xlm <- c(0,280)
ylm <- c(-1.2,1.2)

plot(co35,deltanat5,xlim=xlm,ylim=ylm,ylab=expression(paste(delta^13*C[Cibnat]-delta^13*C[DICnat]," (\u2030)")),
     xlab=expression(paste(CO[3]^'2-'," ",(mu*M))), cex = 0.3, pch = 21, col="orange", bg = "orange")
abline(0.76,-0.006,col="lightblue")
abline(1.55,-0.014,col="lightblue")
abline(0.15,-0.001,col="lightgreen")
par(new=T)
abline(lm.deltanat.co3All,lwd=1)
par(new=T)
abline(lm.deltanat.co3LH,lwd=1)
par(new=T)
plot(co35,deltanat5,xlim=xlm,ylim=ylm,xlab="",ylab="", cex = 0.3, pch = 21, col="orange", bg = "orange")
par(new=T)
plot(co38,deltanat8,xlim=xlm,ylim=ylm,col="blue",xlab="",ylab="",pch=21,bg="blue",cex=0.3)
par(new=T)
plot(co36,deltanat6,xlim=xlm,ylim=ylm,col="green",xlab="",ylab="",pch=21,bg="green",cex=0.3)
par(new=T)
plot(co31,deltanat1,xlim=xlm,ylim=ylm,col="black",xlab="",ylab="",pch=21,bg="black",cex=0.3)
par(new=T)
plot(co32,deltanat2,xlim=xlm,ylim=ylm,col="purple",xlab="",ylab="",pch=21,bg="purple",cex=0.3)
par(new=T)
plot(co37,deltanat7,xlim=xlm,ylim=ylm,col="cyan",xlab="",ylab="",pch=21,bg="cyan",cex=0.3)
par(new=T)
plot(co33,deltanat3,xlim=xlm,ylim=ylm,col="red",xlab="",ylab="",pch=21,bg="red",cex=0.3)
par(new=T)
plot(co39,deltanat9,xlim=xlm,ylim=ylm,col="turquoise",xlab="",ylab="",pch=21,bg="turquoise",cex=0.3)
par(new=T)
plot(co310,deltanat10,xlim=xlm,ylim=ylm,col="orchid1",xlab="",ylab="",pch=21,bg="orchid1",cex=0.3)
par(new=T)
plot(co311,deltanat11,xlim=xlm,ylim=ylm,col="brown",xlab="",ylab="",pch=21,bg="brown",cex=0.3)
par(new=T)
plot(co312,deltanat12,xlim=xlm,ylim=ylm,col="green3",xlab="",ylab="",pch=21,bg="green3",cex=0.3)
par(new=T)
plot(co313,deltanat13,xlim=xlm,ylim=ylm,col="seagreen",xlab="",ylab="",pch=21,bg="seagreen",cex=0.3)
par(new=T)
plot(co314,deltanat14,xlim=xlm,ylim=ylm,col="red3",xlab="",ylab="",pch=21,bg="red3",cex=0.3)
par(new=T)
plot(co315,deltanat15,xlim=xlm,ylim=ylm,col="yellow3",xlab="",ylab="",pch=21,bg="yellow3",cex=0.3)
par(new=T)
plot(co316,deltanat16,xlim=xlm,ylim=ylm,col="orange3",xlab="",ylab="",pch=21,bg="orange3",cex=0.3)
par(new=T)
plot(co317,deltanat17,xlim=xlm,ylim=ylm,col="orange2",xlab="",ylab="",pch=21,bg="orange2",cex=0.3)
par(new=T)
plot(co318,deltanat18,xlim=xlm,ylim=ylm,col="yellow2",xlab="",ylab="",pch=21,bg="yellow2",cex=0.3)

text(290,1.2,"Spero et al. (1997)",pos=2,col="lightblue")

par(mai=c(0.8,0.6,0.2,0.3))
xlm <- c(0,5.6)
plot(dep5/(1000),deltanat5,xlim=xlm,ylim=ylm,ylab="",
     xlab=expression("depth (km)"), cex = 0.3, pch = 21, col="orange", bg="orange")

abline(0.03,-0.02,col="lightgreen")
abline(0.05,-0.03,col="lightgreen")
par(new=T)
plot(dep5/(1000),deltanat5,xlim=xlm,ylim=ylm,xlab="",ylab="", cex = 0.3, pch = 20, col="orange", bg="orange")

depbkmAll = depb/1000
lm.deltanat.depAll <- lm(deltanatAll ~ depbkmAll)
abline(lm.deltanat.depAll,lwd=1)
depbkmLH = depLHb/1000
lm.deltanat.depLH <- lm(deltanatLH ~ depbkmLH)
abline(lm.deltanat.depLH,lwd=1)
text(5.8,1.2,"Hesse et al. (2014)",pos=2,col="lightgreen")

par(new=T)
plot(dep8/1000,deltanat8,xlim=xlm,ylim=ylm,col="blue",xlab="",ylab="",pch=21,bg="blue",cex=0.3)
par(new=T)
plot(dep6/1000,deltanat6,xlim=xlm,ylim=ylm,col="green",xlab="",ylab="",pch=21,bg="green",cex=0.3)
par(new=T)
plot(dep1/1000,deltanat1,xlim=xlm,ylim=ylm,col="black",xlab="",ylab="",pch=21,bg="black",cex=0.3)
par(new=T)
plot(dep2/1000,deltanat2,xlim=xlm,ylim=ylm,col="purple",xlab="",ylab="",pch=21,bg="purple",cex=0.3)
par(new=T)
plot(dep7/1000,deltanat7,xlim=xlm,ylim=ylm,col="cyan",xlab="",ylab="",pch=21,bg="cyan",cex=0.3)
par(new=T)
plot(dep3/1000,deltanat3,xlim=xlm,ylim=ylm,col="red",xlab="",ylab="",pch=21,bg="red",cex=0.3)
par(new=T)
plot(dep9/1000,deltanat9,xlim=xlm,ylim=ylm,col="turquoise",xlab="",ylab="",pch=21,bg="turquoise",cex=0.3)
par(new=T)
plot(dep10/1000,deltanat10,xlim=xlm,ylim=ylm,col="orchid1",xlab="",ylab="",pch=21,bg="orchid1",cex=0.3)
par(new=T)
plot(dep11/1000,deltanat11,xlim=xlm,ylim=ylm,col="brown",xlab="",ylab="",pch=21,bg="brown",cex=0.3)
par(new=T)
plot(dep12/1000,deltanat12,xlim=xlm,ylim=ylm,col="green3",xlab="",ylab="",pch=21,bg="green3",cex=0.3)
par(new=T)
plot(dep13/1000,deltanat13,xlim=xlm,ylim=ylm,col="seagreen",xlab="",ylab="",pch=21,bg="seagreen",cex=0.3)
par(new=T)
plot(dep14/1000,deltanat14,xlim=xlm,ylim=ylm,col="red3",xlab="",ylab="",pch=21,bg="red3",cex=0.3)
par(new=T)
plot(dep15/1000,deltanat15,xlim=xlm,ylim=ylm,col="yellow3",xlab="",ylab="",pch=21,bg="yellow3",cex=0.3)
par(new=T)
plot(dep16/1000,deltanat16,xlim=xlm,ylim=ylm,col="orange3",xlab="",ylab="",pch=21,bg="orange3",cex=0.3)
par(new=T)
plot(dep17/1000,deltanat17,xlim=xlm,ylim=ylm,col="orange2",xlab="",ylab="",pch=21,bg="orange2",cex=0.3)
par(new=T)
plot(dep18/1000,deltanat18,xlim=xlm,ylim=ylm,col="yellow2",xlab="",ylab="",pch=21,bg="yellow2",cex=0.3)
dev.off()