quartz(width = 12, height = 8)
layout(matrix(c(1,2,3,4,5,6), 2, 3, byrow = TRUE))
opar<-par(mar=c(5,5,4,1)+0.1)
lo<- 50; up<-100 
ApCas<-(as.matrix(read.table("JK apatite calib_1.txt", skip=32)))[lo:up,lo:up]
ApPs<-(as.matrix(read.table("JK apatite calib_2.txt", skip=32)))[lo:up,lo:up]
#plot(c(0,2000),c(0,4000), typ='n', main = 'apatite (r), clam (gr)', cex=0.25)

ClCas<-(as.matrix(read.table("JK clam calib_1.txt", skip=32)))[lo:up,lo:up]
ClPs<-(as.matrix(read.table("JK clam calib_2.txt", skip=32)))[lo:up,lo:up]

C3Cas<-(as.matrix(read.table("JK Ca-P calib_1.txt", skip=32)))[lo:up,lo:up]
C3Ps<-(as.matrix(read.table("JK Ca-P calib_4.txt", skip=32)))[lo:up,lo:up]

#lox<- 200; upx<-280; loy<- 20; upy<-100   # was 218; upx<-278; loy<- 15; upy<-95
lox<- 25; upx<-150; loy<- 185; upy<-315   


CanFs<-(CanF<-as.matrix(read.table("JK depression b1_1.txt", skip=33)))[lox:upx,loy:upy]
lpx<-length(t(CanF)[,1]); lpy<-length(t(CanF)[1,])
CanCas<-(CanCa<-as.matrix(read.table("JK depression b1_2.txt", skip=33)))[lox:upx,loy:upy]
CanPs<-(CanP<-as.matrix(read.table("JK depression b1_3.txt", skip=33)))[lox:upx,loy:upy]
CanCls<-(CanCl<-as.matrix(read.table("JK depression b1_4.txt", skip=33)))[lox:upx,loy:upy]
CanCa5s<-(CanCa5<-as.matrix(read.table("JK depression b1_5.txt", skip=33)))[lox:upx,loy:upy]
image(t(CanP), col = grey(seq(0,1,0.01)), main="P")
lines(c(loy,upy)/lpx, c(lox,lox)/lpy,col='yellow')
lines(c(loy,upy)/lpx, c(upx,upx)/lpy,col='yellow')
lines(c(loy,loy)/lpx, c(lox,upx)/lpy,col='yellow')
lines(c(upy,upy)/lpx, c(lox,upx)/lpy,col='yellow')
image(t(CanCa), col = grey(seq(0,1,0.01)), main="Ca")
lines(c(loy,upy)/lpx, c(lox,lox)/lpy,col='yellow')
lines(c(loy,upy)/lpx, c(upx,upx)/lpy,col='yellow')
lines(c(loy,loy)/lpx, c(lox,upx)/lpy,col='yellow')
lines(c(upy,upy)/lpx, c(lox,upx)/lpy,col='yellow')
image(t(CanF), col = grey(seq(0,1,0.01)), main="F")
lines(c(loy,upy)/lpx, c(lox,lox)/lpy,col='yellow')
lines(c(loy,upy)/lpx, c(upx,upx)/lpy,col='yellow')
lines(c(loy,loy)/lpx, c(lox,upx)/lpy,col='yellow')
lines(c(upy,upy)/lpx, c(lox,upx)/lpy,col='yellow')
image(t(CanCl), col = grey(seq(0,1,0.01)), main="Cl")
lines(c(loy,upy)/lpx, c(lox,lox)/lpy,col='yellow')
lines(c(loy,upy)/lpx, c(upx,upx)/lpy,col='yellow')
lines(c(loy,loy)/lpx, c(lox,upx)/lpy,col='yellow')
lines(c(upy,upy)/lpx, c(lox,upx)/lpy,col='yellow')
image(t(CanCa5), col = grey(seq(0,1,0.01)), main="Ca")
lines(c(loy,upy)/lpx, c(lox,lox)/lpy,col='yellow')
lines(c(loy,upy)/lpx, c(upx,upx)/lpy,col='yellow')
lines(c(loy,loy)/lpx, c(lox,upx)/lpy,col='yellow')
lines(c(upy,upy)/lpx, c(lox,upx)/lpy,col='yellow')
plot(c(0,1500),c(0,3500), typ='n', 
                        main = 'apatite (r), clam (gr), CaPO4 (blu)',
                        xlab='[P]' , ylab='[Ca]',
                        cex.lab=1.5,
                        cex=0.25)
points(ApPs,ApCas, typ='p', col='red', cex=0.25)
points(ClPs,ClCas, typ='p', col='green', cex=0.25)
points(C3Ps,C3Cas, typ='p', col='blue', cex=0.25)
points(CanPs,CanCas, typ='p', col='purple', cex=0.25)
mtext(paste("xlo,xhi=",lox,", ", upx," /  ylo,yhi=", loy,", ", upy,sep=''), line=-2)

browser()
quartz(width = 12, height = 4)
layout(matrix(c(1,2,3,4), 1, 4, byrow = TRUE))
opar<-par(mar=c(1,2,3,2)+0.1)

image(t(CanPs), col = grey(seq(0,1,0.01)), main="P 2")
image(t(CanCas), col = grey(seq(0,1,0.01)), main="Ca 3")
image(t(CanFs), col = grey(seq(0,1,0.01)), main="F 1")
image(t(CanCls), col = grey(seq(0,1,0.01)), main="Cl 4")
#image(t(CanCa5s), col = grey(seq(0,1,0.01)), main="Ca 5")