rm(list=ls())
load("Test_18.Dta")
Vo<-as.matrix(Vo)
nx<-length(Vo[1,])
ns<-length(Vo[,1])


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

# Now look for adjusting 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

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
	  	      #quartz(width=12, height=7.5)
	  	      #plot(Xs,Vsn[i,],typ='l',main=paste('wavelength',i, 'linear'))
	  	      #lines(Xs,V1[3:(nx-2)], col='red', lwd=3)
	  	      #lines(Xs,V4[3:(nx-2)], col='green')
	  	      #browser()
	  	      #quartz(width=12, height=7.5)
	  	      #plot(Xs,Vsn[i,],typ='l',main=paste('wavelength',i, 'quadratic'))
	  	      #lines(Xs,V2[3:(nx-2)], col='red', lwd=3)
	  	      #lines(Xs,V5[3:(nx-2)], col='green')
	  	      #browser()
	  	      #quartz(width=12, height=7.5)
	  	      #plot(Xs,Vsn[i,],typ='l',main=paste('wavelength',i, 'cubic'))
	  	      #lines(Xs,V3[3:(nx-2)], col='red', lwd=3)
	  	      #lines(Xs,V6[3:(nx-2)], col='blue')
	  	      #browser()
	  	      #quartz(width=12, height=7.5)
	  	      #plot(Xs,Vsn[i,],typ='l',main=paste('wavelength',i, 'cubic'))
	  	      #lines(Xs,V1p[3:(nx-2)], col='red', lwd=3)
	  	      #lines(Xs,V2p[3:(nx-2)], col='green')
	  	      #lines(Xs,V3p[3:(nx-2)], col='blue')
	  	      #lines(Xs,V5p[3:(nx-2)], col='purple')
	  	      #browser()
	  	      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)])
	  	      	                                          }
	  	      	                                      }
	  	      	                                     }                 
	  	      	                }
	  	      #plot(Xo,Vo[i,],typ='l',main=paste('spectrum',i), lwd=2)
	  	      #lines(Xo,Vn[i,],typ='l',col='red')
	  	      #browser()	          
	  	      cat('.')
	          }
              cat('\n')
#	  	           }  # end for i