#dev.new(width=5.5,height=9) # only needed once
pdf("f10_res_vs_ftd_species.pdf", width = 5.5, height = 9 ,encoding="WinAnsi")
par(mfrow=c(2,1))

infile <- "data_dx1000_dz5_1.1.dat"
#infile <- "data_dx500_dz10_1.1.dat"
source("read_mapped_data.R")
lm.deltanat.co3.depAll <- lm(c13sednatb - c13wcnatb ~ co3b + depb)
ftd <- fitted(lm.deltanat.co3.depAll)
ftd <- c13wcnatb
res <- resid(lm.deltanat.co3.depAll)

par(mai=c(0.4,0.9,0.4,0.2))
print("LA1")
datasetplot <- datasetb
xlm = c(-0.4,1.8)
ylm = c(-1.4,1.4)
plot(ftd,res,xlim=xlm,ylim=ylm,col="black",pch=21,bg="black",cex=0.7,
     ylab="Residuals (\u2030)",
     xlab=expression(paste(delta^13*C[DICnat], " (\u2030)"))
    )
par(new=T)
reswue <- res[ispecwueb]
plot(ftd[ispecwueb],res[ispecwueb],xlim=xlm,ylim=ylm,xlab="",ylab="",col="green",pch=21,bg="green",cex=0.7)
print(paste("Cib All: ",mean(res^2, na.rm=T)^.5,mean(res,na.rm=T),length(res[!is.na(res)])))
print(paste("Cib Wue: ",mean(reswue^2, na.rm=T)^.5,mean(reswue,na.rm=T),length(reswue[!is.na(reswue)])))

par(new=T)
reskul <- res[ispeckulb]
plot(ftd[ispeckulb],res[ispeckulb],xlim=xlm,ylim=ylm,xlab="",ylab="",col="red",pch=21,bg="red",cex=0.7)
print(paste("Cib Kul: ",mean(reskul^2, na.rm=T)^.5,mean(reskul,na.rm=T),length(reskul[!is.na(reskul)])))

par(new=T)
reslob <- res[ispeclobb]
plot(ftd[ispeclobb],res[ispeclobb],xlim=xlm,ylim=ylm,xlab="",ylab="",col="blue",pch=21,bg="blue",cex=0.7)
print(paste("Cib Lob: ",mean(reslob^2, na.rm=T)^.5,mean(reslob,na.rm=T),length(reslob[!is.na(reslob)])))

par(new=T)
respac <- res[ispecpacb]
plot(ftd[ispecpacb],res[ispecpacb],xlim=xlm,ylim=ylm,xlab="",ylab="",col="orange",pch=21,bg="orange",cex=0.7)
print(paste("Cib Pac: ",mean(respac^2, na.rm=T)^.5,mean(respac,na.rm=T),length(respac[!is.na(respac)])))

text(-0.05,-0.7,expression(sigma),pos=2,cex=1.2)
text(-0.,-0.85,"All 0.29",pos=2)
text(-0.,-1.,"Wue 0.21",pos=2,col="green")
text(-0.,-1.15,"Kul 0.35",pos=2,col="red")
text(-0.,-1.3,"Lob 0.57",pos=2,col="blue")
text(-0.,-1.45,"Pac 0.30",pos=2,col="orange")
text(1.9,-1.45,"LA1",pos=2)

print("LA5")
#lm.deltanat.0All <- lm(c13sednatb - c13wcnatb ~ 0)
#ftd <- c13wcnatb
lm.deltanat.0All <- lm(c13sednat - c13wcnat ~ 0)
ftd <- c13wcnat
res <- resid(lm.deltanat.0All)

par(mai=c(0.8,0.9,0.1,0.2))
plot(ftd,res,xlim=xlm,ylim=ylm,col="black",pch=21,bg="black",cex=0.7,
     ylab="Residuals (\u2030)",
     xlab=expression(paste(delta^13*C[DICnat], " (\u2030)"))
    )
par(new=T)
#reswue <- res[ispecwueb]
reswue <- res[ispecwue]
#plot(ftd[ispecwueb],res[ispecwueb],xlim=xlm,ylim=ylm,xlab="",ylab="",col="green",pch=21,bg="green",cex=0.7)
plot(ftd[ispecwue],res[ispecwue],xlim=xlm,ylim=ylm,xlab="",ylab="",col="green",pch=21,bg="green",cex=0.7)
print(paste("Cib All: ",mean(res^2, na.rm=T)^.5,mean(res,na.rm=T),length(res[!is.na(res)])))
print(paste("Cib Wue: ",mean(reswue^2, na.rm=T)^.5,mean(reswue,na.rm=T),length(reswue[!is.na(reswue)])))

par(new=T)
#reskul <- res[ispeckulb]
reskul <- res[ispeckul]
#plot(ftd[ispeckulb],res[ispeckulb],xlim=xlm,ylim=ylm,xlab="",ylab="",col="red",pch=21,bg="red",cex=0.7)
plot(ftd[ispeckul],res[ispeckul],xlim=xlm,ylim=ylm,xlab="",ylab="",col="red",pch=21,bg="red",cex=0.7)
print(paste("Cib Kul: ",mean(reskul^2, na.rm=T)^.5,mean(reskul,na.rm=T),length(reskul[!is.na(reskul)])))

par(new=T)
#reslob <- res[ispeclobb]
reslob <- res[ispeclob]
#plot(ftd[ispeclobb],res[ispeclobb],xlim=xlm,ylim=ylm,xlab="",ylab="",col="blue",pch=21,bg="blue",cex=0.7)
plot(ftd[ispeclob],res[ispeclob],xlim=xlm,ylim=ylm,xlab="",ylab="",col="blue",pch=21,bg="blue",cex=0.7)
print(paste("Cib Lob: ",mean(reslob^2, na.rm=T)^.5,mean(reslob,na.rm=T),length(reslob[!is.na(reslob)])))

par(new=T)
#respac <- res[ispecpacb]
respac <- res[ispecpac]
#plot(ftd[ispecpacb],res[ispecpacb],xlim=xlm,ylim=ylm,xlab="",ylab="",col="orange",pch=21,bg="orange",cex=0.7)
plot(ftd[ispecpac],res[ispecpac],xlim=xlm,ylim=ylm,xlab="",ylab="",col="orange",pch=21,bg="orange",cex=0.7)
print(paste("Cib Pac: ",mean(respac^2, na.rm=T)^.5,mean(respac,na.rm=T),length(respac[!is.na(respac)])))

text(-0.05,-0.7,expression(sigma),pos=2,cex=1.2)
text(-0.,-0.85,"All 0.33",pos=2)
text(-0.,-1.,"Wue 0.25",pos=2,col="green")
text(-0.,-1.15,"Kul 0.36",pos=2,col="red")
text(-0.,-1.3,"Lob 0.50",pos=2,col="blue")
text(-0.,-1.45,"Pac 0.20",pos=2,col="orange")
text(1.9,-1.45,"LA5",pos=2)

dev.off()