% Cálculo de Potência: Dener A. L. Brandão - UFMG 2019
R = 37.5; % Raio das pás: 37,5 m
P_turb = 1.5e6; % Potência da turbina: 1,5 MW
p = 1.225; % Densidade do ar (1,225 kg/m^3 a 15ºC e 1 atm
A = pi*R^2;
beta = 0; % Ângulo de ataque igual a 0
i = 0;
j = 0;
for V_ven = 5.5:1:14.5 % Variação da velocidade do vento;
    j = j+1;
    fator_rad = pi/30;
    w_turb = 0:3:30; % Variação de w_turb de 0.1 a 30 rpm
    w_turb1 = fator_rad.*w_turb;
    lambda = w_turb1.*R/V_ven;
    lambda_i = 1./(1./(lambda + 0.08*beta) - 0.035/(beta^3 + 1)); % Fórmula de lambda_i
    Cp = 0.22.*(116./lambda_i - 0.4*beta - 5).*exp(-12.5./lambda_i); %Fórmula de Cp.
    P = p*A*V_ven^3.*Cp/2;
    handle = plot(w_turb, P/1000); %Plota cada vetor de potência
    ylabel('Potência (kW)','FontSize',15);
    xlabel('V_{Turbina} (RPM)','FontSize',15);
    title('Gráfico Potência x Velocidade da Turbina','FontSize',15);
    set(gca,'FontSize',15); % Define tamanho da fonte dos eixos
    set(handle,'LineWidth',[2]);
    grid on;
    hold on;
    legend on
    legend('5,5 m/s', '6,5 m/s','7,5 m/s','8,5 m/s','9,5 m/s','10,5 m/s','11,5 m/s','12,5 m/s','13,5 m/s','14,5 m/s','Max_P');
    axis([0 30 0 4e3]);
end
for i = 1:1:length(P)
        Max_P(i) = max(P(i));
        if (Max_P(i) >= 1.5e6)
            Max_P(i) = 1.5e6;
        end
    end
handle2 = plot(w_turb, Max_P/1000,'m');
set(handle2,'LineWidth',[2]);
hold on