rm(list=ls())
print(Sys.time())
print(Ta<-proc.time())
filnam<-"BL20_03_16_32_26_25x25" # was "M3a_canalAF-Ram_11_10_1"
nrows=25
ncols=25

## BEGIN OF FUNCTION DEFFINITION ##

RamSpecEval<- function(X0, Y0, FF, NC) {

Yi<- lowess(X0, Y0,f=FF)$y
Ylo<- Y0 - Yi
nx<-length(X0)
for(i in 1:NC) {
for(j in 1:nx) { Ylo[j]<- min(Y0[j],Yi[j])
	              }	             
	           Yi<- lowess(X0, Ylo,f=FF)$y
                }
Vy<-Vx<- Y0-Yi
cat('.')       # Added by Joe
Vy
}
## END OF FUNCTION DEFFINITION ##

XX<-read.table(paste(filnam,".txt",sep=''), nrows=nrows*ncols+1)        # 2296   {1684 works}
Lx<-length(as.matrix(XX)[1,])
Ly<-length(as.matrix(XX)[,3])
Xo<- as.numeric(XX[1,3:Lx])
XYloc<- XX[2:Ly,1:2]
XlocUnique<-unique(sort(XYloc[,2]))
YlocUnique<-unique(sort(XYloc[,1]))
Vo<-Yo<- XX[2:Ly,3:Lx]
Xind<- 170:3001
Xo<-Xo[Xind]
Yo<-Vo<-Vo[ ,Xind]
Lx<-Lx-2
Lx<-length(Xo)
Ly<-Ly-1

cat('Begin RamSpecEval Function\n') 
for (i in 1:Ly) {Vo[i,]<- RamSpecEval(Xo, Yo[i,], 0.1, 10 ) }     ##     Yo[i,] - lowess(Xo, Yo[i,],f=0.1)$y}
# This script is currently experimental
cat('Ended RamSpecEval Function\n') 

mset<-100
Yadj<- 100       #1300
for (i in 2:(Ly-1)){
for (j in 3:(Lx-2)){if (Vo[i,j]>mset) {if (Vo[i+1,j]<(mset/5)) {if(Vo[i-1,j]<(mset/5)){Vo[i,j]<-0}}}
} # end for j
	  	      cat('.')       # Added by Joe
	          }# end for i
              cat('\n')

minX<- min(Xo); maxX<- max(Xo)
minY<- min(Vo); maxY<- max(Vo)
save(filnam,Vo,Xo,XlocUnique,YlocUnique, file=paste(filnam,"_PE.Dta", sep=''))
quartz(width=10.5, height=8.5)
plot(c(minX, maxX+3*Ly), c(minY, maxY+Yadj), typ='n',main=paste(filnam,' spectra 1 to', Ly,'mset =', mset, 'Yadj =',Yadj))
for (i in 1:Ly){lines(Xo+3*(Ly-i), Vo[i,]+3*(Ly-i), col=rainbow(Ly)[i])
if(0==((i/ncols)-floor(i/ncols))){lines(Xo+3*(Ly-i),Vo[i,]+3*(Ly-i),col='black', lwd=2)
text(Xo[Lx]+3*(Ly-i), Vo[i,Lx]+3*(Ly-i), i,col='black', pos=4, cex=0.6)
}
	  	      cat('.')       # Added by Joe
	          }# end for i
              cat('\n')
# Computations of SVD
Vsvd<-svd(Vo)
VV<-Vsvd$v[,1:16]
VD<-Vsvd$d[1:16]
VDs<- sum(Vsvd$d)
VDp<-100*VD/VDs

VV<-Vsvd$v[,1:16]
cat('Top 16 SVDs percent:',100*Vsvd$d[1:16]/sum(Vsvd$d),'\n')
quartz()
Ev<-t(VV)%*%t(Vo)
PCV1<-matrix(Ev[1,], nrows,ncols, byrow=TRUE)
filled.contour(YlocUnique, XlocUnique, PCV1, main=paste(filnam, "PCV1"), color = topo.colors)
quartz()
PCV2<-matrix(Ev[2,], nrows,ncols, byrow=TRUE)
filled.contour(YlocUnique, XlocUnique, PCV2, main=paste(filnam, "PCV2"), color = topo.colors)
quartz()
PCV3<-matrix(Ev[3,], nrows,ncols, byrow=TRUE)
filled.contour(YlocUnique, XlocUnique, PCV3, main=paste(filnam, "PCV3"), color = topo.colors)
quartz()
PCV4<-matrix(Ev[4,], nrows,ncols, byrow=TRUE)
filled.contour(YlocUnique, XlocUnique, PCV4, main=paste(filnam, "PCV4"), color = topo.colors)
quartz(width=10.5, height=6)
plot(c(min(Xo), max(Xo)),c(-0.5,0.5), typ='n', main=paste(filnam,'Top 6 SVDs'))
lines(Xo,-1*VV[,1], typ='l', col='black')
lines(Xo,VV[,2], typ='l', col='red')
lines(Xo,VV[,3], typ='l', col='green')
lines(Xo,VV[,4], typ='l', col='blue')
lines(Xo,VV[,5], typ='l', col='purple')
lines(Xo,VV[,6], typ='l', col='brown')
Apaum<-c(430, 449, 581, 590, 609, 960,1047, 1070); yup<-0.4; di<-0.02
for(i in 1:8) { lines(c(Apaum[i],Apaum[i]),c(-0.45,yup),col='brown',lwd=2, lty=2)
text(Apaum[i], yup+di*(i%%3), Apaum[i], pos=3, cex=0.66)}
Ap<-1318
lines(c(Ap,Ap), c(-0.45,0.01), col='purple', lwd=2, lty=2)
text(Ap, -0.45, Ap, pos=1, cex=0.66)

for (i in 1:6) {text(5000,0.5-i*0.05,paste("SDV",i,'=',round(VDp[i],2)))}
print(proc.time()-Ta)