clear;

load /sirius2/zhup/wrk/tower_data/Everglads/out_10m_dat_1.mat;

%---12:00-01:15
u_d=u(:,1:5);
v_d=v(:,1:5);
w_d=w(:,1:5);

%---03:45-05:00
u_n=u(:,64:68);
v_n=v(:,64:68);
w_n=w(:,64:68);

u_al_d=u_d*cos(30*pi/180)-v_d*sin(30*pi/180);
v_cr_d=u_d*sin(30*pi/180)+v_d*cos(30*pi/180);
u_al_n=u_n*cos(30*pi/180)-v_n*sin(30*pi/180);
v_cr_n=u_n*sin(30*pi/180)+v_n*cos(30*pi/180);


std_w_n=std(w_n,0,1);
std_w_d=std(w_d,0,1);
std_v_n=std(v_cr_n,0,1);
std_v_d=std(v_cr_d,0,1);
std_u_n=std(u_al_n,0,1);
std_u_d=std(u_al_d,0,1);

for nn=1:4
w1=squeeze(w_n(:,nn));

N=length(w1);

dt=0.05;
df = 1/(dt*N);

nyq=floor(N/2);
Period=N*dt;
for n=1:nyq
f(n)=n/Period;
end

%-----window
window=blackhar(N)';
DNT2=fft(window.*detrend(w1)',N)/N;
PNT2=DNT2(2:N).* conj(DNT2(2:N));
if ((N/2)-floor(N/2)) == 0
PDS2(1:nyq-1)=2*PNT2(1:nyq-1);
PDS2(nyq)=PNT2(nyq);
else
PDS2(1:nyq)=2*PNT2(1:nyq);
end
PDS2_w_n(:,nn)=PDS2/df./(norm(window)^2/N); % power density S

end; %nn

for nn=1:4
w1=squeeze(w_d(:,nn));

N=length(w1);

dt=0.05;
df = 1/(dt*N);

nyq=floor(N/2);
Period=N*dt;
for n=1:nyq
f(n)=n/Period;
end

%-----window
window=blackhar(N)';
DNT2=fft(window.*detrend(w1)',N)/N;
PNT2=DNT2(2:N).* conj(DNT2(2:N));
if ((N/2)-floor(N/2)) == 0
PDS2(1:nyq-1)=2*PNT2(1:nyq-1);
PDS2(nyq)=PNT2(nyq);
else
PDS2(1:nyq)=2*PNT2(1:nyq);
end
PDS2_w_d(:,nn)=PDS2/df./(norm(window)^2/N); % power density S

end; %nn

%------------------------
for nn=1:4
w1=squeeze(u_al_n(:,nn));
%w1=squeeze(u_n(:,nn));

N=length(w1);

dt=0.05;
df = 1/(dt*N);

nyq=floor(N/2);
Period=N*dt;
for n=1:nyq
f(n)=n/Period;
end

%-----window
window=blackhar(N)';
DNT2=fft(window.*detrend(w1)',N)/N;
PNT2=DNT2(2:N).* conj(DNT2(2:N));
if ((N/2)-floor(N/2)) == 0
PDS2(1:nyq-1)=2*PNT2(1:nyq-1);
PDS2(nyq)=PNT2(nyq);
else
PDS2(1:nyq)=2*PNT2(1:nyq);
end
PDS2_u_n(:,nn)=PDS2/df./(norm(window)^2/N); % power density S

end; %nn

for nn=1:4
w1=squeeze(u_al_d(:,nn));
%w1=squeeze(u_d(:,nn));

N=length(w1);

dt=0.05;
df = 1/(dt*N);

nyq=floor(N/2);
Period=N*dt;
for n=1:nyq
f(n)=n/Period;
end

%-----window
window=blackhar(N)';
DNT2=fft(window.*detrend(w1)',N)/N;
PNT2=DNT2(2:N).* conj(DNT2(2:N));
if ((N/2)-floor(N/2)) == 0
PDS2(1:nyq-1)=2*PNT2(1:nyq-1);
PDS2(nyq)=PNT2(nyq);
else
PDS2(1:nyq)=2*PNT2(1:nyq);
end
PDS2_u_d(:,nn)=PDS2/df./(norm(window)^2/N); % power density S

end; %nn

%------------------------
for nn=1:4
w1=squeeze(v_cr_n(:,nn));
%w1=squeeze(v_n(:,nn));

N=length(w1);

dt=0.05;
df = 1/(dt*N);

nyq=floor(N/2);
Period=N*dt;
for n=1:nyq
f(n)=n/Period;
end

%-----window
window=blackhar(N)';
DNT2=fft(window.*detrend(w1)',N)/N;
PNT2=DNT2(2:N).* conj(DNT2(2:N));
if ((N/2)-floor(N/2)) == 0
PDS2(1:nyq-1)=2*PNT2(1:nyq-1);
PDS2(nyq)=PNT2(nyq);
else
PDS2(1:nyq)=2*PNT2(1:nyq);
end
PDS2_v_n(:,nn)=PDS2/df./(norm(window)^2/N); % power density S

end; %nn

for nn=1:4
w1=squeeze(v_cr_d(:,nn));
%w1=squeeze(v_d(:,nn));

N=length(w1);

dt=0.05;
df = 1/(dt*N);

nyq=floor(N/2);
Period=N*dt;
for n=1:nyq
f(n)=n/Period;
end

%-----window
window=blackhar(N)';
DNT2=fft(window.*detrend(w1)',N)/N;
PNT2=DNT2(2:N).* conj(DNT2(2:N));
if ((N/2)-floor(N/2)) == 0
PDS2(1:nyq-1)=2*PNT2(1:nyq-1);
PDS2(nyq)=PNT2(nyq);
else
PDS2(1:nyq)=2*PNT2(1:nyq);
end
PDS2_v_d(:,nn)=PDS2/df./(norm(window)^2/N); % power density S

end; %nn

%------------------------

m_PDS2_w_d=squeeze(mean(PDS2_w_d,2));
m_PDS2_w_n=squeeze(mean(PDS2_w_n,2));
m_PDS2_v_d=squeeze(mean(PDS2_v_d,2));
m_PDS2_v_n=squeeze(mean(PDS2_v_n,2));
m_PDS2_u_d=squeeze(mean(PDS2_u_d,2));
m_PDS2_u_n=squeeze(mean(PDS2_u_n,2));

m_std_w_d=mean(std_w_d);
m_std_w_n=mean(std_w_n);
m_std_v_d=mean(std_v_d);
m_std_v_n=mean(std_v_n);
m_std_u_d=mean(std_u_d);
m_std_u_n=mean(std_u_n);

mm=19;

save data/spectrum_d_n.mat f m_PDS2_w_d m_PDS2_w_n m_PDS2_v_d m_PDS2_v_n m_PDS2_u_d m_PDS2_u_n m_std_w_d m_std_w_n m_std_v_d m_std_v_n m_std_u_d m_std_u_n;

subplot(2,2,1)
h=loglog(f,m_PDS2_w_n/m_std_w_n^2,'b');set(h,'linewidth',0.5);
hold on;
h=loglog(f,m_PDS2_w_d/m_std_w_d^2,'r');set(h,'linewidth',0.5);
x=[0.001 10.0];
y=[1000.0 0.0001];
loglog(x,y,'g--');
xlabel('Frequency (1/s)');
ylabel('Power density (m^2/s^3)');
hold off;

subplot(2,2,2)
h=loglog(f,m_PDS2_v_n/m_std_v_n^2,'b');set(h,'linewidth',0.5);
hold on;
h=loglog(f,m_PDS2_v_d/m_std_v_d^2,'r');set(h,'linewidth',0.5);
x=[0.001 10.0];
y=[1000.0 0.0001];
loglog(x,y,'g--');
xlabel('Frequency (1/s)');
ylabel('Power density (m^2/s^3)');
hold off;

subplot(2,2,3)
h=loglog(f,m_PDS2_u_n/m_std_u_n^2,'b');set(h,'linewidth',0.5);
hold on;
h=loglog(f,m_PDS2_u_d/m_std_u_d^2,'r');set(h,'linewidth',0.5);
x=[0.001 10.0];
y=[1000.0 0.0001];
loglog(x,y,'g--');
xlabel('Frequency (1/s)');
ylabel('Power density (m^2/s^3)');
hold off;


