%nd2obslsf
% ND2OBS.m is a script that is an example of using frequency and amplitude output 
% from ADVOBS to clean a data record (or selected bursts) from station ND, Dep 2, 
% 2001 Grays Harbor Experiment, using OBSCLLSF.m.
%
% J.Lacy, USGS, November 2001
%
% Pick some bursts:
nmin = 6;
nmax = 200;
freq = 0.0047;
amp1 = 173;
amp2 = 180;
brange = [nmin:nmax];
fn = '/work3/jlacy/GH2001/GH2001d2/ND/ADV_B45B/ND2002.adr';
sta = input('enter station name in single quotes ')

a1 = zeros(length(brange),1);
a2 = zeros(length(brange),1);
phase1 = zeros(length(brange),1);
phase2 = zeros(length(brange),1);
error1 = zeros(length(brange),1);
error2 = zeros(length(brange),1);
mn1 = zeros(length(brange),1);
mn2 = zeros(length(brange),1);
std1 = zeros(length(brange),1);
std2 = zeros(length(brange),1);

%obscl = zeros(length(brange),2400,2);

%A = adr2mat(fn,[nmin,nmax]);
for i = nmin:nmax;
  ii = i - nmin+1;
A = adr2mat(fn,[i,i]);
[nr,nc] = size(A.extsensor{1});
del_t = 1/A.samprate(1);
 [obsgood1, amp1lsf, pha1,err1] = obscllsf(A.extsensor{1}(1:nr,1), del_t, freq, amp1);
%  axis ([0 n 200 1600])
  title (['OBS 1, Burst ',num2str(i)])
%  pause;
 [obsgood2, amp2lsf, pha2,err2] =obscllsf(A.extsensor{1}(1:nr,2), del_t, freq, amp2);
%   axis ([0 n 200 1600])
    title (['OBS 2, Burst ',num2str(i)])
%  pause;
   mn1(ii) = mean(A.extsensor{1}(1:nr,1));
   mn2(ii) = mean(A.extsensor{1}(1:nr,2));
   std1(ii) = std(A.extsensor{1}(1:nr,1));
   std2(ii) =  std(A.extsensor{1}(1:nr,2));
   if pha1<0;
     pha1 = pha1+2*pi;
   end
   if pha2<0;
     pha2 = pha2+2*pi;
   end
   phase1(ii) = pha1;
   phase2(ii) = pha2;
   a1(ii) = amp1lsf;
   a2(ii) = amp2lsf;
   error1(ii) = err1;
   error2(ii) = err2;
   obscl(ii,:,1) = obsgood1;
   obscl(ii,:,2)= obsgood2;
end

figure (4)
subplot (3,1,1)
bnum = [nmin:nmax];
plot (bnum,a1,'bo',bnum,a2,'g+')
title ('Amplitudes')
ylabel ('Counts')
subplot (3,1,2)
plot (bnum, error1,'bo',bnum,error2,'g+')
title ('Error')
%ylabel ('s^-1')
subplot (3,1,3)
plot (bnum, phase1, 'bo',bnum, phase2,'g+')
title('Phase')
ylabel ('Radians')
xlabel ('Burst number')
orient tall

figure (5)
subplot (2,2,1)
plot (mn1, a1,'o',mn2,a2,'+')
xlabel ('Burst mean')
ylabel ('Amplitude')
subplot (2,2,2)
plot (std1, a1, 'o',std2, a2, '+')
xlabel ('Burst std')
ylabel ('Amplitude')

ione = find(mn1<2000);
itwo = find(mn2<2000);

eval([sta,'mnamp1 = mean(a1(ione))'])
eval([sta,'mnamp2 = mean(a1(ione))'])

eval([sta,'medamp1 = median(a1(ione))'])
eval([sta,'medamp2 = median(a2(ione))'])


[h,sig] = ttest2(a1(ione), a2(itwo), 0.05, 0)

eval([sta,'amp1 = a1;']);
eval([sta,'amp2 = a2;']);

eval([sta, 'pha1 = phase1;']);
eval([sta, 'pha2 = phase2;']);

eval([sta, 'err1 = error1;']);
eval([sta, 'err2 = error2;']);

eval(['save ', sta, 'obslsf2 ', sta, '* mn1 mn2 ione itwo bnum freq'])
