%% parametri generali

M=0.22% mach
h =2000 %quota desiderata [m]
W=600*9.81; %[N]
Pa_sl= 90; %potenza motore see level [hp]

%ala
S=9.51; %[m^2] superficie alare
S1=102.36; %[ft2] superficie alare
DHDADI= 4; %[deg] angolo diedrio
CLalpha_h=5.73; %polare profilo
A=6.3; 
tip_ratio=0.34; % lambda
b= 12.7265*2; % [ft]
lambda_quartiw=6.5; % media angoli freccia ala misurati

%tail orrizzontale
Sh= 4.282; %[m2] superficie orrizontal tail
Sh1= 45.43; %[ft2]


%tail verticale
Sv=2.14; %[m2]
lambda_mezziv=50;%[deg]
CLalpha_v=6.33;% del profilo del tail verticale 
lamda_quartiv=55; %[deg]
zv= 2.25 %[ft] distanza asse x e fuoco tail verticale misurata

% Valutazione della temperatura e della densita' al variare della quota (fino a tropopausa)

rho_sl  = 1.225;               % Densità aria a livello del mare [kg/m^3]
T_sl    = 288.15;                  % Temperatura a livello del mare [K]
a       = -0.0065;             % Gradiente termico [K/m]
T  = T_sl + a*h;
rho = rho_sl*(T/T_sl)^4.25;

% Cy_bdot p 2076  p 1938

% eta_h pag 1961 
U1=M*sqrt(1.4*287*T); %[m/s] vorrà forse averlo in ft/s?
Pav=(rho/rho_sl)*Pa_sl; %[hp]
q=0.5*rho*U1^2;% [Pa]
Dp=1.67; %[m] diametro propulsore
Dp1= 5.84; %[ft]
Sh_slip= 2.862;% [m2]pagina 1960 baganta dal propulsore
Sh_slip1= 30.8;% [ft2]

eta_h=1+(Sh_slip/Sh)*((2200*Pav)/(q*U1*pi*Dp1^2)) %(p269 o 271 del libro, non del file) da ricontrollare, viene più di 5
eta_V=1+(1.98/Sv)*((2200*Pav)/(q*U1*pi*Dp1^2)) 

%%%Vh%%%
amc= 1252; %[mm]
xach= 4081.5/amc; % posizione fuoco adim.
xcg= 280.35/amc;% posizione xg adim. ala rispetto bordo attacco

Vh= (xach-xcg)*(Sh/S)

%clalphav p 2076  p 1938
lv=12.14; % [ft]distanza ta xcg aereo e ac tail verticale (10.27 p 2086)
%  Avf/Av= 1.3 Av=(1.18^2/2.14) kvh xv=0.65 xv/bv=0.55 x/xcr=0.3
% avhf/avf=0.85 kvh=1.1
Aveff= 1.58

%kv grafici
kv=0.9

%CLalphav
k=CLalpha_v/((1-M^2)^0.5)/(2*pi) %CLalpha_v è coeff di portanza del vt
beta=(1-M^2)^0.5
CLalphav=2*pi*Aveff/(2+ (Aveff^2*beta^2/k^2*( 1+(tand(lambda_mezziv))^2 )+4)^0.5) %(p 1938)


% 1+ de sigma/ de beta per eta zw/zf=1
zf=1;% diametro secondo napolitano
zw=0.35;% [m] distanza asse x e punta ala in fusoliera 
uno_deltasigmabeta_eta= 0.724+ 3.06*((Sv/S)/(1+cosd(lamda_quartiv))) + 0.4*zw/zf+0.009*A



cy_betav=-kv*CLalpha_v*uno_deltasigmabeta_eta*Sv/S
cn_beta=-cy_betav*lv/b




%de sigma/de beta 
% zv/bv=0.58
sigmabetagamma=-0.61;
sigmabaetat=-0.005;
st=3*pi/180;
sigmabetawingf=0.045;
sigma_beta= sigmabetagamma*(DHDADI/57.3)-sigmabaetat*st+sigmabetawingf
lp= 12.74; %[ft] distanza ra i qaari ala e i quarti tail verticale 

Cybdot=2*CLalphav*sigma_beta*Sv/S*lp/b

% Calcolo Cyr (p 2118)  
lv=12.14; % [ft]distanza ta xcg aereo e ac tail verticale (10.27 p 2086)
Cybetav =-kv*CLalphav*uno_deltasigmabeta_eta*(Sv/S); %ricavabile da 10.28 (p 2076)
Cyr = -2*(Cybetav)*(lv/b) %equazione valida per alpha nullo



% Cmde (p 2125)
kb = 0.8;%da p 259(del libro, non del file)che rimanda a fig 8.51 ma coi dati di 8.52 (p 1950 del file)
cldelta_cldeltath = 0.8;%da fig 8.15 (p 1920) necessario il clalpha_clalphath SONO CL
cldeltath = 4;%da fig 8.14 p 1918 è CL
kprimo = 0.8;%da fig 8.13 (ha veramente senso tenerlo in conto?)
alphadeltaCl_alphadeltacl = 1.2;% da fig 8.53
%NOTA: tutti i valori presi da figura appena annotati sono presi un po' a caso e da verificare 
Clalphah = CLalpha_v; %assunto uguale a quello del vtail perchè credo abbiano lo stesso NACA. Credo.
alphade = kb*(cldelta_cldeltath)*cldeltath*(kprimo/Clalphah)*(alphadeltaCl_alphadeltacl); %(da 10.94 p 2127)

Cmih = -Clalphah*eta_h*Vh;%da eq 10.92 (p 2126)
Cmde = alphade*Cmih%(p 2128) il valore finale dovrebbe aver senso

% Cyda
Cyda = 0; %per la maggior parte dei velivoli

% Clda: procedura descritta a p 2132 e ss
SPANFI= 7.7; % [ft]
SPANFO= 12.217; % [ft]
 
%etai= SPANFI/(b/2); etao= SPANFO/ (b/2); k_beta=(CLalpha_h/((1-M^2)^0.5))*beta/(2*pi);  % 10.54 pag 2109 
% beta*A/k_beta; lambda_beta= atand(tand(lambda_quartiw)/beta); valori
% utili per trobvare beta_cldelta/k =0.35
k_beta=(CLalpha_h/((1-M^2)^0.5))*beta/(2*pi) ;

beta_cldelta_k =0.35; % dal grafico

Cl_delta_primo= (k_beta/beta)*beta_cldelta_k;
CL_alphaa= 5.73; % polare alettone che per ora è ugyuale a quello dell'ala

cl_delta1= 0.85*4; % 0.85 relatico cldelta teorico trovato da grafico, 4 è cldelta teorico da grafico pag 1918
alpha_deltaa=cl_delta1/CL_alphaa;

Cl_delta_a=(alpha_deltaa)*Cl_delta_primo

% Cn delta a %%%
ka=-0.16;
Clw=W/(0.5*rho*U1^2*S); %dipede dalla velocità

Cn_delta_a=Cl_delta_a*ka*Clw
%% Cy delta r pag 2151%%%

kprimo_r=0.8; % pag
kb_r= 0.91;%figura 8.52 pag 1950
 cfr_c=0.33;
cldeltar_cldeltarteorico=0.86; %pag 1920
cldeltar_teorico=4.8;%pag 1920

Cy_deltar= CLalphav*kprimo_r*kb_r*cldeltar_cldeltarteorico*cldeltar_teorico*Sv/S

%% Cl delta r pagina 2151%%%

zv=2.25;%[ft] distanza xg aereo e ac tail verticale

Cl_deltar= zv*cosd(0)/b*Cy_deltar

%% Cn delta r pagina 2152%%%

Cn_deltar=-Cy_deltar*lv/b

%% Cm alpha pag 2071%%

%pag 1962 CL_alpha_roskam
df=1.33*2;% diametro sul piano più largo
Kwf= 1+0.025*(df/b)-0.25*(df/b)^2
k_w=CLalpha_h/((1-M^2)^0.5)/(2*pi); %CLalpha_h è della wing
beta=(1-M^2)^0.5;
A1=b^2/S1;
lambda_mezziw=2;%deg
CLalphaW=2*pi*A1/(2+ (A1^2*beta^2/k_w^2*( 1+(tand(lambda_mezziw))^2/beta^2 )+4)^0.5) %(p 1938)
Cl_alphawf=CLalphaW*Kwf

KA=(1/A1)-1/(1+A^(1.7));
Klamda=(10-3*tip_ratio)/7;
Hh=1.32; %[ft]
lh=12.29;%[ft] distanza xaw e xatail orrizontale
Kh=(1-Hh/b)/(2*lh/b)^(1/3)
da_dalpha=4.44*(KA*Klamda*Kh*(cosd(lambda_quartiw)^0.5))^1.19*CLalphaW/CLalpha_h;
CL_alpharoskam=Cl_alphawf+CLalpha_v*eta_h*(Sh/S)*(1-da_dalpha)
cm_alpha=(-0.05)*CL_alpharoskam





