% Computes U relative to ocean surface
% Assumes HiWinGS_2013_10minMeansNFluxes_r1.mat is loaded as
% array MMF, containing cspd and cdir @ 10 min resolution.
% Also assumes true U and dir exist in workspace

CN = -MMF.cspd10m(ib).*cos((MMF.cdir10m(ib))./180.*pi);  %'toward' North
CE = -MMF.cspd10m(ib).*sin((MMF.cdir10m(ib))./180.*pi);  %'toward' East
% UN = MMF.U10.*cos(MMF.wdir./180.*pi);
% UE = MMF.U10.*sin(MMF.wdir./180.*pi); 
% UrN = MMF.Ur10.*cos(MMF.wdirR./180.*pi);
% UrE = MMF.Ur10.*sin(MMF.wdirR./180.*pi); 
% 
% figure;subplot(2,1,1);plot(MMF.yday,UrN-UN,'r.',MMF.yday,CN,'g');title('Northward current speed');grid;axis([274 303 -.5 .5]);legend('UrN-UN','CN','Location','SouthWest')
% subplot(2,1,2);plot(MMF.yday,UrE-UE,'r.',MMF.yday,CE,'g');title('Eastward current speed');grid;axis([274 303 -0.1 1]);legend('UrE-UE','CE','Location','SouthWest')
rdcon = pi/180;
% Ul=sqrt((UN+CN).^2+(UE+CE).^2);figure;plot(MMF.yday,MMF.Ur10,'k',MMF.yday,Ul,'r');figure;plot(MMF.yday,MMF.Ur10-Ul,'r')
% dirl=mod(atan2((UE+CE),(UN+CN))/rdcon+360,360);figure;plot(MMF.yday,MMF.wdirR,'k',MMF.yday,dirl,'r');figure;plot(MMF.yday,MMF.wdirR-dirl,'r')

UN = U.*cos(dir./180.*pi); % North
UE = U.*sin(dir./180.*pi); % East

Ul = U; dirl = dir;
Ul = sqrt((UN+CN).^2+(UE+CE).^2);  % Wind speed (m/s) relative to water
dirl  =mod(atan2((UE+CE),(UN+CN))/rdcon+360,360);  % Wind direction (deg) from relative to water

figure; plot(MMF.yday10m,MMF.Ur1010m,'k',jdy,Ul,'r');
figure; plot(MMF.yday10m,MMF.wdirR10m,'k',jdy,dirl,'r');

% ucN = MMF.cspd.*cos((MMF.cdir)./180.*pi);
% ucE = MMF.cspd.*sin((MMF.cdir)./180.*pi);  
% ucNl = MMF.cspd.*cos((MMF.cdir+180)./180.*pi);
% ucEl = MMF.cspd.*sin((MMF.cdir+180)./180.*pi);  


