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  
lox<- 80; upx<-160; loy<- 270; upy<-390   

CanFs<-(CanF<-as.matrix(read.table("JK bristle 1_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 bristle 1_2.txt", skip=33)))[lox:upx,loy:upy]
CanPs<-(CanP<-as.matrix(read.table("JK bristle 1_3.txt", skip=33)))[lox:upx,loy:upy]
CanCls<-(CanCl<-as.matrix(read.table("JK bristle 1_4.txt", skip=33)))[lox:upx,loy:upy]
CanCa5s<-(CanCa5<-as.matrix(read.table("JK bristle 1_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)

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

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")