% make_currfigs.m
%
% 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
load_wasp
load_cm95
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;

% 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);

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;

% Use Matlab PSD and CSD functions to get spectral power densities G11,G22 
% and cross spectral power denisty G12
%[Pxx, Pxxc, F] = PSD(X,NFFT,Fs,WINDOW,NOVERLAP,P)
% Note size of v1int and v2int is 1775x1
% he default values for the parameters are NFFT = 256 (or LENGTH(X), whichever is smaller), NOVERLAP = 0, WINDOW = HANNING(NFFT), Fs = 2, P = .95, and DFLAG = 'none'.
lengthseries=length(v1int);
nfft=256
windowwidth=256
fs=1;
over=0;  % ~30% overlap
[g11,g11conf,f11]=psd(v1int,nfft,fs,windowwidth,over);
[g22,g22conf,f22]=psd(v2int,nfft,fs,windowwidth,over);
[g12,g12conf,f12]=csd(v1int,v2int,nfft,fs,windowwidth,over);

% zoom in on tidal and subtidal frequencies
flow=find(f12/fs<=0.1);
lowf=f12(flow); 
lowp=phase(flow);
%NOTE: Hanning window applied to entire data series yields
% EDOF=(length(timeseries)/windowhalfwidth)*8/3=5.33
% So, 95% confidence level = gamma_95 = 1-(0.05)^(1/(5.33-1))
edof=(8/3)*(length(v1int)/windowwidth/2)
gamma_95=1-(0.05)^(1/(edof-1))  % = 0.499
c95=find(coh(flow)>gamma_95);

figure(40)
clf;
subplot(2,1,1)
semilogy(lowf,g11(flow),'b');
hold on;
semilogy(lowf,g22(flow),'r-.');
%axis([0 0.1 1e-5 1e0]);
%set(gca,'YTick',[1e-5 1e0])
grid on;
title('Power spectra of along-axis current component')
ylabel('spectral density')
legend('N mooring', 'S mooring')

% calculate coherence and phase between S and N meter
coh= abs(g12).^2 ./ (g11.*g22);
phase= atan2(-imag(g12),real(g12)) * 180/pi;

subplot(2,1,2)
%set(gca,'Position',[0.1300    0.1100    0.7750    0.55]);
[ax,hc,hp]=plotyy(lowf/fs,coh(flow),lowf(c95)/fs,lowp(c95));
set(hc,'Color',[0 0 1],'LineStyle','-');
set(hp,'MarkerFaceColor',[0 1 0],'LineStyle','none','marker','*');
hold on;
grid
plot([min(f12) max(f12)],[gamma_95 gamma_95],'k','LineWidth',[1])
title('Coherence and phase (>95%), S vs N moorings')
xlabel('Frequency, hr^-^1')
set(get(ax(1),'YLabel'),'string','Coherence^2');
set(get(ax(2),'YLabel'),'string','Phase');
y1tls=['0.0';'0.2';'0.4';'0.6';'0.8';'1.0'];
set(ax(1),'YTick',[0:0.2:1],'YTickLabel',y1tls);
y2tls=['-180';'-135';' -90';' -45';'   0';'  45';'  90';' 135';' 180'];
set(ax(2),'YTick',[-180:45:180],'YTickLabel',y2tls);
text(0.005,.95,sprintf('N=%g; NFFT=%g; Windowwidth=%g; EDOF=%g; overlap=%g',length(v1int),nfft,windowwidth,edof,over));
print -depsc2 power_coh_phase_NS.eps2 


% Create time-series figure
%

% load or rename predicted, meg, N (wasp2000) and S (fm2000/ois) 
un=u2; vn=v2;
us=u1; vs=v1;
load /shuksan/scottv/snap/flowmow/currents/predicted/pred.mat
load /shuksan/scottv/snap/flowmow/currents/tivey/meg_data.mat

umeg=;
vmeg=;

% transform to along/across axis components

% get the abe50 interval for comparison
time_table;
mjdi=tm.strd50+tm.mor150(1)/86400;
mjdf=tm.strd50+tm.top250(2)/86400;

figure(31)
clf
%subplot 2,1,1;
hold on;
grid;
plot(pred.mjd(predind),pred.v(predind),'k'); 
plot(mab50.day(mab50ind),mab50.v(mab50ind)*100,'r'); 
plot(wasp2161.mjd(waspind),wasp2161.v(waspind)*100,'m-.');
plot(meg.mjd(megind),meg.v(megind),'g');
plot(curr.mjd,curr.v*100,'c--');
%plot(curr.mjd,curr.v3*100,'b--');
legend('predicted','FM50','WASP','meg','ABE50')
title('Comparison of V from FM50mab, WASP, Tivey MAVS, ABE50, and predicted currents, and pressure')
ylabel('N/S Velocity, V [cm/s]')
xlabel('MJD')
relabelx;

if 0,
subplot 2,1,2;
hold on;
grid;
plot(pres.mjd(MPind),pres.dbars(MPind),'g');
xlabel('MJD')
relabelx;

end;

print -depsc2 grand_comp_abe50.eps2
