%%% Momento beccheggio velivolo completo %%%

%% Velivolo Jetstar %%

clear all
clc

%% Dati %%

% ala
S = 542.5 %[ft^2]
b = 53.75 %[ft]
mgc = 10.93 %[ft]
sweptc4 = 30 %[deg]
iw = 2 %[deg]
e = 0.9
AR = b^2/S

% sezioni alari
cmac = -0.04
claw = 0.105 %[/deg]
a0L = -1 %[deg]
Xacw = 0.25

% piano orizzontale di coda
Sh = 149 %[ft^2]
bh = 24.75 %[ft]
mgch = 6.5 %[ft]
lh = 23.6 %[ft]
sweptc4h = 30 %[deg]
ih = -1 %[deg]
ARh = bh^2/Sh

% sezioni piano di coda
clah = 0.10 %[/deg]
cmde = -0.81 %[/rad]

% fusoliera
lf = 58.6 %[ft]
dmax = 7.2 %[ft]
Xcg = 0.25

% altri dati velivolo (M = 0.20 sea level)
cL = 0.737
cma = -0.80 %[/rad]
K2K1 = 0.92

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

cLaw = claw*57.3/(1 + claw*57.3/(pi*AR*e)) % variazione del cL con alpha per l'ala [/rad]
cLah = clah*57.3/(1 + clah*57.3/(pi*ARh*e)) %... per piano di coda [/rad]
e0 = 2*cL/(pi*AR*e) % epsilon 0 [rad]
depsw = 2*cLaw/(pi*AR*e) % variazione dell'angolo di downwash con alpha sull'ala [/rad]
cL0wf = cLaw/57.3*(iw-a0L) % coeff. di portanza a zero alpha
Vh = Sh/S*lh/mgc % rapporto volumetrico del piano di coda

% discretizzazione fusoliera (metodo delle strisce)
DX = [2.6 2.6 2.2 2.2 3.4 3.4 3.4 3.4 14.6 5.2 5.2 5.2 5.2]; %[ft]
wf = [2.5 4.2 5.5 6.3 6.6 7.2 7.2 7.2 7.2 6.6 5.4 3.8 2.1]; %[ft]
ifs = [-3 -3 -10 -10 0 0 0 0 0 0 0 0 -4]; %[deg]

% contributo della fusoliera al coeff. di momento
cmof = zeros(13,1);
for j = 1:13
    cmof(j) = K2K1/(36.5*S*mgc)*wf(j)^2*(a0L+ifs(j))/57.3*DX(j);
end
cm0f = sum(cmof)

% upwash e downwash
% distanze centroidi tronchi fusoliera dal bordo d'attacco delle ali
xbdaa = sum(DX(1:8)); % avanti il bordo d'attacco
xbdad = sum(DX)-sum(DX(1:8)); % dietro il bordo d'attaccp

% distanze dei tronchi dal bordo d'attacco alare
x = zeros(1,13);
for j = 1:8
    x(j) = xbdaa - (DX(j)/2 + sum(DX(1:j-1)));
end
for j = 10:13
    x(j) = xbdad + DX(j)/2 - sum(DX(j:end));
end

xratio = x/DX(9) % rapporto distanze tronchi/corda di radice

% stima di depsilon/dalpha per upwash
epsa = [1.11 1.13 1.15 1.17 1.2 1.3 2.86 3.9 0 0 0 0 0];

% stima di depsilon/dalpha per downwash
for j = 10:13
    epsa(j) = x(j)/lh*(1-depsw);
end
epsa % depsilon/dalpha)

plot(xratio(1:6),epsa(1:6),'k-*')
hold on
plot(xratio(7:8),epsa(7:8),'r-*')
hold on
plot(xratio(10:13),epsa(10:13),'g-*')
title('d\epsilon/d\alpha lungo la fusoliera')
xlabel('Distanza x/mgc dei centroidi dei tronchi di fusoliera dal bordo di attacco alare')
ylabel('d\epsilon/d\alpha')
legend('upwash on stations 1:6','strong upwash on stations 7,8',...
   'downwash on stations 10:13')
grid on

% contributo della fusoliera alla variazione del coeff. di momento con alpha
for j = 1:13
    cmAf(j) = 1/(36.5*S*mgc)*wf(j)^2*epsa(j)*DX(j);
end
cmaf = sum(cmAf)

figure
plot(1:8,cmof(1:8),'k+-')
hold on
plot(10:13,cmof(10:13),'k+-')
hold on
plot(1:8,cmAf(1:8),'r+-')
hold on
plot(10:13,cmAf(10:13),'r+-')
legend('cmof','','cm\alphaf','')
title('cmo e cm\aplha dovuti ai singoli tronchi di fusoliera')
xlabel('Suddivisione della fusoliera')
ylabel('cmo ; cm\alpha')
grid on

% effetto MUNK
Xacwf = Xacw - cmaf/cLaw

% contributi parziali al momento di beccheggio
cm0 = cmac + cm0f + cL0wf*(Xcg-Xacwf) + Vh*cLah*e0

cma = cLaw*(Xcg-Xacwf) - Vh*(1-depsw)*cLah

cmih = - Vh*cLah

cmde % è un dato

alpha = [0:4:20];
de = [-20:5:20];

for j = 1:6
    cmf(j) = cm0f + cmaf*alpha(j)/57.3;
    cmw(j) = cmac + cLaw*(Xcg-Xacwf)*alpha(j)/57.3;
    cmh(j) = - Vh*cLah*(ih/57.3 - e0 + (1-depsw)*alpha(j)/57.3);
    cm(j) = cmf(j) + cmw(j) + cmh(j);
    for k = 1:9
        cmCG(j,k) = cm0 + cma*alpha(j)/57.3 + cmih*ih/57.3 + cmde*de(k)/57.3;
    end
end

%% NOTA: riga 145 e riga 147 sono due modi di esprimere il cm!
%% infatti la curva cmCG figura 3 grafico 4 è identica alla curva del cmCG
%% di figura 4 con deltae = 0

figure
subplot(2,2,1)
plot(alpha,cmf,'b')
grid on
title ('Coefficiente di momento dovuto alla fusoliera')
xlabel ('\alpha [deg]')
ylabel ('cmf')

subplot(2,2,2)
plot(alpha,cmw,'r')
grid on
title ('Coefficiente di momento dovuto all''ala')
xlabel ('\alpha [deg]')
ylabel ('cmw')

subplot(2,2,3)
plot(alpha,cmh,'g')
grid on
title ('Coefficiente di momento dovuto al piano di coda')
xlabel ('\alpha [deg]')
ylabel ('cmh')

subplot(2,2,4)
plot(alpha,cmf,'b')
hold on
plot(alpha,cmw,'r')
hold on
plot(alpha,cmh,'g')
hold on
plot(alpha,cm,'k')
grid on
title ('Coefficiente di momento - Ala, fusoliera, piano di coda, totale')
xlabel ('\alpha [deg]')
ylabel ('cm')
legend ('cm fuselage', 'cm wing','cm tail','cm CG')

figure
plot(alpha,cmCG)
title ('Coefficiente di momento del velivolo totale al variare della deflessione dell''equilibratore')
xlabel ('\alpha [deg]')
ylabel ('cm CG')
legend ('\deltae = -20°','\deltae = -15°','\deltae = -10°','\deltae = -5°',...
    '\deltae = 0°','\deltae = 5°','\deltae = 10°','\deltae = 15°','\deltae = 20°')
grid on

% posizione del punto neutro (approssimata)
Xn = Xacwf + Vh*cLah/cLaw*(1-depsw)