% cross_corr_curr.m
%
% Computes and plots cross-correlation function for 2 time series
%  -- currents in this case
% 1) interpolates the 2 series so that they are on the same time base
% 2) uses xcorr to calculate the function
% 3) plots the result
%
% SRV, created 8/21/01, updated 05/29/02, Flow Mow

% assumes wasp and cm50 structured arrays are loaded...
% load /shuksan/scottv/snap/flowmow/currents/mooring/no_filter.hourly/all.mat
% load /shuksan/scottv/snap/flowmow/currents/wasp/no_filter.hourly/wasp2161.mat

% define variables and project v and u into along/across axis components (20oN)
rot=20;     % approx overall strike of valley
%rot=17.35; % the trend of the MEF "box" (aligned with W valley fissures at MEF)
phi=atan2(mab50.v,mab50.u);
mag=sqrt(mab50.u.^2+mab50.v.^2);
u1=cos(phi+rot*pi/180).*mag;
v1=sin(phi+rot*pi/180).*mag;
t1=mab50.mjd;

phi=atan2(wasp2161.v,wasp2161.u);
mag=sqrt(wasp2161.u.^2+wasp2161.v.^2);
u2=cos(phi+rot*pi/180).*mag;
v2=sin(phi+rot*pi/180).*mag;
t2=wasp2161.mjd;

% establish temporal overlap and appropriate time base
maxt=min([max(t1) max(t2)])
mint=max([min(t1) min(t2)])
tres=1/24  % 1 hour resolution of data 
tbase=[mint:tres:maxt]';

v1int=interp1q(t1,v1,tbase);
v2int=interp1q(t2,v2,tbase);

% detrend
p1=polyfit(tbase,v1int,1);
p2=polyfit(tbase,v2int,1);

% check with a quick plot
 figure(20)
 hold; grid;
 %plot(t1,v1,'c')
 plot(tbase,v1int)
 %plot(t2,v2,'y')
 plot(tbase,v2int,'r')
 title('cm50 and WASP v on a 1hr t-base with detrend lines')
 plot(tbase,polyval(p1,tbase),'k')
 plot(tbase,polyval(p2,tbase),'k-.')
 %legend('t1int','t2int','t1fit','t2fit')
 %legend('S data','N data','S trend','N trend')
 axis([mint maxt -0.1 0.15])
 text(mint+(maxt-mint)/5, -0.085, 'S mean = +1.5cm/s')
 text(mint+(maxt-mint)/5, -0.065, 'N mean = +4.7cm/s')
 xlabel('time, Modified Julian Day')
 ylabel('Along-axis velocity, m/s')
 relabelx;
 print -depsc2 cm50_wasp_v_alongaxis.all.eps2
 %print -djpeg cm50_wasp_v_alongaxis.all.jpg

% Actual detrending calculation and mean subtraction
v1int=v1int-polyval(p1,tbase);
v2int=v2int-polyval(p2,tbase);
% the means after detrending are basically zero (e-11-13), so skip subtraction 
%v1int=v1int-mean(v1int);
%v2int=v2int-mean(v2int);

% Print out a file of the t base and 2 v series
fid = fopen('nsdat','w');
fprintf(fid,'%10.5f	%5.4f	%5.4f\n',[tbase'; v1int'; v2int']);
fclose(fid);

acor1=xcorr(v1int,v1int);
acor2=xcorr(v2int,v2int);
vcorr=xcorr(v1int,v2int);

figure(21);
clf; hold; grid;
plot(vcorr)
%plot(acor1,'r');
%plot(acor2,'g');
%legend('both','series 1', 'series 2')
title('Cross correlation function for S&N along-axis v: xcorr(v1,v2)')
print -depsc2 crosscorr_v_SvsN.eps2

coh=(vcorr.^2)./(acor1.*acor2);
figure(31);
clf; hold; grid;
plot(coh)
%plot(acor1,'r');
%plot(acor2,'g');
%legend('both','series 1', 'series 2')
title('Coherence function for S&N along-axis v, order (v1,v2)')
print -depsc2 coh_v_NvsS.eps2


if 0,
acov1=xcov(v1int,v1int);
acov2=xcov(v2int,v2int);
vcov=xcov(v1int,v2int);

figure(22);
clf; hold; grid;
plot(vcov)
plot(acov1,'r');
plot(acov2,'g');
legend('both','series 1', 'series 2')
title('Auto/covariance functions')

%find the power spectra of the autocorrelated data:
psvcorr=fft(vcorr);
npsvcorr=psvcorr .* conj(psvcorr)/(length(vcorr)*length(vcorr));
nyquist=1/(2*tres);  % days
fv=nyquist*(0:length(vcorr)/2)/(length(vcorr)/2);

psacor1=fft(acor1); 
npsacor1=psacor1 .* conj(psacor1)/(length(acor1)*length(acor1)); 
nyquist=1/(2*tres);  % days 
fa1=nyquist*(0:length(acor1)/2)/(length(acor1)/2);

psacor2=fft(acor2); 
npsacor2=psacor2 .* conj(psacor2)/(length(acor2)*length(acor2)); 
nyquist=1/(2*tres);  % days 
fa2=nyquist*(0:length(acor2)/2)/(length(acor2)/2);

figure(23);
clf; hold on; grid;
plot(fv*24,npsvcorr(1:((length(vcorr)-1)/2)+1,:))
plot(fa1*24,npsacor1(1:((length(acor1)-1)/2)+1,:),'r')
plot(fa2*24,npsacor2(1:((length(acor2)-1)/2)+1,:),'g')
axis([0 25 0 0.001])
xlabel('hours')
ylabel('~variance^2')
%legend('both','series 1', 'series 2')
legend('cross','cm50 auto', 'wasp auto')
title('Power spectra of detrended auto/cross correlation functions')
 print -dpsc power_corr_cm50_wasp.ps
 print -djpeg power_corr_cm50_wasp.jpg

end;


