# Transform into script to plot data from a SVD
inplot<- function(ix,iy,cx,cy, rsq) {points(t(CanPs)[ix,iy],t(CanCas)[ix,iy], typ='p', col='purple', cex=0.25)
	                                }
lo<- 50; up<-250 
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]
loy<-87; upy<-162; lox<- 325; upx<-425


CanCas<-(CanCa<-as.matrix(read.table("JK elongate tube 1_2.txt", skip=32)))[loy:upy,lox:upx]
CanCa2s<-(CanCa2<-as.matrix(read.table("JK elongate tube 1_5.txt", skip=32)))[loy:upy,lox:upx]
lpx<-length(t(CanCa)[,1])
lpy<-length(t(CanCa)[1,])
CanPs<-(CanP<-as.matrix(read.table("JK elongate tube 1_3.txt", skip=32)))[loy:upy,lox:upx]
CanFs<-(CanF<-as.matrix(read.table("JK elongate tube 1_1.txt", skip=32)))[loy:upy,lox:upx]
CanCls<-(CanCl<-as.matrix(read.table("JK elongate tube 1_4.txt", skip=32)))[loy:upy,lox:upx]
IonM<-matrix(0,Lem<-length(CanCas),5)
IonM[,1]<- CanCas[1:Lem]
IonM[,2]<- CanCa2s[1:Lem]
IonM[,3]<- CanPs[1:Lem]
IonM[,4]<- CanCls[1:Lem]
IonM[,5]<- CanFs[1:Lem]
SVD5<-svd(IonM)
nrows<-length(CanCas[,1])
ncols<-length(CanCas[1,])

VV<-SVD5$v
# VV[,3]<- -VV[,3]
VD<-SVD5$d
Nsvd<-length(VD)
VDs<- sum(SVD5$d)
VDp<-100*VD/VDs
Xwid<- length(CanCas[1,])
Ywid<- length(CanCas[,1])
cat('Top',Nsvd,'SVDs percent:',VDp,'\n')
Ev<-t(VV)%*%t(IonM)
Ilab<- c("Ca+Ca", "Ca-Ca","P", "Cl", "F")
for (i in 1:Nsvd){
quartz(width=1+ 3*Xwid/Ywid, height=3)
opar<- par(mar=c(1,1,1,4))
PCV<-matrix(Ev[i,], nrows,ncols, byrow=FALSE)  # [nrows:1,]  # [nrows:1,ncols:1]
filled.contour(1:ncols, 1:nrows, t(PCV), main=paste("PCV",i, Ilab[i]), color = topo.colors)
	            }

browser()

quartz(width = 9, height = 7)     # quartz(width = 8, height = 7.85)
layout(matrix(c(1,1,1,1,2,2,2,2,
                1,1,1,1,2,2,2,2,
                1,1,1,1,2,2,2,2,
                7,7,3,3,6,6,6,6,
                7,7,3,3,6,6,6,6,
                4,4,5,5,6,6,6,6,
                4,4,5,5,6,6,6,6), 7, 8, byrow = TRUE))

#opar<- par(mar=c(2,1,1,1))
opar<- par(mar=c(0.2,0.2,0.2,0.5))

image(t(CanP), col = grey(seq(0,1,0.01)), xaxt='n', yaxt='n')
lines(c(lox,lox)/lpx,c(loy,upy)/lpy, col='yellow')
lines(c(upx,upx)/lpx,c(loy,upy)/lpy, col='yellow')
lines(c(lox,upx)/lpx,c(loy,loy)/lpy, col='yellow')
lines(c(lox,upx)/lpx,c(upy,upy)/lpy, col='yellow')
image(t(CanCa), col = grey(seq(0,1,0.01)), xaxt='n', yaxt='n')
lines(c(lox,lox)/lpx,c(loy,upy)/lpy, col='yellow')
lines(c(upx,upx)/lpx,c(loy,upy)/lpy, col='yellow')
lines(c(lox,upx)/lpx,c(loy,loy)/lpy, col='yellow')
lines(c(lox,upx)/lpx,c(upy,upy)/lpy, col='yellow')

opar<- par(mar=c(0.2,0.2,0.2,0.5))
image(t(CanCas), col = grey(seq(0,1,0.01)), xaxt='n', yaxt='n')
lox<- 1; upx<-lpx<-length(t(CanCa)[,1])
loy<- 1; upy<-lpy<-length(t(CanCa)[1,])
PCV2<-matrix(Ev[3,], nrows,ncols, byrow=FALSE)
image(t(CanPs), col = grey(seq(0,1,0.01)), xaxt='n', yaxt='n')
image(t(PCV2), col = grey(seq(0,1,0.01)), xaxt='n', yaxt='n')

opar<- par(mar=c(2,1.5,1,1))
plot(c(0,1600),c(0,3800), typ='n', main = 'apatite (r), clam (g), Ca Phosphate (b)', cex=0.25)
points(ApPs,ApCas, typ='p', col='red', cex=0.25)
MApPs<- mean(ApPs); MApCas<- mean(ApCas)
points(c(0,MApPs), c(0,MApCas))
lines(c(0,MApPs), c(0,MApCas), col='red',  lty=4, lwd=3)
text(MApPs,MApCas, 1.67, col='red', pos=3, cex=2, offset=2.5)
Rdydx<- MApCas/MApPs

points(ClPs,ClCas, typ='p', col='green', cex=0.25)

points(C3Ps,C3Cas, typ='p', col='blue', cex=0.25)
MC3Ps<- mean(C3Ps); MC3Cas<- mean(C3Cas)
Bdydx<- MC3Cas/MC3Ps
points(c(0,MC3Ps), c(0,MC3Cas))
lines(c(0,MC3Ps), c(0,MC3Cas), col='blue',  lty=4, lwd=3)
text(MC3Ps,MC3Cas, round(1.67*Bdydx/Rdydx,2), col='blue', pos=3, offset=2.5, cex=2)
rat<- 2.400 # 2.09996 produces a 3.50 line for CAp 
s<-0.9
points(s*c(0,MApPs/rat), s*c(0,MApCas), col='black')
text(s*MApPs/rat, s*MApCas, round(1.67*rat,2), col='black', pos=3, offset=1, cex=2)
lox<- 1; upx<-lpx<-length(t(CanCas)[,1])
loy<- 1; upy<-lpy<-length(t(CanCas)[1,])
Cx<-(lox+upx)/(2*lpx); Cy<- (loy+upy)/(2*lpy); Rx<- (upx-lox)/(2*lpx)
for (i in 1:lpx) {
	for(j in 1:lpy){if (PCV2[j,i] < -25) {inplot(i,j,Cx,Cy, Rx^2)}		 
		           }
	             }
lines(s*c(0,MApPs/rat), s*c(0,MApCas), lty=2, lwd=3, col='black')
rat<- 4.191617   # 1.199998 for 2.0;  1.6 for 2.67;  2.399995 for 4.0 # 4.191617 for a 7.0 line for CAp
s<-0.81
text(s*MApPs/rat, s*MApCas, round(1.67*rat,2), col='black', pos=3, offset=1, cex=2)
lines(s*c(0,MApPs/rat), s*c(0,MApCas), lty=2, lwd=3, col='black')

# browser()

SVi<- rep(1:Nsvd,Nsvd)	            
VVi<- rep(VV[1:(Nsvd*Nsvd)],1)
opar<- par(mar=c(1,1,2,1))
plot(SVi,VVi, typ='n', main=paste("PCs"))
for (i in 1:Nsvd) {lines(SVi[1:Nsvd], VVi[Nsvd*(i-1)+1:Nsvd], col=rainbow(5)[i], lwd=6-i)
	            text(SVi[i], VVi[Nsvd*(i-1)+i], Ilab[i], col=rainbow(5)[i], pos=3)
	            }            
lines(SVi[1:Nsvd],rep(0,Nsvd), lty=3, lwd=5)