clear all; close all;
set (gca,'fontsize',16,'fontname','Times New Roman');  

% The choice of 'a' corresponds to the choice of pitching direction as
% applied in G. R. Pirrung and H. Aa. Madsen, 'Dynamic inflow effects in
% measurements and high fidelity computations'
a=0.5;
%a=0.0;

% velocity at hub normalized by radius!
% because radius = 5 m and v_free= 5 m/s
V_tr=(1.0-1.5*a);

xRmax=20;
%% Compute downwash due to cylindrical vortex sheet 

phi=0:0.01:2*pi;
rR=[0.3 0.47 0.63 0.80 0.95];
dx=0.001;
xR=0:dx:xRmax;

% Numerical integration of the induced velocity due to a cylindrical wake
for l=1:length(xR)
    for i=1:length(rR)
        phi_int(i,l)=(trapz(phi,(1-rR(i).*cos(phi))./(xR(l)^2+1+rR(i)^2-2*rR(i).*cos(phi)).^(3/2)));
    end
end

for l=2:length(xR)
    for i=1:length(rR)
        ui(i,l)=trapz(xR(1:l),phi_int(i,1:l));
    end
end
ui(:,1)=0;

for i=1:length(rR)
    ui_norm(i,:)=ui(i,:)./ui(i,length(xR));
end

set (gca,'fontsize',16,'fontname','Times New Roman'); 
 plot(xR, ui_norm(1,:),'-k',xR, ui_norm(2,:),'-b',xR, ui_norm(3,:),'-g',xR, ui_norm(4,:),'-r',xR, ui_norm(5,:),'-m','LineWidth',2)
 xlim([0 2.5]);
hleg=legend('0.3 R','0.47 R', '0.63 R', '0.8 R', '0.95 R','Location', 'Southeast');
xlabel('x/R [-]')
line([0.35 0.35], [0.0 1.0],'Color','black','LineStyle','--')
line([0.35/4.0 0.35/4.0], [0.0 1.0],'Color','black','LineStyle','--')
ylabel('ui/ui_\infty [-]')
xlim([0 2.5])
ylim([0.0 1.0])
 %%
 % find the time constants according to the values of induced velocities
 % ECN did it to the forces, assuming they behave the same way in time as
 % the induced velocities
for i=1:length(rR)
    for l=1:length(xR);
       tau_en(i,l)=-(xR(l))/log(1-ui_norm(i,l));
    end
end
tau_en(:,:)=tau_en(:,:);
figure;
plot(xR, tau_en(1,:), '-k', xR, tau_en(2,:), '-b', xR, tau_en(3,:), '-g', xR, tau_en(4,:), '-r', xR, tau_en(5,:), '-m','LineWidth',2)
set (gca,'fontsize',16,'fontname','Times New Roman'); 
hleg=legend('0.3 R','0.47 R', '0.63 R', '0.8 R', '0.95 R','Location','Southeast');
xlabel('x/R [-]')
ylabel('\tau [s]')
xlim([0 2.5])
ylim([0.0 1.0])
%% what if we only look at the induced velocity from 0.35 seconds onwards?
t_start=0.35;
for i=1:length(rR)
    ui_norm_tstart(i,:)=ui(i,:)-ui(i,floor(t_start*V_tr/dx));
    ui_norm_tstart(i,:)=ui_norm_tstart(i,:)./ui_norm_tstart(i,length(xR));
end
figure;
set (gca,'fontsize',20,'fontname','Times New Roman'); 
plot(xR/V_tr, ui_norm_tstart(1,:),'-k',xR/V_tr, ui_norm_tstart(2,:),'-b',xR/V_tr, ui_norm_tstart(3,:),'-g',xR/V_tr, ui_norm_tstart(4,:),'-r',xR/V_tr, ui_norm_tstart(5,:),'-m','LineWidth',2)
hleg=legend('0.3 R','0.47 R', '0.63 R', '0.8 R', '0.95 R','Location', 'Southeast');
xlabel('Time [s]')
ylabel('(ui-ui(t=0.35s)/(ui_\infty-ui(t=0.35s)) [-]')
xlim([0.0 2.5])
if a==0.5 
    ylim([0.0 0.7]);
elseif a==0.0
    ylim([0.0 1.0]);
end