%
% set the window 
% load the data
%
%
Pstart = 2200;
Pend = Pstart + 1000;
% surface  Z,N,E
load HLZ.dat;
load HLN.dat;
load HLE.dat;
% surface - Z,N,E
load HL7.dat;
load HL8.dat;
load HL9.dat;
% 91 m borehole - Z,N,E 
load HL1.dat;
load HL2.dat;
load HL3.dat;
% 45 m borehole - Z,N,E 
load HL4.dat;
load HL5.dat;
load HL6.dat;
%
%
HL2 = HL2 - mean(HL2);
HL3 = HL3 - mean(HL3);
%
HL5 = HL5 - mean(HL5);
HL6 = HL6 - mean(HL6);
%
% filter with a 3rd order butterworth lowpass
% matlab normalizes to Nyquist
% load into new array 
%
Wn = [0.01 0.02];
[B,A] = butter(3,Wn);
%
tmp = HLN - mean(HLN);
RefNgram = filter(B,A,tmp);
tmp = HLE - mean(HLE);
RefEgram = filter(B,A,tmp);
clear tmp;
%
%
%
rad = 180.0/pi;
factor = 1;
for j = 1:2
	if (j == 1)
		Pgram(:,1) = filter(B,A,HL2);
		Pgram(:,2) = filter(B,A,HL3);
	end
	if (j == 2) 
		Pgram(:,1) = filter(B,A,HL5);
		Pgram(:,2) = filter(B,A,HL6);
	end
%
	for i = factor:factor:360
		angle = i/rad;
%
% set up rotation matrix
%
		rot(1,1) = cos(angle);
		rot(1,2) = -1*sin(angle);
		rot(2,1) = sin(angle);
		rot(2,2) = cos(angle);
%
		rnew = Pgram*rot;
		CHAN1 = rnew(:,1);
		CHAN2 = rnew(:,2);
%
% now calculate the covariance between the rotated seismogram
% and the original known seismogram
%
		c = cov(CHAN1(Pstart:Pend),RefNgram(Pstart:Pend));
%
		[evec,eval] = eig(c);
		c11(i) = c(1,1);
		c22(i) = c(1,1);
		lam1(i) = max(eval(1,1),eval(2,2));
		lam2(i) = min(eval(1,1),eval(2,2));
		flin(j,i) = (abs(lam2(i))/abs(lam1(i)));
		theta(i) = rad*atan((evec(2,1))/(evec(1,1)));
%
		error1N(j,i) = sum(abs(CHAN1(Pstart:Pend) - RefNgram(Pstart:Pend)));
		error1E(j,i) = sum(abs(CHAN2(Pstart:Pend) - RefEgram(Pstart:Pend)));
%
		error2N(j,i) = sum(abs(CHAN2(Pstart:Pend) - RefNgram(Pstart:Pend)));
		error2E(j,i) = sum(abs(CHAN1(Pstart:Pend) - RefEgram(Pstart:Pend)));
%
		e1(j,i) = error1N(j,i) + error1E(j,i);
		e2(j,i) = error2N(j,i) + error2E(j,i);
	end
end
%
 [y,i] = min(e1(1,:))
 [y,i] = min(e2(1,:))
%
 [y,i] = min(e1(2,:))
 [y,i] = min(e2(2,:))
%
plot(e1(2,:),'green');
hold
plot(e2(2,:),'red');