close all
clear all
% ***** read u
load u_v_time_4nodes_re1000.dat
fil=u_v_time_4nodes_re1000;
% plot u
u1=fil(:,1);
v1=fil(:,2);
u2=fil(:,3);
v2=fil(:,4);
u3=fil(:,5);
v3=fil(:,6);
u4=fil(:,7);
v4=fil(:,8);

% time step
dt=8.177E-04;
n=length(u1);
% compute time array
t=dt:dt:n*dt;

umean=mean(u3);

% subtract the mean
u3_fluct=u3-umean;

% number of points in the fft

% compute RMS
n=length(u3_fluct);
urms=0;
for i=1:n
   urms=urms+u3_fluct(i)^2/n;
end
urms=urms^0.5;

% f_morm is normalized frequency [0 pi]
[pu3 f_norm]=pwelch(u3_fluct);
% the true frequency
f= f_norm/pi/dt;

plot(f,pu3,'b-','linewidth',2)
hold

% add line with -5/3 slope
x=[1 100];
ynoll=1;
y(1)=ynoll;
y(2)=y(1)*(x(2)/x(1))^(-5/3);
plot(x,y,'r--','linew',4)


axis([0.1 1500 1e-6 100])
title('spectrum of u','fontsize',20)
xlabel('f [1/s]','fontsize',20)
ylabel('E(u) [m^2/s]','fontsize',20)
handle=gca;
set(handle,'yscale','log','xscale','log')
xticks([1 10 100 1000])
%set(h,'xscale','log')
set(handle,'fontsi',[20])
print u_spectum.ps -deps
%



