load winds_jfk_low;
tm1=tm;
load eta_bat;
tm2=tm;
H=hlow;
dt=1/24;
time=max(tm1(1),tm2(1)):dt:min(tm1(end),tm2(end));
ue=interp1(tm1,ue,time);
un=interp1(tm1,un,time);
H=interp1(tm2,H,time);
i=sqrt(-1);
U=ue+i*un;
c=corrcoef(U,H)
% Rotate winds
rot=angle(c(2,1));
U=U*exp(-i*rot);
subplot(2,1,1);
u=real(U);
plot(time,H/std(H),time,u/std(u));
datetick('x',6,'keepticks');


% Estimate time decorelation time scale;
lag1=0;
lag2=72;
for ilag=lag1+1:lag2+1;
    X=u(1:end-ilag+1);
    Y=u(ilag:end);
    c=corrcoef(X,Y);
    r1(ilag)=c(2,1);
end

for ilag=lag1+1:lag2+1;
    X=H(1:end-ilag+1);
    Y=H(ilag:end);
    c=corrcoef(X,Y);
    r2(ilag)=c(2,1);
end

subplot(2,1,1);
plot(lag1:lag2,r1,'b',lag1:lag2,r2,'r');
grid on





% Find  lagged correlation
subplot(2,1,2);
for ilag=lag1+1:lag2+1;
    X=u(1:end-ilag+1);
    Y=H(ilag:end);
    c=corrcoef(X,Y);
    r(ilag)=c(2,1);
end
    
plot(lag1:lag2,r);
grid on