# Read in callibration data
layout(matrix(c(1,2), 1, 1, byrow = TRUE))
opar<- par(mar=c(5, 3, 4, 3) + 0.1)
# Adjust skip to get Conc and V(Conc) headers.

Z<- read.table("EXPOR016.TXT", header=TRUE, skip=7, nrows=1, sep=',')

# Adjust multiplier to make 
iConc<-7:9
iVconc<-10:12
Conc<-as.real(Z[iConc])/1000
Vconc<-as.real(Z[iVconc])
Lconc<-log10(Conc)
ELconc<-predict(lm(Lconc~Vconc))
plot(Vconc, ELconc, ty='b', main="-pH vs mV CALLIBRATION")
points(Vconc, Lconc,col='red')
points(90, -5.25, col='red')
text(93, -5.25, "Obs[-pH]", col='red', pos=4)
points(90, -5.35, col='black')
text(93, -5.35, "E[-pH]", col='black', pos=4)
points(90, -5.45, col='blue', cex=2)
text(93, -5.45, "New E[-pH]", col='blue', pos=4)

NewVC <- data.frame(Vconc = c(100, 120, 140, 160, 180, 200))
points(t(NewVC), predict(lm(Lconc~Vconc), newdata=NewVC), cex=2, col='blue')
browser()

layout(matrix(c(1:8), 2, 4, byrow = TRUE))

filnam<-"EXPOR015.TXT"
X<-read.table(filnam, header=TRUE, skip=0, nrows=-1, sep=',')
print(names(X))
attach(X)
minAll<-min(min(AvgO.XuV),min(AvgO.YuV),min(AvgO.ZuV))
maxAll<-max(max(AvgO.XuV),max(AvgO.YuV),max(AvgO.ZuV))
Ypos<- ProbeY- sum(ProbeY)/length(ProbeY)
Xpos<- ProbeX- sum(ProbeX)/length(ProbeY)
Zpos<- ProbeZ- sum(ProbeZ)/length(ProbeY)
plot(c(min(Ypos),max(Ypos)), c(minAll,maxAll), 
      ty='n', 
      main=paste(filnam, 'Obs[uVdiff]'))
      mtext("Observed uV difference (xyz) = (rgb)", line=-1.5)
lines(Ypos, AvgO.XuV, col="red")
lines(Ypos, AvgO.YuV, col="green")
lines(Ypos, AvgO.ZuV, col="blue")

minAll<-min(Origin.1.mV)
maxAll<-max(Origin.1.mV)
# plot(c(min(Ypos),max(Ypos)), c(minAll,maxAll), 
#       ty='n', 
#       main=paste(filnam, 'Obs[mV]'))
#       mtext("Observed mV at probe origin", line=-1.5)
# lines(Ypos, Origin.1.mV, col="black")

NewV0 <- data.frame(Vconc = Origin.1.mV)
plot(Ypos, predict(lm(Lconc~Vconc), newdata=NewV0), 
      ty='b', 
      main=paste(filnam, 'E[-pH]'))
      mtext("E[-pH] = by callibration", line=-1.5)

Oconc<- 10^(predict(lm(Lconc~Vconc), newdata=NewV0))
plot(Ypos, Oconc, 
      ty='b', 
      main=paste(filnam, 'E[H+]'))
      mtext("E[H+] = 10^(E[-pH]) at probe origin", line=-1.5)

Xconc<- 10^(predict(lm(Lconc~Vconc), newdata=NewV0+AvgO.XuV/1000))
Yconc<- 10^(predict(lm(Lconc~Vconc), newdata=NewV0+AvgO.YuV/1000))
Zconc<- 10^(predict(lm(Lconc~Vconc), newdata=NewV0+AvgO.ZuV/1000))
dXH<- Xconc-Oconc      
dYH<- Yconc-Oconc      
dZH<- Zconc-Oconc 
minAll<-min(min(dXH),min(dYH),min(dZH))
maxAll<-max(max(dXH),max(dYH),max(dZH))
dHtot<-(dXH^2 + dYH^2 + dZH^2)^0.5
Hflux<- dHtot*10^19     
plot(c(min(Ypos),max(Ypos)), c(minAll,maxAll), 
      ty='n', 
      main=paste(filnam, '[H+]dxyz'))
      mtext("[H+]dxyz", line=-1.5, col='red')
lines(Ypos, dXH, col="red", ty='b')
lines(Ypos, dYH, col="green", ty='b')
lines(Ypos, dZH, col="blue", ty='b')

detach(X)
require(scatterplot3d)
 # mul<-2*10^8
 mul<-600/maxAll
 
 Miny<- min(c(Ypos,Ypos+mul*dYH))
 Minz<- min(c(Zpos,Zpos+mul*dZH))
 Minx<- min(c(Xpos,Xpos+mul*dXH))
 Maxy<- max(c(Ypos,Ypos+mul*dYH))
 Maxz<- max(c(Zpos,Zpos+mul*dZH))
 Maxx<- max(c(Xpos,Xpos+mul*dXH))

 
 
 # Left Scat 3D                
s3d<-scatterplot3d(c(Miny,Maxy), c(Minz,Maxz), c(Minx,Maxx), highlight.3d=TRUE,
      col.axis="blue", col.grid="lightblue",
      main="scatterplot3d - L", pch=16,
      cex.symbols=2,
      angle= 43,
      typ='n',
      axis=FALSE,
      mar=c(0.5,0.5,0.5,0.5),
      scale.y=0.25)
for(i in 1:121) {
s3d$points3d(Ypos[i],
             Zpos[i],
             Xpos[i],
      col="black", type="b", pch=16)
                 }
for(i in 1:121) {
s3d$points3d(c(Ypos[i],Ypos[i]+mul*dYH[i]),
             c(Zpos[i],Zpos[i]+mul*dZH[i]),
             c(Xpos[i],Xpos[i]+mul*dXH[i]),
      col="blue", type="l", pch=16)
                 }
for(i in 1:121) {
s3d$points3d(Ypos[i]+mul*dYH[i],
             Zpos[i]+mul*dZH[i],
             Xpos[i]+mul*dXH[i],
      col="red", type="b", pch=16)
                 }
 # Right Scat 3D 
               
s3d<-scatterplot3d(c(Miny,Maxy), c(Minz,Maxz), c(Minx,Maxx), highlight.3d=TRUE,
      col.axis="blue", col.grid="lightblue",
      main="scatterplot3d - R", pch=16,
      cex.symbols=2,
      angle= 37,
      typ='n',
      axis=FALSE,
      mar=c(0.5,0.5,0.5,0.5),
      scale.y=0.25) # mar(b,l,t,r)
for(i in 1:121) {
s3d$points3d(Ypos[i],
             Zpos[i],
             Xpos[i],
      col="black", type="b", pch=16)
                 }
for(i in 1:121) {
s3d$points3d(c(Ypos[i],Ypos[i]+mul*dYH[i]),
             c(Zpos[i],Zpos[i]+mul*dZH[i]),
             c(Xpos[i],Xpos[i]+mul*dXH[i]),
      col="blue", type="l", pch=16)
                 }
for(i in 1:121) {
s3d$points3d(Ypos[i]+mul*dYH[i],
             Zpos[i]+mul*dZH[i],
             Xpos[i]+mul*dXH[i],
      col="red", type="b", pch=16)
                 }
text(125,10, "Y")
opar<- par(mar=c(5, 3, 4, 3) + 0.1)                 
Yax<- seq(-150, 150, 30)
Xax<-Yax
Flux<-matrix(Hflux, 11,11)      
persp(Flux, theta = 330, phi = 30, expand = 1.2, col = "lightblue", zlab="H-flux", xlab="Z-axis", ylab="Y-axis")
persp(Flux, theta = 333, phi = 30, expand = 1.2, col = "lightblue", zlab="H-flux", xlab="Z-axis", ylab="Y-axis")