% make_currfigs.m
% 9/26/02 NOTE old.make_currfigs.m creates figures that may be good for thesis, 
%      especially the phase/coh plot, excluded from fm1 paper here...
%
% Computes and plots coherence and phase as functions 
% of frequency for 2 scalar time series  -- v component currents in this case
%
% 1) interpolates the 2 series so that they are on the same time base
% 2) uses psd to calculate detrend, window, and dft each series 
% 3) calculates cross and auto energy density estimates (G12,G11,G22)
% 4) calculates and plots coh and phase from G12, G11, and G22 
%
% SRV, created 05/31/02, Flow Mow
% SRV, 7/5/02: adapted from /flowmow/processing/currents/coh/cohphase.m 
% to create stats and figures for horizontal flux paper
%
% 1) compute means for vectors in location figure
% 2) generate cohphase figure (with PSD in upper subplots?)
% 3) t-series example (quivers N vs S, plus v-component comparison N/S/abe/tivey) 
%

local=1;

% load wasp and fm S mooring structured arrays...
load_fmcm; % m/s 
load_wasp; % m/s
load_cm95; % speed in m/s
load /snap/flowmow/ABE/abe50/abe50_currents.mat; % 10sec data in m/s 
load /snap/flowmow/ABE/abe51/abe51_currents.mat; % 10sec data in m/s
cm2175.u=cm2175.spd_mps.*cos(pi/180*(90-cm2175.dirT));
cm2175.v=cm2175.spd_mps.*sin(pi/180*(90-cm2175.dirT));
cm1900.v=cm1900.spd_mps.*sin(pi/180*(90-cm1900.dirT));
cm1900.u=cm1900.spd_mps.*cos(pi/180*(90-cm1900.dirT)); 

% calculate means and for plotting in gmt script: ~/thesis/fm1/figs/make_fmloc.pl
meancm2175u=mean(cm2175.u); meancm2175v=mean(cm2175.v);
meancm1900u=mean(cm1900.u); meancm1900v=mean(cm1900.v);
meanmab50u=mean(mab50.u); meanmab50v=mean(mab50.v);
meanmab250u=mean(mab250.u); meanmab250v=mean(mab250.v);
meanwasp2161u=mean(wasp2161.u); meanwasp2161v=mean(wasp2161.v);
meancm2175dir=90-(180*atan2(meancm2175v,meancm2175u)/pi) 
meancm2175mag=sqrt(meancm2175u.^2+meancm2175v.^2)
meancm1900dir=90-(180*atan2(meancm1900v,meancm1900u)/pi)
meancm1900mag=sqrt(meancm1900u.^2+meancm1900v.^2)
meanmab50dir=90-(180*atan2(meanmab50v,meanmab50u)/pi)
meanmab50mag=sqrt(meanmab50u.^2+meanmab50v.^2)
meanmab250dir=90-(180*atan2(meanmab250v,meanmab250u)/pi)
meanmab250mag=sqrt(meanmab250u.^2+meanmab250v.^2)
meanwasp2161dir=90-(180*atan2(meanwasp2161v,meanwasp2161u)/pi) 
meanwasp2161mag=sqrt(meanwasp2161u.^2+meanwasp2161v.^2)


% project v and u into along/across axis components 
% (17.35N to be orthogonal to the heat flux survey N and S surfaces)
%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;

phi=atan2(curr50.vsmooth,curr50.usmooth);
mag=sqrt(curr50.usmooth.^2+curr50.vsmooth.^2);
u50=cos(phi+rot*pi/180).*mag;
v50=sin(phi+rot*pi/180).*mag;
t50=curr50.mjd;

if 1, 
phi=atan2(curr51.vsmooth,curr51.usmooth);
mag=sqrt(curr51.usmooth.^2+curr51.vsmooth.^2);
u51=cos(phi+rot*pi/180).*mag;
v51=sin(phi+rot*pi/180).*mag;
t51=curr51.mjd;
end;


% 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);
u1int=interp1q(t1,u1,tbase);
u2int=interp1q(t2,u2,tbase);

if 0,
tbase50=tbase(find(tbase>min(t50) & tbase<max(t50)) );
v50int=interp1q(t50,v50,tbase50);
u50int=interp1q(t50,u50,tbase50);
%tbase51=find(tbase>min(t51) & tbase<max(t51));
%v51int=interp1q(t51,v51,tbase51);
end;

if 0, 
% 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);
end;

% Create time-series figure for fm1 

% load or rename predicted, mav40, N (wasp2000) and S (fm2000/ois) 
un=u2*100; vn=v2*100; us=u1*100; vs=v1*100; % NOTE m/s to cm/s conversion here!
u50cm=u50*100; v50cm=v50*100;

load /snap/flowmow/currents/predicted/pred.mat;  % v (true N) in cm/s
load /snap/flowmow/currents/tivey/mav40proc.mat;  % v (true N) in cm/s

% transform meg's mavs to along/across axis components
phi=atan2(mav40.v,mav40.u);
mag=sqrt(mav40.u.^2+mav40.v.^2);
mav40.ua=cos(phi+rot*pi/180).*mag;
mav40.va=sin(phi+rot*pi/180).*mag;

% transform abe currents here

% get the abe50/51 interval for comparative time-series plotting efficiency
time_table;
mjdi=tm.strd50+tm.mor150(1)/86400;
mjdf=tm.strd51+tm.top651(2)/86400;
ni=find(wasp2161.mjd>=mjdi & wasp2161.mjd<=mjdf);
si=find(mab50.mjd>=mjdi & mab50.mjd<=mjdf);
%mav40i=find(mav40.mjd>=mjdi & mav40.mjd<=mjdf);
predi=find(pred.mjd>=mjdi & pred.mjd<=mjdf);
mi=find(mav40.mjd>=mjdi & mav40.mjd<=mjdf);

%smooth meg's mavs velocities with a low-pass filter and interpolate to 1 hr base
[b,a] = butter(10,.05); % 10th order, cut off frequency = 0.05*(samplerate/2)
% with sampling rate =2/min=0.033/s, this gives a cutoff freq of ~1300s, or period~22min
mav40.uas=filtfilt(b,a,mav40.ua);
mav40.vas=filtfilt(b,a,mav40.va);
% maybe better to create hourly averages centered on the tbase than to interpolate...
if 0,
% maybe better to create hourly averages centered on the tbase than to interpolate...
if isempty(mav40.vahr),
mav40.vahr=[]; mav40.uahr=[]; mav40.thr=[];
for tint=mint:tres:maxt,  % this is a very slow loop...
 if (mav40.mjd(length(mav40.mjd))>tint-tres/2),
  if (mav40.mjd(1)<(tint+tres/2)),
   block=find(mav40.mjd>tint-tres/2 & mav40.mjd<=tint+tres/2);
   mav40.vahr=[mav40.vahr; mean(mav40.vas(block))];
   mav40.uahr=[mav40.uahr; mean(mav40.uas(block))];
   mav40.thr=[mav40.thr; tint];
  end;
 end;
end;
end;
%% ARGH. Waste of time. 
%% The above averaging scheme is equivalent to the following interpolation! (and much slower)
mav40.uasi=interp1q(mav40.mjd,mav40.uas,tbase);
mav40.vasi=interp1q(mav40.mjd,mav40.vas,tbase);
%mav40.uasi=interp1q(mav40.mjd,mav40.uas,wasp2161.mjd(ni));
%mav40.vasi=interp1q(mav40.mjd,mav40.vas,wasp2161.mjd(ni));
end;

%% Plot fm1 version of comparison feather plot
figure(41)
clf;

subplot(2,1,1);
%title('Comparison of near-bottom velocities: N mooring, S mooring, and at MEF');
quiver(24*(tbase-mjdi),zeros(size(tbase)),u2int*100,v2int*100,0,'b.')
grid
axis equal
axis([0    65   -12    12]);
text(5,-7,'North meter: year 2000, depth 2161 m, elevation 15 mab')
set(gca,'XTickLabel',[]);
ylabel('Velocity [cm/s]');

subplot(2,1,2);
quiver(24*(tbase-mjdi),zeros(size(tbase)),u1int*100,v1int*100,0,'r.')
grid
axis equal
axis([0    65   -12    12]);
%set(gca,'XTickLabel',[]);
text(5,-7,'South meter: year 2000, depth 2142 m, elevation 50 mab')
ylabel('Velocity [cm/s]');

xlabel('Time [hr]')
%xlabel('Elapsed time [hr]')

print -depsc2 fm1_quiver.eps2
print -djpeg  fm1_quiver.jpg


% Make FM1 plot of zero-mean oscillations
if 1,
figure(42)
clf

subplot(2,1,1);
hold on;
grid
plot((pred.mjd(predi)-mjdi)*24,pred.v(predi),'k');
plot((wasp2161.mjd(ni)-mjdi)*24,vn(ni)-mean(vn(ni)),'b--');
plot((mab50.mjd(si)-mjdi)*24,vs(si)-mean(vs(si)),'g');
axis([0 75 -12 12])
legend('Predicted','North','South')
%title('Zero-mean along-axis V comparison: N/S moorings and MEF');
ylabel('v [cm/s]')

subplot(2,1,2);
hold on;
grid
plot((curr50.mjd-mjdi)*24,100*(curr50.vsmooth-mean(curr50.vsmooth)),'b--');
plot((curr51.mjd-mjdi)*24,100*(curr51.vsmooth-mean(curr51.vsmooth)),'b-.');
plot((mav40.mjd(mi)-mjdi)*24,mav40.vas(mi)-mean(mav40.vas(mi)),'g');
plot((pred.mjd(predi)-mjdi)*24,pred.v(predi),'k');
axis([0 75 -12 12])
legend('ABE50','ABE51','Tivey')
%legend boxoff
ylabel('v [cm/s]')
xlabel('Time [hr]')

print -depsc2 fm1_v-tseries.pred_n_s_tivey_50_51.eps2
print -djpeg fm1_v-tseries.pred_n_s_tivey_50_51.jpg
end;



if 0,
%% Plot time series figure for thesis
figure(51)
clf;

subplot(2,1,1);
plot((pred.mjd(predi)-mjdi)*24,pred.v(predi),'k');
hold on;
grid
plot((wasp2161.mjd(ni)-mjdi)*24,vn(ni)-mean(vn(ni)),'b');
plot((mab50.mjd(si)-mjdi)*24,vs(si)-mean(vs(si)),'r');
legend('predicted','South','North');
title('Comparison of near-bottom along-axis V');
ylabel('Along-axis v [cm/s]')

subplot(2,1,2);
plot((pred.mjd(predi)-mjdi)*24,pred.v(predi),'k');
hold on;
grid
plot((curr50.mjd-mjdi)*24,curr50.vsmooth-mean(curr50.vsmooth),'b');
plot((curr51.mjd-mjdi)*24,curr51.vsmooth-mean(curr51.vsmooth),'b');
plot((mav40.mjd(mi)-mjdi)*24,mav40.vas(mi)-mean(mav40.vas(mi)),'r');
legend('predicted','ABE50','ABE51','mav40')
ylabel('Along-axis v [cm/s]')
xlabel('Time [hr]')

end;
