# Uses wall width filnam.log files and filnam.log#.txt 
# files and fits gausian distributions to the edges 
# Computes and plots the tip velocity on the same plot with wall thickness.
# Enter variable parameters here:
npeak<-0
Perc<- 0.85
dx<-25

# End of parameter entry
cat("This R-script is meant to process the output file of Wall_1_prep.R which sets initial best guesses for\n")
cat(" peak locations.  Abort if you have not used Wall_1_prep.R on the log file yet!\n")
cat("The F)latten key function has been eliminated and is now automated. The function keys have been commented out.\n")
cat("The i)nherit key function works well to set up the next contour using the last contour's SDs and peak adjustments.\n")

# Function Deffinitions -----------------------------

ask<- function(text) {
cat(text,": ", sep="")
indata<-readLines(con=stdin(),n=1,ok=TRUE)
                     }
                     
Graphic<- function(A=Out, FIL=filnam, IREP=irep) {
par(mar=c(4, 4, 4, 4) + 0.1)
layout(matrix(c(1:10), 1, 1, byrow = TRUE))
NF<-length(A[,1])
plot(A[,1], A[,3]-A[,2], main = paste(FIL,"DIFF", "rep =",IREP), ty='l')
lines(A[,1], A[,3]-A[,2], typ='b', col='orange', lwd=2)
minOut3m2<- min(A[,3]-A[,2])
maxOut3m2<- max(A[,3]-A[,2])
dOut3m2<-maxOut3m2-minOut3m2
minV<-min(A[2:NF,2]-A[1:(NF-1),2])
maxV<-max(A[2:NF,2]-A[1:(NF-1),2])
dV<-maxV-minV
V<- minOut3m2 + dOut3m2*(A[2:NF,2]-A[1:(NF-1),2] -minV)/dV
lines(A[2:NF,1]-0.5, V, col='blue', lwd=2)
ind<- seq(0,length(A[,1]),10)
ind[1]<-1
OutX<-A[ind,1]
OutY<- A[ind,3]- A[ind,2]
text(OutX, OutY, labels= OutX, pos=3)
mtext(" Blue = Velocity", side = 3, adj=0, line= -2, col = 'blue', cex=2)
mtext("Orange = Wall Thickness ", side = 3, adj=1, line= -2, col = 'orange', cex=2)
browser()
par(mar=c(2.5, 0.1, 2, 0.1) + 0.1)
layout(matrix(c(1:10), 1, 10, byrow = TRUE))
}



MovePeak<- function(X,iA, A=Out, da=dx, gauss=Gauss, FIL=filnam) {
layout(matrix(c(1,2,3), 1, 3, byrow = TRUE))
while(iA>0) {
        imax<- A[iA, 9]; jmax<- A[iA, 11];
        N<-length(A[,1])
        Nl<- round(length(X[,1])/N)
        Y<-as.matrix(X[(4:Nl)+(iA-1)*Nl,])
        Xa<-spline(as.real(Y[,1]), as.real(Y[,3]))
        minx<- max(c(1, imax-da));      maxx<- min(c(length(Xa$y), imax+da))
        Xax<- Xa$x[minx:maxx];          Xay<- Xa$y[minx:maxx]
        imax<-which.max(Xay);           imin<-which.min(Xay)
        sig<- mi<- 1; dm<- 0.01; mag<-dm + 1;  F<- 1; D<-1
        if (A[iA, 8] ==1) { F<- A[iA, 4]; D<- A[iA, 5]}
        sig<-sig1<- A[iA, 10];  u1<- A[iA, 2]; sig2<- A[iA, 12];  u2<- A[iA, 3]                        
        dx<- ds<- 0.004; di<- 1    # Move peaks
           while (!ds==0) {  
           	                 plot(Xax, Xay, ty='p', lty= 1, main=paste(FIL,iA))
           	                 
                             lines(c(u1, u1), c(Xay[imin], Xay[imax]), col='black', lty=3)
                             lines(c(u2, u2), c(Xay[imin], Xay[imax]), col='red', lty=2)
                             Ey1<- dnorm(Xax, u1, sig1); Ey2<- dnorm(Xax, u2, sig2)
                             Ey0<- Ey1 + D*Ey2; Mult<- sum(Xay)/sum(Ey0)
                             F<- (Xay[imax] - Xay[imin])/(Mult*Ey0[imax])
                             Ey<- Xay[imin] + Mult*Ey0*F
                             Ey1<- Xay[imin] + Mult*Ey1*F
                             Ey2<- Xay[imin] + Mult*D*Ey2*F
                             lines(Xax, Ey1, col='green', lty=2)
                             lines(Xax, Ey2, col='green', lty=2)
                             lines(Xax, Ey, col='green')
                             mtext(paste(" (u1,s1)=",round(u1,3), round(sig1,3)) , line= 0.5, adj=0, cex=1, col = 'blue')                              
                             mtext(paste(" (F,D)",round(F,3), round(D,3)) , line=-1.5, adj=0, cex=1, col = 'blue')                              
                             mtext(paste(" mag=",round(mag,3)) , line=-3, adj=0, cex=1, col = 'blue')                              
                             mtext(paste(" dx=",round(dx,4)) , line=-4.5, adj=0, cex=1, col = 'blue')                              
                             mtext(paste("(u2,s2)=",round(u2,3), round(sig2,3)) , line=0.5, adj=1, cex=1, col = 'blue')
                             mtext(" Active " , line=-9.5, adj=mi-1, cex=1, col = 'red') 
       cat(" prox/dist || -/+dx | stop || x2 / div 2 || browse | mag/min | dMag | Diff | Inherit |abort |\n")
 char<-ask("  <  / >   ||  -/+  |   =  ||  2 /   5   ||   b    |  S   s  |  p/P |  d/D |    i/I  |  *   |")
                            if (char=='b'){ browser()}
                            if (char=='-'){ di<- -1;  cat("minus",dx,". \n")
                                            if (mi==1) {u1<-u1+di*dx; sig<- sig1} else {u2<-u2+di*dx; sig<- sig2}
                                          }
                            if (char=='+'){ di<-  1;  cat("plus",dx,". \n")
                                            if (mi==1) {u1<-u1+di*dx; sig<- sig1} else {u2<-u2+di*dx; sig<- sig2}
                                          }
                            if (char=='='){ di<-  1;  cat("plus",dx,". \n")
                                            if (mi==1) {u1<-u1+di*dx; sig<- sig1} else {u2<-u2+di*dx; sig<- sig2}
                                          }
                            if (char=='2'){ dx<- 2*dx;  cat("dx doubled to",dx,".\n")}
                            if (char=='5'){ dx<- dx/2;  cat("dx halved to",dx,".\n")}
                            if (char=='*'){ ds<-0; cat("Exit this frame! \n")}
                            if (char=='<'){ mi<-1; sig<- sig1} # dx<- ds; removed
                            if (char=='>'){ mi<-2; sig<- sig2} # dx<- ds; removed
                            if (char=='S'){ sig<- mag*sig; if (mi==1) sig1<-sig else sig2<- sig }
                            if (char=='s'){ sig<- sig/mag; if (mi==1) sig1<-sig else sig2<- sig }
                            if (char=='P'){dm<- 2* dm; mag<- 1+dm}
                            if (char=='p'){dm<- dm/2; mag<- 1+dm}
                            # if (char=='F'){F<- mag*F}
                            # if (char=='f'){F<- F/mag} # F<- (Xay[imax] - Xay[imin])/(Ey[imax] - Ey[imin])
                            # if (char=='g') F<- (Xay[imax] - Xay[imin])/(Mult*(Ey0[imax] - Ey0[imin]))                          
                            if (char=='D'){D<- mag*D}
                            if (char=='d'){D<- D/mag}                            
                            if (char=='i'){load(file="inheritance"); sig<- sig1
                            	           F<- (Xay[imax] - Xay[imin])/(Mult*(Ey0[imax] - Ey0[imin]))}
                            if (char=='I'){load(file="inheritance"); sig<- sig1
                            	           F<- (Xay[imax] - Xay[imin])/(Mult*(Ey0[imax] - Ey0[imin]))}                            
                            }    # end of while !ds==0   
save(F,D,dm,sig1,sig2, file="inheritance")
 A[iA,2]<- u1; A[iA,3]<- u2; A[iA,4]<- F; A[iA,5]<- D; A[iA,10]<- sig1; A[iA,12]<- sig2; A[iA,8]<- 1                            
# A[iA,]<- as.real(c(A[iA,1], u1, u2, F, D, Diff, 3, 1, minx-1+imax, sig1, jmax, sig2 ))

save(A, file="Asaved")
iA<- as.real(ask("New frame number? {0 exits} "))
} # end while iA > 0            
layout(matrix(c(1:10), 1, 10, byrow = TRUE)) 
A }
                                       
 # End of Function Deffinitions #----------------------------------------
         
filnam<-"junk"
notuseit<- TRUE
if(fexist<- file.exists("LastFile")) {load("LastFile"); cat("Last file name used =", filnam,".")
	                         notuseit<-('n'== ask("Use it?(y/n)")) }
	                         if(notuseit) {
                                          print(list.files("./","log"))
                                          filnam<- print(ask("Enter log file name"))
	                                      } 	                         
XX<-readLines(filnam)
save(filnam, file="LastFile")

if(cexist<- file.exists(Cfile<-paste(filnam,".comment.txt",sep=''))) {Ctext<-readLines(Cfile); 
	        cat("\nComment file exists:", Cfile,".\n\n")
	        print(Ctext)
	                }
	                         	                         


Nframes<-length(grep("Image Plane", XX))   # Automatic frame counting

XX<-read.csv(filnam, header=FALSE)
Nlines<- round(length(XX[,1])/Nframes)

cat("This file has", Nframes, "Image Frames and has", Nlines, "pixel width per Frame.\n")
iOff<- 0
ibrowse<-0
# Three types of startups could be possible:
# (1) Startup with a new log data set is not possible because this needs a txt output file from the prior processing by Wall_all.R 
# (2) Startup with a previous Out file saved as a textfile -or- 
# (3) Restart after a crash using current Out matrix
# NB- You can use the same or smaller output index if you wish to overwrite the output with that index.

# {
cat("Available files: \n")
print(list.files("./",filnam))
iin<- as.real(ask("Enter the input index to load"))
if('n'== ask("Use the Out matrix in memory?")) { Out<-read.table(file=paste(filnam,iin,".txt", sep=''))} 
# Always give opportunity to set a new output suffix number                                             
irep<- as.real(ask("Enter the output index to append"))
# irep is the numeric suffix used to label the Output of the script.
# irepeat choses whether to analyze the data in another cycle
# Nframes are the number of image frames analyzed
# secondset

irepeat<-1
# Analysis Loops 
 while( irepeat > 0) {
par(mar=c(2.5, 0.1, 2, 0.1) + 0.1)
layout(matrix(c(1:10), 1, 10, byrow = TRUE)) 

i<-0
while( i < Nframes+10) {
i<-i +1
if( i<Nframes+1) {
secondset<- FALSE; nosecond<- TRUE
YY<-as.matrix(XX[(4:Nlines)+(i-1)*Nlines,])
Xo<-spline(as.real(YY[,1]), as.real(YY[,3]))
Yo<-spline(as.real(YY[,1]), as.real(YY[,2]))
Yos<-Yo$y
Gauss<- Out[i,8]
npeak<-Out[i,7]
if(abs(Out[i,7])< 2) {imax<- Out[i, 9]; jmax<- Out[i, 9]; dij<- jmax-imax}                
if(abs(Out[i,7])==3) {secondset<-TRUE; nosecond<- FALSE; imax<- Out[i, 9]; jmax<- Out[i, 11]; dij<- jmax-imax}
if(abs(Out[i,7])==2) {secondset<-TRUE; nosecond<- FALSE; imax<- Out[i, 9]; jmax<- Out[i, 11]; dij<- jmax-imax} 
minx<- max(c(1, sx<-imax-dx));      maxx<- min(c(length(Xo$y), imax+dx))
Xox<- Xo$x[minx:maxx];          Xoy<- Xo$y[minx:maxx]
Yoy<- Yos[minx:maxx];           XoyCrit<-Perc*max(Xoy)

if(sx< 1) {imax<- dx+1+sx; jmax<- imax + dij} else { imax<- dx+1; jmax<- imax + dij } # Works for all Out[,7] types

imin<-which.min(Xoy)  
if(!Gauss) {plot(Xox, Xoy, ty='l', lty= 1, main=i, lwd=1.5)
            lines(c(Xox[imax], Xox[imax]), c(Xoy[imin], Xoy[imax]), col='black', lty=3)
            lines(c(Xox[jmax], Xox[jmax]), c(Xoy[imin], Xoy[jmax]), col='green', lty=4)
           }
if (Gauss) { F<- Out[i, 4]; D<- Out[i, 5]        
            sig1<- Out[i, 10];  u1<- Out[i, 2]; sig2<- Out[i, 12];  u2<- Out[i, 3]                        
            plot(Xox, Xoy, ty='p', lty= 1, main=i)
            lines(c(u1, u1), c(Xoy[imin], Xoy[imax]), col='blue', lty=3)
            lines(c(u2, u2), c(Xoy[imin], Xoy[imax]), col='red', lty=3)
            Ey1<- dnorm(Xox, u1, sig1); Ey2<- dnorm(Xox, u2, sig2)
            Ey<- Ey1 + D*Ey2; Mult<- sum(Xoy)/sum(Ey)
            Ey<- Xoy[imin] + Mult*Ey*F
            Ey1<- Xoy[imin] + Mult*Ey1*F
            Ey2<- Xoy[imin] + Mult*D*Ey2*F
            lines(Xox, Ey, col='black', lty=1, lwd=1)
            lines(Xox, Ey1, col='blue', lty=2, lwd=2)
            lines(Xox, Ey2, col='red', lty=2, lwd=2)
            mtext(paste("u,s1=",round(u1,2), round(sig1,2)) , line=-1.15, adj=0, cex=0.8, col = 'blue')                              
            mtext(paste("u,s2=",round(u2,2), round(sig2,2)) , line=-2.25, adj=1, cex=0.8, col = 'red')
            }

lines(Xox[c(1,maxx-minx)],c(XoyCrit, XoyCrit), col="orange", lty=2, lwd=2)
# maxflat<-FALSE; #  maxi<-length(Xox); # lowi<-imax; hii<-imax  

mtext(npeak, line=-6, adj=0, cex=1.5, col = 'blue') 
if(npeak==1) Out[i,11] <- Out[i,9]
cat(i, "Max ix,y:",Out[i,9], Out[i,11],
      "| type:", Out[i,7], Out[i,2], Out[i,3], "dif(X)=", Difx<- abs(Out[i,3]-Out[i,2]), "\n"); 
     # Out[i,2:6]<- c(Xox[imax], Xox[jmax], Difx, Dify, Diff)  # defeats Objective
        } # end of "if i < Nframes" allowing truncating profiles at end of frames 
         
 if (i==ibrowse) {browser(); ibrowse<-0}
 if (((i+iOff)/10)== round((i+iOff)/10)) {mtext(filnam, line=-1, adj=0, cex=1, col = 'red', side=1)
            cat(" back | fwd | 98% dx | 102% dx | tog-5 | modify | goto | Graphic | browse |  wid | nar | abort | load A |\n")
      char<-ask("  <   |  >  |   -    |    +    |   5   |   m    |   g  |    G    |   b    |   w  |  n  |   *   |    a   |")
                            if (i>10) { if (char=='<'){i<-i-20}                                      }
                            if (char=='a'){load(file="Asaved"); Out<- A;i<-i-10}
                            if (char=='b'){ibrowse<- as.real(ask("Enter frame number to browse at"));i<-i-10}
                            if (char=='g'){i<- as.real(ask("Goto frame number "));i<-i-5; cat("Go to frame",i,"as center frame. \n")}
                            if (char=='G'){Graphic(); i<-i-10}
                            if (char=='-'){ Perc<-0.98*Perc; i<-i-10; cat("Criterion increased to",Perc,"of peak. \n")}
                            if (char=='5'){ if(iOff==0) {iOff<-5} else {iOff<-0}; i<-i-5}
                            if (char=='+'){ Perc<-Perc/0.98; i<-i-10; cat("Criterion decreased to",Perc,"of peak. \n")}
                            if (char=='w'){ dx<-dx+1; i<-i-10; cat("dx widened to",dx,"\n")}
                            if (char=='n'){ dx<-dx-1; i<-i-10; cat("dx narrowed to",dx,"\n")}
                            if (char=='*'){ stop("Abort this program")}
                            if (char=='m'){Out<- MovePeak(XX, fn<-as.real(ask("Enter frame number to move peak")));i<-i-10
                                           layout(matrix(c(1:10), 1, 10, byrow = TRUE))
                                          }
                            }  # end i/10
                                     }   # end of 'while i < NFrames + 1'
                                     
# mtext(filnam, side = 1, line=-1, cex=0.9, col = 'red') #
write.table(Out, file=paste(filnam,irep,".txt", sep=''))
browser()
par(mar=c(4, 4, 4, 4) + 0.1)
layout(matrix(c(1:10), 1, 1, byrow = TRUE))

plot(Out[,1], Out[,3]-Out[,2], main = paste(filnam,"DIFF", "rep =",irep), ty='l')
lines(Out[,1], Out[,3]-Out[,2], typ='b', col='orange', lwd=2)
minOut3m2<- min(Out[,3]-Out[,2])
maxOut3m2<- max(Out[,3]-Out[,2])
dOut3m2<-maxOut3m2-minOut3m2
minV<-min(Out[2:Nframes,2]-Out[1:(Nframes-1),2])
maxV<-max(Out[2:Nframes,2]-Out[1:(Nframes-1),2])
dV<-maxV-minV
V<- minOut3m2 + dOut3m2*(Out[2:Nframes,2]-Out[1:(Nframes-1),2] -minV)/dV
lines(Out[2:Nframes,1]-0.5, V, col='blue', lwd=2)

ind<- seq(0,length(Out[,1]),10)
ind[1]<-1
OutX<-Out[ind,1]
OutY<- Out[ind,3]- Out[ind,2]
text(OutX, OutY, labels= OutX, pos=3)
mtext(" Blue = Velocity", side = 3, adj=0, line= -2, col = 'blue', cex=2)
mtext("Orange = Wall Thickness ", side = 3, adj=1, line= -2, col = 'orange', cex=2)

        irepeat<- as.real(ask("Enter 0 to end | 1 to add another analysis"))
if (irepeat>0) irep<-irep+1
         }  # end of while irepeat > 0