clear all;close all; clc
tic
%%%Volba hodnot základních parametrů%%%
% V_out=        náběžná rychlost větru; [m/s]
% lambda=       TSR; [-]
% R=            poloměr rotoru; [m]
% c=            délka tětivy profilu; [m]
% N=            počet lopatek; [-]
% h=            výška rotoru(výška lopatky);[m]
% ro=           hustota vzduchu; [kg/m^3]
% ny=           kinematická viskozita vzduchu; [m^2/s]
% alfa0=        úhel nastavení; [°]

%%%Načtení charakteristik profilu%%%
%Načtení matice charakteristik porfilu%
% profil=importdata('');    matice charakteristik porfilu
% x=[(profil(:,1))];        vektor úhlů náběhu                  
% y=[(profil(:,2))];        vektor CL
% yy=[(profil(:,3))];       vektor CD

%Použití interpolační funkce - načtení databáze CL,CD,Re profilu (kapitola 7.4.1)%
% Avektor=load('');   vektor CL
% Bvektor=load('');   vektor Re
% Cvektor=load('');   vektor CD
  
%%%Volba parametrů ovlivňující časovou náročnost výpočtu a konvergenci%%%
% delta_theta=  krok azimutálního úhlu [°]
% presnost=     parametr přesnost [-]
% r=            relaxacni faktor pro upstream
% r_dw=         relaxacní faktor pro downstream

%%%Výpočet ostatních parametrů%%%
sigma=(N*c)/(2*R);          %hustota lopatkové mříže [-]
omega=(lambda*V_out)/R;     %úhlová rychlost [1/s]
Vobv=omega*R;               %obvodová ryclhost [m/s]
otacky=(omega*30)/(pi);     %otacky [1/s]
Re_ave=(c*omega*R)/ny;      %průměrné Re [-]

%%%Rozdělení prostoru rotoru%%%
alfa0_rad=deg2rad(alfa0);
delta_theta_rad=deg2rad(delta_theta);
A=((-pi/2)+delta_theta_rad);
B=((pi/2)-delta_theta_rad);
theta=A:delta_theta_rad:B;
theta_deg=rad2deg(theta);
C=((0.5*pi)+delta_theta_rad);
D=((1.5*pi)-delta_theta_rad);
thetadw=C:delta_theta_rad:D;
thetadw_deg=rad2deg(thetadw);

%%%algoritmus%%%
u_new=ones(size(theta));    %Inicializace vektoru počátečních hodnot faktoru u
chyba=ones(size(theta));    %Inicializace vektoru pro kontrolu konvergence výpočtu
counter=0;                  %Počet iterací
CFx1=1.816;                 %Koeficient přímky Glauertovi korekce
CFxT=2*sqrt(CFx1)-CFx1;     %Přechodná hodnota koeficientu tahu                     %rovnice [46]
ut=0.5*sqrt(CFx1);          %Přechodná hodnota faktoru u

%%%Výpočet proudových trubic v sekci upstream%%%
    for i = 1:numel(theta)
disp(['theta: ',num2str(theta(i))])    
    while ( chyba(i) > presnost) 
    u(i)=u_new(i);
    Vau(i)=u(i).*V_out;                                                             %rovnice [17]
    Wup(i)=V_out.*sqrt((lambda+u(i).*sin(theta(i))).^2+(u(i).*cos(theta(i))).^2);   %rovnice [24]
    alpha_up(i)=atan((u(i).*cos(theta(i)))./(lambda+u(i).*sin(theta(i))));          %rovnice [23]
    alpha_up_deg(i)=rad2deg(alpha_up(i));                                           
    Re_up(i)=(Wup(i).*c)./ny;                                                       %rovnice [61]
    
    CL(i)=interp1(x,y,alpha_up_deg(i)-alfa0);                                       %interpolace charakteristik profilu
    CD(i)=interp1(x,yy,alpha_up_deg(i)-alfa0);                                      %interpolace charakteristik profilu
   
%   [CL(i),CD(i)]=getCLCD(Re_up(i),alpha_up_deg(i)-alfa0,Avektor,Bvektor,Cvektor)   %při uvažování interpolační funkce    
    
    CN(i)=CL(i).*cos(alpha_up(i))+CD(i).*sin(alpha_up(i));                                  %rovnice [25]
    CT(i)=CL(i).*sin(alpha_up(i))-CD(i).*cos(alpha_up(i));                                  %rovnice [26]
    CFa(i)=-CT(i).*sin(theta(i))+CN(i).*cos(theta(i));                                      %rovnice [27]
    CFx(i)=(1/(2*pi)).*(((Wup(i)./V_out).^2).*((N*c)./(R)).*(CFa(i)./abs(cos(theta(i)))));  %rovnice [31]
    
      if (CFx(i)>CFxT)                              %Glauertova korekce
      u_new(i)=(-CFx1+CFx(i))/(-4*(sqrt(CFx1)-1));                                          %rovnice [47]
      else                                          %Hybnostní teorie
      u_new(i)=0.5+0.5*sqrt(1-CFx(i));                                                      %rovnice [34]
      end
      
    chyba(i)=abs(u(i)-u_new(i));                    %Kontrola konverengece výpočtu
        if chyba(i)>presnost                        
            u_new(i)=u(i)+r.*(u_new(i)-u(i));       %Použití relaxačního faktoru            %rovnice [63]
        end
    disp(u(i))
    counter=counter+1
    end
        u_new(i)=u(i);
        i=i+1;                                      %Výpočet následující produové trubice
    end     
    
 %%%Výpočet proudových trubic v sekci downstream%%%
u_new_flip=flip(u_new);     %Převrácení vektoru u_new kvůli geometrii
udw_new=flip(u_new);        %Inicializace vektoru počátečních hodnot faktoru u
counterdw=0;                %Počet iterací                 
chybadw=ones(size(theta));  %Inicializace vektoru pro kontrolu konvergence výpočtu
  for j = 1:numel(thetadw)
      disp(['thetadw: ',num2str(thetadw(j))])
  while ( chybadw(j) > presnost)
    udw(j)=udw_new(j);
    Vdw(j)=udw(j).*(2.*u_new_flip(j)-1).*V_out;                                                                                         %rovnice [18;22]
    Wdw(j)=(2.*u_new_flip(j)-1).*V_out.*sqrt((((lambda)/(2*u_new_flip(j)-1))+udw(j).*sin(thetadw(j))).^2+(udw(j).*cos(thetadw(j))).^2); %rovnice [36]
    alpha_dw(j)=atan((udw(j).*cos(thetadw(j)))./(((lambda)/(2*u_new_flip(j)-1))+udw(j).*sin(thetadw(j))));                              %rovnice [35]
    alpha_dw_deg(j)=rad2deg(alpha_dw(j));                           
    Re_dw(j)=(Wdw(j).*c)./ny;                                                                                                           %rovnice [62]
    
    CLdw(j)=interp1(x,y,alpha_dw_deg(j)-alfa0);                                     %interpolace charakteristik profilu
    CDdw(j)=interp1(x,yy,alpha_dw_deg(j)-alfa0);                                    %interpolace charakteristik profilu
    
%   [CLd(j),CDd(j)]=getCLCD(Re_dw(j),alpha_dw_deg(j)-alfa0,Avektor,Bvektor,Cvektor) %při uvažování interpolační funkce 
   
    CNdw(j)=+CLdw(j).*cos(alpha_dw(j))+CDdw(j).*sin(alpha_dw(j));                                   %rovnice [37]
    CTdw(j)=CLdw(j).*sin(alpha_dw(j))-CDdw(j).*cos(alpha_dw(j));                                    %rovnice [38]
    CFadw(j)=-CTdw(j).*sin(thetadw(j))+CNdw(j).*cos(thetadw(j));                                    %rovnice [39]
    V_e(j)=(2.*u_new_flip(j)-1).*V_out;                                                             %rovnice [22]
    CFxdw(j)=(1/(2*pi)).*(((Wdw(j)./V_e(j)).^2).*((N*c)./(R)).*(CFadw(j)./abs(cos(thetadw(j)))));   %rovnice [42]
 
    if (CFxdw(j)>CFxT)                                      %Glauertova korekce
        udw_new(j)=(-CFx1+CFxdw(j))/(-4*(sqrt(CFx1)-1));                                            %rovnice [47]
    else                                                    %Hybnostní teorie      
        udw_new(j)=0.5+0.5*sqrt(1-CFxdw(j));                                                        %rovnice [34]
    end
       
     chybadw(j)=abs(udw(j)-udw_new(j));                     %Kontrola konverengece výpočtu
     
        if chybadw(j) > presnost
            udw_new(j)=udw(j)+r_dw.*(udw_new(j)-udw(j));    %Použití relaxačního faktoru            %rovnice [63]
        end
   counterdw=counterdw+1;
   disp(udw_new(j))
   end
    
   if V_e(j)<=0         %%%Případ záporné rychlosti v ekvilibiru - nulování%%%
       udw_new(j)=NaN;
       V_e(j)=0;
       Wdw(j)=R*omega;
       alpha_dw_deg(j)=alfa0;
       alpha_dw(j)=alfa0_rad;
       Vdw(j)=0;
       
       CLdw(j)=interp1(x,y,alpha_dw_deg(j));     
       CDdw(j)=interp1(x,yy,alpha_dw_deg(j));
       
%      [CLdw(j),CDdw(j)]=getCLCD(Re_dw(j),alpha_dw_deg(j),Avektor,Bvektor,Cvektor)
       
       CNdw(j)=+CLdw(j).*cos(alpha_dw(j))+CDdw(j).*sin(alpha_dw(j));
       CTdw(j)=CLdw(j).*sin(alpha_dw(j))-CDdw(j).*cos(alpha_dw(j));
       CFadw(j)=-CTdw(j).*sin(thetadw(j))+CNdw(j).*cos(thetadw(j));  
   end
   j=j+1;                                                   %Výpočet následující produové trubice    
  end
  
%%%Výpočty sil a točivého momentu%%%
%Rovnice [48-60]
  for l = 1:numel(theta)
      FN_up(l)=0.5.*CN(l).*ro.*c.*h.*Wup(l).*Wup(l);
      FT_up(l)=0.5.*CT(l).*ro.*c.*h.*Wup(l).*Wup(l);
  end

  for k = 1:numel(thetadw)
      FN_dw(k)=0.5.*CNdw(k).*ro.*c.*h.*Wdw(k).*Wdw(k);
      FT_dw(k)=0.5.*CTdw(k).*ro.*c.*h.*Wdw(k).*Wdw(k);
  end
 
FT_up_ave=(1/(2*pi))*trapz(theta, FT_up);
Torque_up=N*FT_up_ave*R;
FT_dw_ave=(1/(2*pi))*trapz(thetadw, FT_dw);
Torque_dw=N*FT_dw_ave*R;
PT_up=Torque_up*omega;
PT_dw=Torque_dw*omega;
P_wind=0.5*2*R*h*ro*V_out^3;
CP_up=PT_up/P_wind;
CP_dw=PT_dw/P_wind;
CP=(CP_up+CP_dw);

%%%
pocet_iteraci=counter+counterdw;
toc