# 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_2_dGauss.R which helps fit 2 Gaussian peaks to the data.\n")
cat("locations.  Abort if you have not used Wall_all on the log file yet!\n")

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

ask<- function(text) {
cat(text,": ", sep="")
indata<-readLines(con=stdin(),n=1,ok=TRUE)
                     }

                                       
 # End of Function Deffinitions #----------------------------------------
         
filnam<-"junk"

print(list.files("./","log"))
filnam<- print(ask("Enter log file name"))

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

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"))
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
Nframes<- length(Out[,1])
browser()
par(mar=c(4, 4, 4, 4) + 0.1)
layout(matrix(c(1:10), 2, 1, byrow = TRUE))
Wall<- Out[,3]-Out[,2]
Wx<-Out[,1]
plot(Out[,1], Wall, main = paste(filnam,"DIFF", "rep =",irep), typ='n')
lines(Out[,1], Wall, col='orange', lwd=2, typ='p')
minOut3m2<- min(Wall)
maxOut3m2<- max(Wall)
dOut3m2<-maxOut3m2-minOut3m2
Veloc<- Out[2:Nframes,2]-Out[1:(Nframes-1),2]
Vx<- Out[2:Nframes,1]-0.5
minV<-min(Veloc)
maxV<-max(Veloc)
dV<-maxV-minV
V<- minOut3m2 + dOut3m2*(Veloc -minV)/dV
lines(Vx, V, col='blue', lwd=2, typ='p')

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)
lines(Wl<-lowess(Wx,Wall, f=0.02), col='orange',typ='l', lwd=3)
lines(Wlb<-lowess(Wx,Wall, f=0.2), col='orange',typ='l', lwd=2)
lines(Vl<-lowess(Vx,V, f=0.02), col='blue',typ='l', lwd=3)
lines(Vlb<-lowess(Vx,V, f=0.2), col='blue',typ='l', lwd=2)

# Corrected

minX<- min(Out[,1]); maxX<- max(Out[,1]); maxY<- max(abs(Wall-Wlb$y))
plot(c(minX,maxX), c(-maxY,maxY), main = paste(filnam,"DIFF", "rep =",irep), typ='n')
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)
lines(Wx<-Out[,1], Wy<-Wall-Wlb$y, col='orange', lwd=2, typ='p')
lines(Vx, Vy<-V-Vlb$y, col='blue', lwd=2, typ='p')
Wlx<-Wl$x; Wly<-Wl$y-Wlb$y
lines(Wl$x, Wl$y-Wlb$y, col='orange',typ='l', lwd=3)
Vlx<-Vl$x; Vly<-Vl$y-Vlb$y
lines(Vl$x, Vl$y-Vlb$y, col='blue',typ='l', lwd=3)

browser()
AutoConv<- convolve(rev(Vy), rev(Vy),  type="o")  # Gives CORRECT Orientation
AutoIonv<- convolve(rev(Vy), rev(-Vy),  type="o") # Gives CORRECT Orientation
minX<- min(1:length(AutoConv)); maxX<- max(1:length(AutoConv)); maxY<- max(abs(AutoConv))
plot(c(minX,maxX), c(-maxY,maxY), main = paste(filnam,"Autoconvolution", "rep =",irep), typ='n')
Xmid<- (minX+maxX)/2
lines(c(Xmid,Xmid),c(-maxY,maxY))
lines(AutoConv, col="blue", typ='l')
lines(AutoIonv, col="orange")
#spline Wall data onto Velocity time base
     secsper<-1
     junk<-spline(Wx, Wy, n = secsper*length(Vx), method ="fmm", xmin=min(Vx), xmax=max(Vx))     
     Wsy<- junk$y
     Wsx<- junk$x
Conv<- convolve(rev(Vy), rev(Wsy),  type="o")  # Gives CORRECT Orientation
Ionv<- convolve(rev(Vy), rev(-Wsy),  type="o") # Gives CORRECT Orientation
minX<- min(1:length(Conv)); maxX<- max(1:length(Conv)); maxY<- max(abs(Conv))
plot(c(minX,maxX), c(-maxY,maxY), main = paste(filnam,"Convolution", "rep =",irep), typ='n')
lines(c(Xmid,Xmid),c(-maxY,maxY))
lines(Conv, col="blue", typ='l')
lines(Ionv, col="orange")

browser()
AutoConv<- convolve(rev(Vly), rev(Vly),  type="o")  # Gives CORRECT Orientation
AutoIonv<- convolve(rev(Vly), rev(-Vly),  type="o") # Gives CORRECT Orientation
minX<- min(1:length(AutoConv)); maxX<- max(1:length(AutoConv)); maxY<- max(abs(AutoConv))
plot(c(minX,maxX), c(-maxY,maxY), main = paste(filnam,"Autoconvolution of lowess fit", "rep =",irep), typ='n')
Xmid<- (minX+maxX)/2
lines(c(Xmid,Xmid),c(-maxY,maxY))
lines(AutoConv, col="blue", typ='l')
lines(AutoIonv, col="orange")
#spline Wall data onto Velocity time base
     secsper<-1
     junk<-spline(Wlx, Wly, n = secsper*length(Vx), method ="fmm", xmin=min(Vlx), xmax=max(Vlx))     
     Wlsy<- junk$y
     Wlsx<- junk$x
Conv<- convolve(rev(Vly), rev(Wlsy),  type="o")  # Gives CORRECT Orientation
Ionv<- convolve(rev(Vly), rev(-Wlsy),  type="o") # Gives CORRECT Orientation
minX<- min(1:length(Conv)); maxX<- max(1:length(Conv)); maxY<- max(abs(Conv))
plot(c(minX,maxX), c(-maxY,maxY), main = paste(filnam,"Convolution of lowess fit", "rep =",irep), typ='n')
lines(c(Xmid,Xmid),c(-maxY,maxY))
lines(Conv, col="blue", typ='l')
lines(Ionv, col="orange")

browser()
Ind<- 7:(length(Vly)-6)
layout(matrix(c(1:10), 2, 5, byrow = TRUE)) 
plot(Vly[Ind], Wly[Ind+2])
plot(Vly[Ind], Wly[Ind+3])
plot(Vly[Ind], Wly[Ind+4])
plot(Vly[Ind], Wly[Ind+5])
plot(Vly[Ind], Wly[Ind+6])

plot(Vly[Ind], Wly[Ind-2])
plot(Vly[Ind], Wly[Ind-3])
plot(Vly[Ind], Wly[Ind-4])
plot(Vly[Ind], Wly[Ind-5])
plot(Vly[Ind], Wly[Ind-6])

browser()


Ind<- 8:(length(Vly)-7)
Simdat<-matrix(1,length(Ind),12)
Simdat[,1]<-Vy[Ind]
layout(matrix(c(1:10), 2, 5, byrow = TRUE)) 
for (i in 3:7) {
plot(Vy[Ind], Wy[Ind+i], main=filnam)
Cvw<-cor(Vy[Ind], Wy[Ind+i])
mtext(paste("rsq=",round(Cvw^2,4)," "), side = 3, adj=1, line= -2, col = 'blue', cex=1 )
mtext(paste("Wall shift=",i," "), side = 3, adj=1, line= -4, col = 'blue', cex=1 )
Simdat[,i-1]<-Wy[Ind+i]
               }

for (i in -3:-7) {
plot(Vy[Ind], Wy[Ind+i], main=filnam)
Cvw<-cor(Vy[Ind], Wy[Ind+i])
mtext(paste(" rsq=",round(Cvw^2,4)), side = 3, adj=0, line= -2, col = 'blue', cex=1 )
mtext(paste("Wall shift=",i," "), side = 3, adj=0, line= -4, col = 'blue', cex=1 )
Simdat[,4+abs(i)]<-Wy[Ind+i]
                }
write.table(Simdat, file=paste(filnam,".out.csv",sep=''), sep=",",row.names = FALSE, quote = FALSE)