# rm(list=ls())
# load("Test_18.Dta")
print(Sys.time())
Ta<-proc.time()

Vo<-as.matrix(Vo)
nx<-length(Vo[1,])
ns<-length(Vo[,1])

cat('First, sort out cosmic rays by adjacent spectra differences.') 
	mset<-200
    iA<-c(1,2,4,5)
for (i in 2:(ns-1)) {
	for (j in 3:(nx-2)){Vt<-Vo[i,j]
		                if (Vo[i,j]>mset) {if (Vo[i+1,j]<(Vt/5)) {if (Vo[i-1,j]<(Vt/5)){Mo<-median(Vo[i,c(j-2,j-1,j+1,j+2)])
					                                                                        Vo[i,j]<-Mo
					                                                                        }
	         			                                 }
	         			       }  
	              } # end for j
	  	      cat('.')
	          } # end for i
              cat('\n')

# Now look for adjusting within spectra.
cat('Now look for adjusting cosmic rays within spectra. ') 
# Ss<-Xm<-3:(ns-2)
Vn<-Vo
V1<-V2<-V3<-V4<-V5<-V6<-Ss<-Xs<-as.numeric(Xo[3:(nx-2)])
for(i in 1:ns){
	  	      # plot(Xo,Vo[i,],typ='l',main=paste('spectrum',i, 'min(So) =', round(min(Ss),2), '|  max(So) =', round(max(Ss),2)))
	  	      # browser()
	  for(j in 3:(nx-2)){
	  	                A<-Vo[i,(j-2):(j+2)]
	  	                So<-var(A[c(1,2,4,5)])^0.5
	  	                Mo<-mean(A[c(1,2,4,5)])
	  	                if(A[3]>Mo+5*So) {dit<-1; Vn[i,j]<-Mo } else {dit<-0}
	  	                Ss[j-2]<-So
	  	                }
	  	      cat('.')
	          }
              cat('\n')
# 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<-2

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

for(i in 1:ns){
	  for(j in 3:(nx-2)){
	  	    A<-Vo[i,(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[i,])
	  	                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[i,(j-2):(j+2)]
	  	      	                                    Vn[i,j]<- mean(A[c(1,2,4,5)])
	  	      	                                    }  else
	  	      	                                    
	  	      	                if (V5p[j]>  Vcrit) {
	  	      	                	if (V1p[j] > V2p[j]) {
	  	      	                                    A<-Vo[i,(j-2):(j+2)]
	  	      	                                    Vn[i,j]<- mean(A[c(1,2,4,5)])
	  	      	                                          }
	  	      	                                      }
	  	      	                                     }                 
	  	      	                }
	  	      cat('.')
	          }
              cat('\n')

print(proc.time()-Ta)