rm(list=ls())
pref<-"Jan 3 2010 preliminary HAP Raman spectra:"

RamSpecEval<- function(filename) {
	XX<-read.table(filename)
attach(XX)
quartz(width=10, height=5)
plot(V1,V2, typ='l', , main=paste(pref,filename, 'raw'), col='red')
# browser()
cat("Zero order is at:",Zo<-V1[which.max(V2)],"\n")
Xo<-V1[V1>149.999]
Yo<-V2[V1>149.999]
detach(XX)
Vn<-Vo<- Yo - lowess(Xo, Yo,f=0.2)$y
nx<-length(Xo)
cat("length(Xo)=", nx, ".\n")
V1<-V2<-V3<-V4<-V5<-V6<-Ss<-Xs<-as.numeric(Xo[3:(nx-2)])
	  for(j in 3:(nx-2)){
	  	                A<-Vo[(j-2):(j+2)]
	  	                So<-var(A[c(1,2,4,5)])^0.5
	  	                Mo<-median(A[c(1,2,4,5)])
	  	                if(A[3]>Mo+5*So) {dit<-1; Vn[j]<-Mo } else {dit<-0}
	  	                Ss[j-2]<-So
	  	                }

# Orthogonal multipliers 
XX<- c(-2, -1, 1, 2,  -1,  1, 1, -1, 1,  -2, 2,  -1)
XX<-matrix(XX, 3, 4, byrow= TRUE)
XX5<- matrix(0,3,5)
XX5[,c(1,2,4,5)]<-XX
XX5[2,]<- c(-2,1,4,1,-2)
SVD5<-svd(XX5)
SVD4<-svd(XX)
SVD4$v<-SVD4$v[,c(1,3,2)]
SVD4$d<-SVD4$d[c(1,3,2)]
SVD4$v[,2]<- -1*SVD4$v[,2]
SVD5$v<-SVD5$v[,c(3,1,2)]
SVD5$d<-SVD5$d[c(3,1,2)]
Perc<-0.2
Vcrit<-1.5

cat('Now look for adjusting linear vs quadratic within spectra. ') 

	  for(j in 3:(nx-2)){
	  	    A<-Vo[(j-2):(j+2)]
	  	    S4<-t(SVD4$v)%*%A[c(1,2,4,5)]/SVD4$d
	  	    S5<-t(SVD5$v)%*%A/SVD5$d
	  	                V1[j]<-S4[1]
	  	                V2[j]<-S4[2]
	  	                V3[j]<-S4[3]
	  	                V4[j]<-S5[1]
	  	                V5[j]<-S5[2]
	  	                V6[j]<-S5[3]
	  	                }
	  	               
	  	                MaxV<-Perc*max(Vn)
	  	                Vsum<-abs(V1)+abs(V2)+abs(V3)
	  	                V1p<-MaxV*V1/Vsum
	  	                V2p<-MaxV*V2/Vsum
	  	                V3p<-MaxV*V3/Vsum
	  	                V5p<-MaxV*V3/Vsum

	  	      for(j in 3:(nx-2)){ if(!is.nan(V5p[j])){
	  	      	                if (V5p[j]< (-Vcrit)) {
	  	      	                                    A<-Vo[(j-2):(j+2)]
	  	      	                                    Vn[j]<- median(A[c(1,2,4,5)])
	  	      	                                    }  else
	  	      	                                    
	  	      	                if (V5p[j]>  Vcrit) {
	  	      	                	if (V1p[j] > V2p[j]) {
	  	      	                                    A<-Vo[(j-2):(j+2)]
	  	      	                                    Vn[j]<- median(A[c(1,2,4,5)])
	  	      	                                          }
	  	      	                                      }
	  	      	                                     }                 
	  	      	                }
quartz(width=10, height=5)
plot(Xo,Vn, typ='l', main=paste(pref,filename, 'fixed'))
mtext(paste("Zero order is at: ", Zo))
}

# End of subroutine deffinition
#-------------------------------------------

fnam="STD_triCaPO4__14_54_105.txt"
RamSpecEval(fnam)
browser()

fnam="STD_triCaPO4__15_05_107.txt"
RamSpecEval(fnam)
browser()

fnam="STD_triCaPO4__15_08_108.txt"
RamSpecEval(fnam)
browser()

fnam="STD_triCaPO4__15_16_110.txt"
RamSpecEval(fnam)
browser()

fnam="STD_triCaPO4__15_33_112.txt"
RamSpecEval(fnam)
browser()

fnam="STD_triCaPO4__15_33_113.txt"
RamSpecEval(fnam)
browser()

fnam="STD_triCaPO4__15_34_114.txt"
RamSpecEval(fnam)
browser()

fnam="STD_triCaPO4__15_38_115.txt"
RamSpecEval(fnam)
browser()

fnam="STD_triCaPO4__15_38_116.txt"
RamSpecEval(fnam)
browser()

fnam="STD_triCaPO4__15_39_117.txt"
RamSpecEval(fnam)