%Neil Lareau
%EFD final project 

clear all
close all

% turbulence spectra and statistics


load redbutte1.mat;
%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
steps=24*60*60*20;


starttime=datenum(2011,04,04,20,00,00);
endtime=datenum(2011,04,05,19,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>=20 & day == 4);
timeindex2=find(hr<20  & day ==5);
timeindex=vertcat(timeindex1,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]);

%what day/time does this start/end

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

figure(2)
subplot(3,1,1)
plot(ux);
subplot(3,1,2)
plot(uy);
subplot(3,1,3)
plot(uz);


%% 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;
thetadir=atand(vmean./umean);

figure(10)
plot(thetadir);


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


%% Now construct the perterbations

uprime=u2-u2mean;
vprime=v1-v1mean;
wprime=w2-w2mean;
tprime=t-permute(tmean,[2 1]);
%% now average over each 30 minutes to come up with the variances and covariances.

%number of 30 minute periods in a day
periods=48;

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)=mean(wprime(first:last).*tprime(first:last));
end


%% 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=48;

pu=zeros(1801,periods);
pv=zeros(1801,periods);
pw=zeros(1801,periods);
pt=zeros(1801,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));

for i=1:periods;
    first=((i-1)*36000)+1;
    last=36000*i;
    L=length(uprime(first:last))./10;
    pu(:,i)=pwelch(uprime(first:last),hanning(L),round(L/2),L,20);
    pv(:,i)=pwelch(vprime(first:last),hanning(L),round(L/2),L,20);
    pw(:,i)=pwelch(wprime(first:last),hanning(L),round(L/2),L,20);
    pt(:,i)=pwelch(tprime(first:last),hanning(L),round(L/2),L,20);
    puu(:,i)=pu(:,i).*pu(:,i);
    pvv(:,i)=pv(:,i).*pv(:,i);
    pww(:,i)=pw(:,i).*pw(:,i);
    puv(:,i)=pu(:,i).*pv(:,i);
    puw(:,i)=pu(:,i).*pw(:,i);
    pvw(:,i)=pv(:,i).*pw(:,i);
    pwt(:,i)=pw(:,i).*pt(:,i);

end



%% TKE spectra for each 30 minute period for April 5th
figure(9)
clf
pcolor(log(puu'+pvv'+pww'));
shading flat
set(gca,'xscale','log','ydir','rev')

figure(10)
clf
loglog(puu(:,end))

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


%% Frozen hypothesis

%here whe relate time to space via t=x/U where u is the mean wind speed (in
%the rotated frame and t is the time step...

%solving for x; x=t/U

%to keep the spacing of x constant we can only use one 30 minute averaging
%period with a constant mean spead.


%we need to recaculate uprime for the first 30 minute window by first
%averaging the 5 minute means into one 30 minute mean

umeanfrozen=mean(umean(1:36000));
uprimefrozen=u2(1:36000)-umeanfrozen;

%now peform the FFT on the frozen values
[pufrozen,f]=pwelch(uprimefrozen(1:36000),hanning(L),round(L/2),L,20);

figure(12)
loglog(f,pufrozen)

%define a new xaxis in terms of wave number
xsize=f./umeanfrozen;

figure(13)
loglog(xsize,pufrozen)

