%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%              Suddivisione masse Raymer 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%suddivisone masse con coefficienti presi dal raymer a pag 599. [kg]
%NB peso e considerazizone all'elica
wtotale=600; %[kg] PESO A VUOTO PRESO DA MANUALE 
wala= 0.16*wtotale;
wtail_o=0.043*wtotale;
wtail_v=0.033*wtotale;
motore=60;
propeller=5;
carrello=0.055*wtotale;

pilota1=90;
pilota2=90;
pilota=pilota1+pilota2;
bagalio= 27;
fuel=43; %nb i tank contengono al max 31l ciascuno


body=wtotale-(wala+wtail_o+ wtail_v+ motore+...
+carrello+propeller+ pilota+ bagalio+ fuel);

W_vuoto=wtotale-(pilota+ bagalio+ fuel);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%              Suddivisione masse Raymer 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%suddivisone masse con coefficienti presi dal raymer a pag 599. [kg]
%NB peso e considerazizone all'elica
wtotale=600; %[kg] PESO A VUOTO PRESO DA MANUALE 
wala= 0.16*wtotale;
wtail_o=0.043*wtotale;
wtail_v=0.033*wtotale;
motore=60;
propeller=5;
carrello=0.055*wtotale;

pilota1=90;
pilota2=90;
pilota=pilota1+pilota2;
bagalio= 27;
fuel=43; %nb i tank contengono al max 31l ciascuno


body=wtotale-(wala+wtail_o+ wtail_v+ motore+ carrello ...
    +propeller+ pilota+ bagalio+ fuel);

W_vuoto=wtotale-(pilota+ bagalio+ fuel);

%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%               suddivisione masse con BS115
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%percentuali rispetto basic empty w.
%ala 12.4 percento, tail verticale 3.37 perc, tail H 3.37 perc, body 39,88
%perc, motore 27.97 perc, propeller 4.63 perc, nouse landig gera 2.4 perc,
%min landig gera 6.06 perc

wempty = 390; %[kg] PESO A VUOTO PRESO DA MANUALE 
wala1 = 55;%0.124*wempty 
wtail_o1 = 15;%0.0337*wempty 
wtail_v1 = 13;%0.0337*wempty
%motore1=0.2798*wempty
motore1 = 67.7;
propeller1 = 18;%0.0463*wempty
carrello_main1 = 24;%0.0606*wempty
carrello_nose1 = 10;%0.0240*wempty
body1 = 31.7680+155.5320;%0.3988*wempty

somma = wala1+wtail_o1+wtail_v1+motore1+propeller1+ ...
    carrello_main1+carrello_nose1+body1;

%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                       Calcolo baricento 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%% PIANO DI CARICO %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
pilota1 = 85;
pilota2 = 85;
pilota=pilota1+pilota2;
bagaglio = 0;
fuel = 40;

densita_fuel = 737.4;% [kg/m3]valor medio ma da cambiare in base al tipo
volume_massimo_serbatoi = 62;%[litri]
massa_massima_fuel = densita_fuel*volume_massimo_serbatoi/1000;
%nb i tank contengono al max 31 litri ciascuno
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%% Massa elementi %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wala = wala1;
wtail_o = wtail_o1;
wtail_v = wtail_v1;
motore = motore1;
propeller = propeller1;
carrello_main = carrello_main1;
carrello_nose = carrello_nose1;
body = body1;

%massimo peso decollo
MTOW = wempty+pilota1+pilota2+bagaglio+fuel;


%%%%%%%%%%%%%%%%%%%%%%%  Posizioni lungo x y z %%%%%%%%%%%%%%%%%%%%%%%%


%posizioni lungo x a partire dalla punta del muso [m]
x_propeller = 0;
x_motore = 0.91;
x_tailo = 6.35;
x_carrello_nose = 0.950;%da manuale di volo
x_carrello_main = 2.688;%da manuale di volo
x_fuel = 2.2; %da manuale di volo
x_tailv = 6.35;

%x baricentro fusoliera approssimata come due coni 
x_troncomuso = 2.9116*3/4;
x_troncocoda = ((6.7936-2.9116)/4)+2.9116;
x_body_troncomuso = (x_troncomuso*body*2/3+x_troncocoda*body*1/3)/(body);

%baricentro fusoliera usando similitudine BS115
x_bodi115=(75.5*1003+3.2*(-1001)+1.3*(-909)+33.3*2122+5*(-1300)+ ...
    4.2*(-1350)+12.7*(-260)+7.5*(-450)+1.3*(150)+13.5*(534)+ ...
    10.7*(-1046)+4.7*(365)+18.5*(273)+11*(45)+2.5*(507)+ ...
    2.5*(1282))/207+800; 

x_bodyBSP=(x_bodi115*7437/7178)/1000;%rispetto al datum

%posizioni lungo x a partire dal datum [m]
x_datum = 1.45;
x_propeller1 = 0-x_datum;
x_motore1 = 0.91-x_datum;
x_ala1 = 0.68+0.35*1.252; %pos del terzo della m.a.c. rispetto al datum
x_tailo1 = 6.35-x_datum;
x_tailv1 = 6.25-x_datum;%misurate 
x_pilota11 = 1; % manuale
x_pilota21 = 1.8;% manuale 
x_carrello_main1 = x_carrello_main-x_datum;
x_carrello_nose1 = x_carrello_nose-x_datum;
x_fuel1 = 2.2-x_datum; %manuale
x_body1 = x_bodyBSP;
x_bagaglio1 = 2.25;%manuale di volo



%posizione lungo y [m]

%wing
hw = 7.76/2;
br = 1.9;% corda al root
bt = 0.67;%corda al tip
bp = 1.3720;% corda al break point
he = 2.7029;% distanza outboard panel

y_ala1e = -he/3*(bp+2*bt)/(bp+bt);% semiala sx negativa outboard panel
y_ala1i = -(hw-he)/3*(br+2*bp)/(br+bp);% semiala sx negativa inboard panel

Si = (br+bp)*(hw-he)/2;%superficie interna trapezio ala
Se = (bt+bp)*(he)/2;%superficie esterna trapezio ala

y_semiala1 = (y_ala1e*Se+y_ala1i*Si)/(Si+Se); %semiala sisnistra negativo
y_semiala2 = -y_semiala1;%semiala destra positivo

%tank fuel 
y_fuel1 = -0.732;% tank  sinistro negativa
y_fuel2 = 0.732;% tank  destro positiva

%tail
h = 1.51;
B = 0.97;
b = 0.43;
y_tailo1 = -h/3*(B+2*b)/(B+b);%sinistra negativa
y_tailo2 = h/3*(B+2*b)/(B+b);%destra positiva
y_bagaglio = 0;
y_motore = 0;
y_pilota1 = 0;
y_pilota2 = 0;
y_carrello = 0;
y_body = 0;
y_tailv = 0;
y_propeller = 0;

%posizione lungo z tutto in metri concorde datum verso l'alto 
z_ala = -0.35; %posizione della radice dell'ala, misurata dal disegno
z_semiala = -0.35+y_semiala2*tand(4); %tenendo conto di angolo diedro
z_fuel = -0.30;
z_motore = -0.10; 
z_pilota1 = 0;
z_pilota2 = 0.1; 
z_tailo = 0.113-y_tailo2*tand(3.5);
z_tailv = 0.7; 
z_carrello_nose = -0.30; 
z_carrello_main = -0.28; 
z_body = -0.1; 
z_propeller = 0;
z_bagaglio = -0.1;

%%%%%%%%%%%%%%%%%%%%%%%%% CALCOLO BARICENTRO %%%%%%%%%%%%%%%%%%%%%%%%%%%%

% calcolo baricentro sfruttando manuale di volo
fuel_momentum = fuel*(30/40);%fuel in kg
pilota1_momentum = pilota1*1;
pilota2_momentum = pilota2*(125/70);
aircraft_momemntum = wempty*(350/390);
bagaglio_momentum = bagaglio*(33.7/15);

somma_momentum=(fuel_momentum+pilota1_momentum +pilota2_momentum+ aircraft_momemntum+bagaglio_momentum);

% BARICENTRO LUnGO X calcolo usando la formula del manuale
x_g_percentuale_mac=(((pilota1_momentum+pilota2_momentum+ ...
    fuel_momentum+bagaglio_momentum+ ...
    aircraft_momemntum)/MTOW)-0.68)/1.252*100;

x_g_mac_datum= (x_g_percentuale_mac*1252/100+680)/1000;  
%[m] posizione di xg a partire da x_g_percentuale_mac a cui si aggiunge 
%la distanza dal datum della corda media a.

xg = x_g_mac_datum;



% Baricentro Zg rispetto al datum
zg = (z_motore*motore+z_semiala*wala+z_tailo*wtail_o+ ...
    z_tailv*wtail_v +z_pilota1*pilota1+z_pilota2*pilota2+ ...
    z_carrello_main*carrello_main+z_carrello_nose*carrello_nose+ ...
    z_fuel*fuel+z_body*body+z_propeller*propeller+ ...
    z_bagaglio*bagaglio)/MTOW;

%posizioni rispetto al baricentro

x_motoreg = xg-x_motore1;
x_alag = xg-x_ala1 ;
x_tailog = xg-x_tailo1;
x_tailvg = xg-x_tailv1 ;
x_pilota1g = xg-x_pilota11;
x_pilota2g = xg-x_pilota21;
x_carrello_maing = xg-x_carrello_main1;
x_carrello_noseg = xg-x_carrello_nose1;
x_fuelg = xg-x_fuel1;
x_bodyg = xg-x_body1;
x_propellerg = xg-x_propeller1;
x_bagagliog = xg-x_bagaglio1;

z_semialag = zg-z_semiala;
z_fuelg = zg-z_fuel ;
z_motoreg = zg-z_motore; 
z_pilota1g = zg-z_pilota1;
z_pilota2g = zg-z_pilota2 ;
z_tailog = zg-z_tailo;
z_tailvg = zg-z_tailv;
z_carrello_maing = zg-z_carrello_main ;
z_carrello_noseg = zg-z_carrello_nose;
z_bodyg = zg-z_body;
z_propellerg = zg-z_propeller;
z_bagagliog = zg-z_bagaglio;

%%%%%%%%%%%%%%%%%%%%%%%%% Momenti di inerzia %%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%Momenti di inerzia BARICENTRICI

Ixx = (motore*(z_motoreg^2+y_motore^2) + ...
    wala*(z_semialag^2+y_semiala2^2) + ...
    fuel*(z_fuelg^2+y_fuel2^2) + ...
    pilota1*(z_pilota1g^2+y_pilota1^2) + ...
    pilota2*(z_pilota2g^2+y_pilota2^2) + ...
    wtail_o*(y_tailo2^2+z_tailog^2) +...
    wtail_v*(z_tailvg^2+y_tailv^2) + ...
    carrello_nose*(y_carrello^2+z_carrello_noseg^2)   +   ...
    carrello_main*(y_carrello^2+z_carrello_maing^2)   +  ...
    body*(z_bodyg^2+y_body^2) + ...
    propeller*(y_propeller^2+z_propellerg^2) + ...
    bagaglio*(z_bagagliog^2+y_bagaglio^2));               

Iyy = (motore*(z_motoreg^2+x_motoreg^2) + ...
    wala*(z_semialag^2+x_alag^2)     + ...
    fuel*(z_fuelg^2+x_fuelg^2) + ...
    pilota1*(z_pilota1g^2+x_pilota1g^2) + ...
    pilota2*(z_pilota2g^2+x_pilota2g^2) + ...
    wtail_o*(x_tailog^2+z_tailog^2 )+ ...
    wtail_v*(x_tailvg^2+z_tailvg^2) + ...
    carrello_nose*(x_carrello_noseg^2+z_carrello_noseg^2) + ...
    carrello_main*(x_carrello_maing^2+z_carrello_maing^2) + ...
    body*(z_bodyg^2+x_bodyg^2) + ...
    propeller*(x_propellerg^2+z_propellerg^2) + ...
    bagaglio*(z_bagagliog^2+x_bagagliog^2));               

Izz = (motore*(y_motore^2+x_motoreg^2) + ...
    wala*(y_semiala1^2+x_alag^2)     + ...
    fuel*(y_fuel1^2+x_fuelg^2) + ...
    pilota1*(y_pilota1^2+x_pilota1g^2)  + ...
    pilota2*(y_pilota2^2+x_pilota2g^2) + ...
    wtail_o*(x_tailog^2+y_tailo1^2) + ...
    wtail_v*(x_tailvg^2+y_tailv^2) + ...
    carrello_nose*(x_carrello_noseg^2+y_carrello^2)   +    ...
    carrello_main*(x_carrello_maing^2+y_carrello^2)   +  ...
    body*(y_body^2+x_bodyg^2) + ...
    propeller*(x_propellerg^2+y_propeller^2) + ...
    bagaglio*(x_bagagliog^2+y_bagaglio^2));               

Jxz = -(motore*z_motoreg*x_motoreg     + ...
    wala*z_semialag*x_alag      +    ...
    fuel*z_fuelg*x_fuelg    +   ...
    pilota1*z_pilota1g*x_pilota1g     +  ...
    pilota2*z_pilota2g*x_pilota2g    +  ...
    wtail_o*x_tailog*z_tailog   +   ...
    wtail_v*x_tailvg*z_tailvg    +  ...
    carrello_nose*x_carrello_noseg*z_carrello_noseg   +  ...
    carrello_main*x_carrello_maing*z_carrello_maing   +  ...
    body*z_bodyg*x_bodyg    + ...
    propeller*x_propellerg*z_propellerg    + ...
    bagaglio*x_bagagliog*z_bagagliog);               

%inerzie aggiuntive rispetto ai singoli baricentri dove ali e tail visti
%come aste e cilindo lungo y, body come cilindro, motore cilindro. 

%ali
Ixx_wing=(1/12)*wala*(hw)^2;% vista come asta
Izz_wing=Ixx_wing;
Iyy_wing=(1/2)*wala*(0.626)^2; %visto come cilindro 

%tail orizzontale 
Ixx_tailo=(1/12)*wtail_o*(1.5091)^2;
Izz_tailo=Ixx_tailo;
Iyy_tailo=(1/2)*wtail_o*(0.35)^2; % media tra c tip e c root

%tail verticale
Ixx_tailv=(1/12)*wtail_o*(0.9787)^2; %visto come asta
Iyy_tailv=Ixx_tailv;
Izz_tailv=(1/2)*wtail_o*(0.55)^2;% visto come cilindro 

%body visto come due coni
Ixx_body=(3/10)*body*2/3*(0.5)^2+(3/10)*body*1/3*(0.5)^2;
Iyy_body=(1/12)*body*(7.178)^2;
Izz_body=Iyy_body;

Ixx_motore=(1/2)*motore*0.3506^2;
Iyy_motore=(1/12)*motore*(0.7112)^2;
Izz_motore=Iyy_motore;

Ixx_tot=Ixx+Ixx_wing+Ixx_tailv+Ixx_tailo+Ixx_body+Ixx_motore;
Iyy_tot=Iyy+Iyy_wing+Iyy_tailv+Iyy_tailo+Iyy_body+Iyy_motore;
Izz_tot=Izz+Izz_wing+Izz_tailv+Izz_tailo+Izz_body+Izz_motore;


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

% controlli 
if(MTOW>600) fprintf ('stop MTOW>600 kg\n')
end

if(MTOW<=600) && (massa_massima_fuel>fuel)
    fprintf ('MTOW ok & FUEL ok\n')
    fprintf('\n coordinate del baricentro rispetto al datum')
    Coord_x_e_z = [xg zg]
     fprintf('xg come percentuale della corda media aerodinamica\n')
    x_g_percentuale_mac
   
    if(x_g_percentuale_mac>= 21.5 && x_g_percentuale_mac<=36)
       fprintf('21.5< xg  MAC <36 soddisfa limiti\n ')
   end
   if(x_g_percentuale_mac< 21.5)
       fprintf('xg  MAC <21.5 non soddisfa limiti \n')
   end
   if(x_g_percentuale_mac> 36)
       fprintf('xg MAC >36 non soddisfa limiti \n')
   end
if (somma_momentum>(MTOW-470)*0.866+420 && somma_momentum<(MTOW-350)*0.890+400)
fprintf(' \n Peso e momentum ok nel Weight and Moment Envelope\n')

else
fprintf('\n  Peso e momentum non ok nel Weight and Moment Envelope\n')
end
end
    fprintf('\n momenti di inerzia Ixx, Iyy, Izz, Jxz')
    Inerzie = [Ixx Iyy Izz Jxz]
   

if(MTOW<460) 
    fprintf ('stop MTOW to light\n ')
end


if(massa_massima_fuel<fuel) 
    fprintf ('stop to much fuel, fuell> max fuel\n ')
end


