clear all
clc

W=300; %kg massa
S=14.1; % m2 Área da asa
b=15; % m envergadura
rho=1.225; %kg/m3 densidade do ar no n´vel do mar
g=9.8; %m/s2 aceleração da gravidade

alpha = [-4 -3 -1 1 3 5 7 9 11 12];
CL = [0.12 0.25 0.49 0.7 0.9 1.08 1.23 1.36 1.46 1.47];
CD = [0.016 0.018 0.022 0.0275 0.035 0.044 0.0535 0.0675 0.0865 0.095];

%CL/CD CL3/CD2

for i=1:length(CL)
    CLCD(i)=CL(i)/CD(i);
    CL3CD2(i)=CL(i)^3/CD(i)^2;
    Vplaneio(i) = (2*W*g/(rho*S*CL(i)))^0.5;
    RD(i) = -(2*W*g/(rho*S)*1/CL3CD2(i))^0.5;
    gamma(i) = -1/(CL(i)/CD(i));
    Vhori(i) = RD(i)/tan(gamma(i));
    V(i) = RD(i)/sin(gamma(i));
end

[CL3CD2max, index1] = max(CL3CD2);
[CLCDmax, index2] = max(CLCD);

menor_ang_planeio_max_CLCD=gamma(index2)
RDminimo_max_CL3CD2_=RD(index1)
Vh_maior_alcance = RD(index2)/sin(gamma(index2))*cos(menor_ang_planeio_max_CLCD)
Vh_maior_autonomia = RD(index1)/sin(gamma(index1))*cos(gamma(index1))
Vstall=min(Vplaneio)

% Maximo tempo de voo - RD mínimo.

% Calculo da atmosfera:
T0=288.15; %K Temp. ISA @MSL
DISA=15; %Delta ISA
H=linspace(5000,0,10);

for i=1:length(H)
    Tstd(i)=288.15-0.0019812*H(i);
    delta(i)=(Tstd(i)/T0)^5.2561;
    teta(i)=(Tstd(i)+DISA)/T0;
    sigma(i)=delta(i)/teta(i);
    Hp(i)=H(i)*((Tstd(i)+DISA)/Tstd(i));
end


tempo_total=0;
for j=1:length(H)-1
    DH(j)=Hp(j)-Hp(j+1);
    RDmin(j)= (2*W*g/(rho*sigma(j)*S)*1/CL3CD2max)^0.5;
    Dt(j) = DH(j)*0.3048/RDmin(j);
    tempo_total=tempo_total+Dt(j);
end

tempo_total_min=tempo_total/60

Alcance = max(Hp)*0.3048*max(CLCD)/1000
%Gráficos 

figure(1)
subplot(3,1,1)
plot(CD,CL,'-b');

%legend('CL vs CD')
xlabel('C_D')
ylabel('C_L')
grid on
grid minor

subplot(3,1,2)
plot(CLCD,CL,'-r');
%legend('CL vs C_L/C_D')
xlabel('C_L/C_D')
ylabel('C_L')
grid on
grid minor

subplot(3,1,3)
plot(CL3CD2,CL,'-r');
%legend('CL vs C_L^3/C_D^2')
xlabel('C_L^3/C_D^2')
ylabel('C_L')
grid on
grid minor

figure(2)
subplot(3,1,1)
plot(alpha,Vplaneio,'-b');
%legend('CL vs CD')
xlabel('\alpha')
ylabel('Vel. de Planeio [m/s]')
grid on
grid minor

subplot(3,1,2)
plot(alpha,RD,'-r');
%legend('CL vs C_L/C_D')
xlabel('\alpha')
ylabel('Razão de descida [m/s]')
grid on
grid minor

subplot(3,1,3)
plot(alpha,gamma,'-r');
%legend('CL vs C_L^3/C_D^2')
xlabel('\alpha')
ylabel('\gamma')
grid on
grid minor

figure(3)
subplot(2,1,1)
plot(V,RD,'-r');
%legend('CL vs C_L/C_D')
xlabel('Velocidade da aeronave [m/s]')
ylabel('Razão de descida [m/s]')
grid on
grid minor

hold on
plot([0 max(V)],[0 max(V)*sin(menor_ang_planeio_max_CLCD)],'--k');

subplot(2,1,2)
plot(Vhori,RD,'-b');
%legend('CL vs C_L/C_D')
xlabel('Velocidade horizontal [m/s]')
ylabel('Razão de descida [m/s]')
grid on
grid minor