# Setup graphics windows() par(mfrow=c(2,2)) par(mar=c(4.5,4.5,1,1)) par(pty="s") srlab <- expression(""^87*Sr/""^86*Sr) ndlab <- expression(""^143*Nd/""^144*Nd) lwd <- 2 my.palette <- colorRampPalette(c("black", "darkgreen", "red"),space = "rgb") # a to b calculation cx1<- 150; ix1 <- 0.715 # end member 1 Sr composition cx2 <- 600; ix2 <- 0.703 # end member 2 Sr composition f1 <- seq(0,1,by=0.05) # proportions of "1" in the mixture cmx <- cx1*f1+(1-f1)*cx2 # Sr concentration in the mix imx <- ix1*cx1*f1/cmx+ix2*cx2*(1-f1)/cmx # 87Sr/86Sr in the mix # a plot(cmx,imx,xlim=c(100,700),ylim=c(0.702,0.7155),xlab="Sr (ppm)",ylab=srlab,type="b",col=my.palette(length(f1)),pch=19,lwd=lwd,xaxs="i",yaxs="i") # End members text(cx1,ix1,"1",font=2,pos=2,col="red") text(cx2,ix2,"2",font=2,pos=4) # b plot(1/cmx,imx,,xlim=c(0.001,0.0077),ylim=c(0.702,0.7155),xlab="1/Sr (ppm)",ylab=srlab,type="b",col=my.palette(length(f1)),pch=19,lwd=lwd,xaxs="i",yaxs="i") # End members text(1/cx1,ix1,"1",font=2,pos=4,col="red") text(1/cx2,ix2,"2",font=2,pos=2) # c calculation cy1 <- 20; cy2 <- 2 # Nd in "1" and "2" iy1 <- 0.511; iy2 <- 0.513 # 143Nd/144Nd, "1" and "2" cmy <- cy1*f1+(1-f1)*cy2 # Nd conc., mix imy <- iy1*cy1*f1/cmy+iy2*cy2*(1-f1)/cmy # 143Nd/144Nd mix # Plot the mixing hyperbola plot(imx,imy,xlab=srlab,ylab=ndlab,type="b",xlim=c(0.702,0.716),ylim=c(0.5109,0.5132),pch=19,col=my.palette(length(f1)),lwd=lwd,xaxs="i",yaxs="i") # End members text(ix1,iy1,"1",font=2,pos=4,col="red") text(ix2,iy2,"2",font=2,pos=2) # Curvature and asymptotes alpha <- (cx2/cy2)/(cx1/cy1) # calc. alpha x0 <- (ix1-alpha*ix2)/(1-alpha) # calc. asymptotes y0 <- (iy2-alpha*iy1)/(1-alpha) abline(v=x0,lty="dashed") # draw asymptotes abline(h=y0,lty="dashed") # d # Input parameters x1<-0.715; x2<-0.703 y1<-0.511; y2<-0.513 alpha <- c(5,2,1,0.5,0.2) Sr1 <- 150 Sr2 <- 600 Nd1 <- 20 Nd2 <- Sr2*Nd1/Sr1/alpha # Parameters to the hyperbola A <- y2*Nd2*Sr1-y1*Nd1*Sr2 B <- Nd1*Sr2-Nd2*Sr1 C <- x1*Nd2*Sr1-x2*Nd1*Sr2 D <- y1*x2*Nd1*Sr2-y2*x1*Nd2*Sr1 x <- seq(x1,x2,length=11) col <- rev(my.palette(length(x))) plot(0,0,xlim=c(0.702,0.716),ylim=c(0.5109,0.5132),type="n",xlab=srlab,ylab=ndlab,add=TRUE,lwd=lwd,xaxs="i",yaxs="i") # Draw hyperbolae, one at time ee<-lapply(1:length(alpha),function(i){ eq <- eval(parse(text=paste("function(x)(",-D[i],"-",A[i],"*x)/(",B[i],"*x+",C[i],")",sep=""))) #curve(eq,from=x1,to=x2,add=TRUE,col="black",lwd=1) y <- eq(x) points(x,y,pch=19,cex=0,type="b",col="black",lwd=lwd) points(x,y,pch=19,col=col,type="p") text(x[5],y[5]+0.0001,alpha[i],font=3) invisible() }) # Asymptotes and curvature x0 <- (x1-alpha*x2)/(1-alpha) y0 <- (y2-alpha*y1)/(1-alpha) q <- x0*y0+(alpha*x2*y1-x1*y2)/(1-alpha) # End members text(x1,y1,"1",font=2,pos=4,col="red") text(x2,y2,"2",font=2,pos=2) dev.copy2pdf(file="PDF output from R/Fig 3 Mix.pdf",colormodel="rgb")