%Neil Lareau
%EFD final project 

clear all
close all

% turbulence spectra and statistics


load date4.mat;


%%plot data to get a quick look at it
figure(1)
plot(sonicdata.Ux_4)

%first lets parse out some of the time data


day=sonicdata.Day-90;
hrmin=sonicdata.HRMIN;
sec=sonicdata.SEC;


%%
%split the HRMIN data;
hr=floor(hrmin(:)./100);
mins=((hrmin(:)./100)-hr)*100;
yr=ones(size(mins)).*2011;
mn=ones(size(mins)).*4;


datevector=datenum(yr,mn,day,hr,mins,sec);
datetest=datestr(datevector(1:30),'yy.mm.dd.HH.MM.SS.FFF')

%% manufacture a time series

%time steps in a day
%total number of hours to consider

steps=47*60*60*20;


starttime=datenum(2011,04,10,15,00,00);
endtime=datenum(2011,04,12,13,59,59.95);

syntime=linspace(starttime,endtime,steps);
start=datestr(syntime(1),'yyyy.mm.dd.HH.MM.SS.FFF')

stop=datestr(syntime(end-20:end),'yyyy.mm.dd.HH.MM.SS.FFF')

%% now define some time window of interest

%one hour of data = 72000 samples

%grab all of april 5th
timeindex1=find(hr>=15 & day == 10);
timeindexa=find(day == 11);
timeindex2=find(hr<14  & day ==12);
timeindex=vertcat(timeindex1,timeindexa,timeindex2);


datevector2=interp1(datevector(timeindex),datevector(timeindex),syntime);

ux=interp1(datevector(timeindex),sonicdata.Ux_4(timeindex),datevector2);
uy=interp1(datevector(timeindex),sonicdata.Uy_4(timeindex),datevector2);
uz=interp1(datevector(timeindex),sonicdata.Uz_4(timeindex),datevector2);
t=interp1(datevector(timeindex),sonicdata.Ts_4(timeindex),datevector2);
ux=permute(ux,[2 1]);
uy=permute(uy,[2 1]);
uz=permute(uz,[2 1]);
t=permute(t,[2 1]);


start=datestr(datevector2(1),'yyyy.mm.dd.HH.MM.SS.FFF')
stop=datestr(datevector2(end),'yyyy.mm.dd.HH.MM.SS.FFF')

figure(2)
clf
subplot(3,1,1)
plot(datevector2,ux);
datetick('x',15)
axis tight
subplot(3,1,2)
plot(datevector2,uy);
datetick('x',15)
axis tight
subplot(3,1,3)
plot(datevector2,uz);
datetick('x',15)
axis tight

%% now start in on the coordinate rotation. This begins by finding the 5 minute and then 30 minute mean values


vmean=zeros(length(datevector2),1);
umean=zeros(length(datevector2),1);
wmean=zeros(length(datevector2),1);

%determine the number of  5 minute segments
segs=length(datevector2)./(20*60*5)


for i=1:segs;
    first=((i-1)*6000)+1;
    last=6000*i;
    
    %means
    umean(first:last)=mean(ux(first:last));% 5 minunte means of Ux
    vmean(first:last)=mean(uy(first:last));% 5 minunte means of Ux
    wmean(first:last)=mean(uz(first:last));% 5 minunte means of Ux
    tmean(first:last)=mean(t(first:last));
    
end




%% plot means an full values
figure(3)
clf
subplot(3,1,1)
plot(ux)
hold on 
plot(umean,'r')
subplot(3,1,2)
plot(uy)
hold on 
plot(vmean,'r')
subplot(3,1,3)
plot(uz)
hold on 
plot(wmean,'r')




%% now we can use the means to define the rotation for each 5 minute window period


theta=atan(vmean./umean);

%calculate the wind direction;

v=-umean;
u=vmean;

% figure(200)
% clf
% subplot(2,1,1)
% plot(datevector2,sqrt(u.^2+v.^2),'r')
% hold on
% %plot(datevector2,v,'k');
% datetick('x',15)
% axis tight

thetadir=atand(u./v);
thetadir(v>=0)=thetadir(v>=0)+180;
thetadir(u>=0 & v<0)=thetadir(u>=0 & v<0)+360;


% figure(200)
% clf
% plot(datevector2,thetadir);
% datetick('x',15)
% axis tight
% hold on

%%

u1=ux.*cos(theta)+ uy.*sin(theta);
v1=-ux.*sin(theta)+ uy.*cos(theta);

u1mean=umean.*cos(theta) + vmean.*sin(theta);
v1mean=-umean.*sin(theta)+ vmean.*cos(theta);



figure(4)
clf
plot(u1,'b')
hold on
plot(u1mean,'cyan','linewidth',3)
plot(v1,'r')
plot(v1mean,'m','linewidth',3)
plot(uz-4,'g');
plot(wmean-4,'y','linewidth',3);
legend('Uwind','Vwind','Wwind (offset)')

%% second rotation to set w to zero

phi=atan(wmean./u1mean);

u2=u1.*cos(phi) + uz.*sin(phi);
u2mean=u1mean.*cos(phi) + wmean.*sin(phi);
w2 = -u1.*sin(phi) + uz.*cos(phi);
w2mean=-u1mean.*sin(phi) + wmean.*cos(phi);

figure(5)
clf
plot(u2,'b');
hold on;
plot(u2mean,'c','linewidth',3)
plot(v1,'r')
plot(v1mean,'m','linewidth',3)
plot(w2-2,'g')
plot(w2mean-2,'y','linewidth',3);
legend('Uwind','Vwind','Wwind (offset -2)')
axis tight

%% plot the mean wind speed and direction where umean is the total mean windspeed
figure(200)
clf
subplot(3,1,1)
[ax,h,h1]=plotyy(datevector2,abs(u2),datevector2,thetadir);
hold on;
plot(ax(1),datevector2,abs(u2mean),'c','linewidth',3)
linkprop(ax,{'Xlim','XTickLabel','Xtick'}); % Link the limits of both axes, and the labels
axes(ax(2));datetick('x',15) % Now use DATETICK
set(h,'linewidth',1,'color','b')
set(h1,'linestyle','none','marker','o','color','k')
set(ax(2),'ytick',[0:20:360]);
set(ax(1),'ytick',[0:1:30])
ylabel(ax(1),'Wind Speed')
ylabel(ax(2),'Wind Direction')
xlabel('Time')
%axis tight
title('Wind Speed and Direction')
%hold on;
%plotyy(ax(1),datevector2,u2mean,datevector2,thetadir)

% subplot(3,1,2)
% plot(datevector2,tmean,'r','linewidth',3)
% datetick('x',15)
% ylabel('Temperature [C]');
% xlabel('Time')


%% Now construct the perterbations

uprime=u2-u2mean;
vprime=v1-v1mean;
wprime=w2-w2mean;
tprime=t-permute(tmean,[2 1]);

tmean=permute(tmean,[2 1]);

umean_=u2mean;
u=u2;
v=v1;
vmean=v1mean;
w=w2;
wmean=w2mean;

syntime=syntime';
save('rotated_0410_0412.mat','umean','u','vmean','v','wmean','w','syntime','t','tmean');

%% calculate the dissipation rate by first finding taylors frozen hypthosis values of 

%we need the gradient (or centered differences of uprime in time);

%spacing for time is 20cycles per second or .05 seconds
duprimedt=gradient(uprime,.05);


%% now average over each 30 minutes to come up with the variances and covariances.

%number of 30 minute periods in a day
periods=47*2;

covmatrix30=zeros(3,3,periods);


for i=1:periods;
    first=((i-1)*36000)+1;
    last=36000*i;
    
    %means
    covmatrix30(1,1,i)=mean(uprime(first:last).^2);% 
    covmatrix30(2,2,i)=mean(vprime(first:last).^2);% 
    covmatrix30(3,3,i)=mean(wprime(first:last).^2);
    covmatrix30(1,2,i)=mean(uprime(first:last).*vprime(first:last));    
    covmatrix30(2,1,i)=covmatrix30(1,2,i);
    covmatrix30(1,3,i)=mean(uprime(first:last).*wprime(first:last));
    covmatrix30(3,1,i)=covmatrix30(1,3,i);
    covmatrix30(2,3,i)=mean(vprime(first:last).*wprime(first:last));
    covmatrix30(3,2,i)=covmatrix30(2,3,i);
    kinheatflux(i)=1004.67*mean(wprime(first:last).*tprime(first:last));
    
    dudxsqr(i)=mean(((1./umean(first:last)).*(duprimedt(first:last))).^2);
    
    dateseries(i)=datevector2(last);
    
end

%% Kinematic heat flux 30 minute means

figure(200)
subplot(3,1,2)
%cla(h2)
%cla(h3)
[ax5,h2,h3]=plotyy(dateseries,kinheatflux,datevector2,tmean);
axis tight
linkprop(ax5,{'Xlim','XTickLabel','Xtick'}); % Link the limits of both axes, and the labels
axes(ax5(2));datetick('x',15) % Now use DATETICK
%datetick('x',15)
%axis tight
set(ax5(1),'ytick',-50:10:300)
set(h2,'linewidth',3,'color','m')
set(h3,'Linewidth',3,'color','r')
ylabel(ax5(1),'Watts / m^2')
ylabel(ax5(2),'Temperature [C]')
title('Temperature and Sensible Heat Flux')
legend('Temperature','Senisble Heat Flux')


subplot(3,1,3)
[ax6,h4,h5]=plotyy(dateseries,.5*(squeeze(covmatrix30(1,1,:))+squeeze(covmatrix30(2,2,:))+squeeze(covmatrix30(3,3,:))),dateseries,15*(1.461*10^-5).*(dudxsqr));
axis tight
linkprop(ax6,{'Xlim','XTickLabel','Xtick'}); % Link the limits of both axes, and the labels
axes(ax6(2));datetick('x',15);
set(ax6(1),'ytick',0:.2:1.6)
set(h4,'linewidth',3,'color','c')
set(h5,'Linewidth',3,'color','b')
ylabel(ax6(1),'TKE (m^2/s^2)')
ylabel(ax6(2),'Dissipation Rate [m^2/s^3]');
%set(ax6(1),'ytick',0:1:7);
title('TKE and Dissipation Rate');
legend('Dissipation Rate','TKE')

%% now plot the covariance matrix

figure(6)
clf
subplot(2,1,1)
flatmat=squeeze(covmatrix30(:,:,1))
pcolor(flatmat);
shading interp
caxis([0 max(covmatrix30(:))])
colorbar
title('covariance matrix (first 30 minutes)');
set(gca,'xtick',[1 2 3],'xticklabel',{'u';'v';'w'}) 
set(gca,'ytick',[1 2 3],'yticklabel',{'u';'v';'w'}) 

subplot(2,1,2)
pcolor(covmatrix30(:,:,2));
shading interp
caxis([0 max(covmatrix30(:))])
set(gca,'xtick',[1 2 3],'xticklabel',{'u';'v';'w'}) 
set(gca,'ytick',[1 2 3],'yticklabel',{'u';'v';'w'}) 
colorbar;
title('covariance matrix (second 30 min.)');

%% standard deviation=  sqrt of variance


%% make a plot of the heat flux and see if it makes sense

figure(7)
plot(tprime,'r');
hold on
plot(wprime);

figure(80)
plot(tprime.*wprime,'g')


%% quick fft for each 30 minutes

%plot of uu
uu=uprime.*uprime;
ww=wprime.*wprime;

figure(8)
clf
plot(uu);




%% best approach is to use a hanning filter and the welch method for a sliding window to make the stats more robust

%to get the TKE spectra we must add the spectra for uu, vv, and ww
%number of 30 minute periods in a day
periods=47*2;

pu=zeros(36000,periods);
%pv=zeros(1801,periods);
pw=zeros(36000,periods);
pt=zeros(36000,periods);
puu=zeros(size(pu));
%pvv=zeros(size(pu));
%pww=zeros(size(pu));
%puv=zeros(size(pu));
%puw=zeros(size(pu));
%pvw=zeros(size(pu));
%pwt=zeros(size(pu));
%co=zeros(18000,periods);
%Puu=zeros(size(pu));
for i=1:periods;
    first=((i-1)*36000)+1;
    last=36000*i;
    L=length(uprime(first:last))./10;

    %do it right
    test1=fft(wprime(first:last));
    test2=fft(tprime(first:last));
    coa(:,i)=real(test1(1:18000)).*real(test2(1:18000)) + (imag(test1(1:18000).*imag(test2(1:18000))));
    co(:,i)=test1(1:18000).*conj(test2(1:18000)) + (test2(1:18000).*conj(test1(1:18000)));
    
    %uu
    pu=fft(uprime(first:last));
    %puu(:,i)=pu(1:1801).*conj(pu(1:18001)) + pu(1:1801).*conj(pu(1:1801));
    puu(:,i)=2*(abs(pu(2:length(pu)/2).^2));
    %pv=fft(vprime(first:last));
    %pvv(:,i)=pv(1:1801).*conj(pv(1:1801)) + pv(1:1801).*conj(pv(1:1801));
    pw=fft(wprime(first:last));
    %pww(:,i)=pw(1:1801).*conj(pw(1:1801)) + pw(1:1801).*conj(pw(1:1801));
    %puv(:,i)=pu(1:1801).*conj(pv(1:1801)) + pv(1:1801).*conj(pu(1:1801));
    %puw(:,i)=pu(1:1801).*conj(pw(1:1801)) + pw(1:1801).*conj(pu(1:1801));
    %pvw(:,i)=pv(1:1801).*conj(pw(1:1801)) + pw(1:1801).*conj(pv(1:1801));
    
    %pwelch puu, pvv, pww terms
    [Puu(:,i),f(:,i)]=pwelch(uprime(first:last),hanning(L),round(L/2),L,20);
    

end




%%


figure(100)
clf
pcolor(puw(2:end,:)')
shading flat
set(gca,'xscale','log','ydir','rev')

%coa(coa>0)=log(coa(coa>0));
%coa(coa<0)=-log(-coa(coa<0));%flip the sign, take the log, flip it back;
figure(102)
clf
pcolor(coa(2:end,:)')
shading flat
set(gca,'xscale','log','ydir','rev')
figure(101)
clf
plot(coa(2:end,63)');
axis tight
set(gca,'xscale','log')

coa(1,63)


%% TKE spectra for each 30 minute period for April 5th
figure(209)
clf
pcolor(log((puu(2:end,:)')));
shading flat
set(gca,'xscale','log','ydir','rev')

figure(210)
clf
pcolor(log((pvv(2:end,:)')));
shading flat
set(gca,'xscale','log','ydir','rev')

figure(211)
clf
pcolor(log((pww(2:end,:)')));
shading flat
set(gca,'xscale','log','ydir','rev')


figure(212)
clf
loglog(puu(2:end,63))

figure(213)
clf
loglog(Puu(2:end,63))


%figure(11)
%clf
%pcolor(log(pwt'))
%shading flat;
%set(gca,'xscale','log','ydir','rev')

%% look at the spectra for the overnight period






