%pc_obs.m
%
%pc_obs finds the frequency and amplitude of the obs-interference
%signal from pcadp .mat files, by calling findfreq.m
%
% J.Lacy, USGS, November 2001
%
fn = input ('enter name of mat file of pcadp data (no ext) in single quotes ') 
%eval (['load /work3/jlacy/GH2001mat/Dep2/MIA/pcadp/',fn])
eval (['load /work3/jlacy/GH2001mat/Dep1/MIA/pc_h40/',fn])
%pick bursts
nmin= 3;
nmax = 200;
fr1a = zeros(nmax,1);
fr2a = zeros(nmax,1);
amponea = zeros(nmax,1);
amptwoa = zeros(nmax,1);
phaonea = zeros(nmax,1);
phatwoa = zeros(nmax,1);
mn1 = zeros(nmax,1);
mn2 = zeros(nmax,1);
std1 = zeros(nmax,1);
std2 = zeros(nmax,1);
for i = nmin:nmax;
   ii = i - nmin+1;
  first = pcadp.bur_ind(i);
  last = pcadp.bur_ind(i)+pcadp.ppb-1;
  del_t = pcadp.avg_int;
  [obsgood1, frmax1, amp1, pha1] = findfreq(pcadp.ad1(first:last), del_t);
%  axis ([0 1200 200 1600])
  title (['OBS 1, Burst ',num2str(i)])
%  pause;
   [obsgood2, frmax2, amp2, pha2] = findfreq(pcadp.ad2(first:last), del_t);
%   axis ([0 1200 200 1600])
    title (['OBS 2, Burst ',num2str(i)])
%   pause;
   fr1a(i) = frmax1/del_t;
   amponea(i) = amp1;
   phaonea(i) = pha1;
   fr2a(i) = frmax2/del_t;
   amptwoa(i) = amp2;
   phatwoa(i) = pha2;
    mn1(i) = mean(pcadp.ad1(first:last));
   mn2(i) = mean(pcadp.ad2(first:last));
   std1(i) = std(pcadp.ad1(first:last));
   std2(i) = std(pcadp.ad2(first:last));
   temp(i) = mean(pcadp.t(first:last));
end
figure 
subplot (3,1,1)
bnum = [1:1:nmax];
plot (bnum,amponea,'bo',bnum,amptwoa,'g+')
title ('Amplitudes')
ylabel ('Counts')
subplot (3,1,2)
plot (bnum, fr1a,'bo',bnum,fr2a,'g+')
title ('Frequencies')
ylabel ('s^-1')
subplot (3,1,3)
plot (bnum, phaonea, 'bo',bnum, phatwoa,'g+')
title('Phase')
ylabel ('Radians')
xlabel ('Burst number')

figure 
subplot (2,2,1)
plot (mn1(nmin:nmax), amponea(nmin:nmax),'o',mn2(nmin:nmax), ...
      amptwoa(nmin:nmax),'+')
xlabel ('Burst mean')
ylabel ('Amplitude')
subplot (2,2,2)
plot (std1(nmin:nmax), amponea(nmin:nmax), 'o',std2(nmin:nmax), amptwoa(nmin:nmax), '+')
xlabel ('Burst std')
ylabel ('Amplitude')
subplot (2,2,3)
plot (mn1(nmin:nmax), fr1a(nmin:nmax),'o',mn2(nmin:nmax),fr2a(nmin:nmax),'+')
ylabel ('Frequency')
xlabel ('Burst mean')
subplot (2,2,4)
plot (temp(nmin:nmax), amponea(nmin:nmax), 'o',temp(nmin:nmax),amptwoa(nmin:nmax),'+')
xlabel('Temperature')
ylabel('Amplitude')

ione = find(mn1(nmin:nmax)<1000)+nmin-1;
itwo = find(mn2(nmin:nmax)<1000)+nmin-1;

mnamp1 = mean(amponea(ione))
mnamp2 = mean(amptwoa(ione))

medamp1 = median(amponea(ione))
medamp2 = median(amptwoa(ione))

mnfr1a = mean(fr1a(ione))
mnfr2a = mean(fr2a(itwo))

medfr1a = median(fr1a(ione))
medfr2a = median(fr2a(itwo))

[h,sig] = ttest2(amponea(ione), amptwoa(itwo), 0.05, 0)