%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                          Magistral Theis:                               %               
%              Operating cost estimation methodology for eVTOL            %
%                        Prof. Davide Ferretto                            %
%                     Candidate Nadia Oteri 295925                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all
close all
clc

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

%% Input Data - Cost estimation

%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DIRECT COST %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%We                                             %[lbs] airplane empty weight
%W1                                             %[lbs] weight of wheels, brakes, tires and tu­bes
%W2                                             %[lbs] weight of engines
%W3                                             %[lbs] weight of starter
%W4                                             %[lbs] weight of cooling fluid
%W5                                             %[lbs] weight of rubber or nylon fuel cells 
%W6                                             %[lbs] weight of batteries, electrical power supply
%W7                                             %[lbs] weight of instruments, avionics and electronic equipment
%W8                                             %[lbs] weight of armament and fire- control system
%W9                                             %[lbs] weight of airconditioning units
%W10                                            %[lbs] weight of APU
%W11                                            %[lbs] weight of trapped fuel and oil
Wto = 3089 * 2.2;                   %[lbs] takeoff weight of the air­plane
%V_max                                          %[kts] V max in keas

% N_rtde                                        %[-] number of airplanes produced for the RDTE phase - commercial programs use: NRDTE = 2 - 8
% F_diff                                        %[-] judgement factor for the difficulty
% F_cad                                         %[-] judgement factor for the effect of computer aided
% Re_r                                          %[USD/mhr] engineering dollar rate per hour
% CEF                                           %[-] Cost Escalation Factor
Rbl = 93 / 1.852;                   %block distance measured in nm -> range tot + loiter 
t_cl = ((85.70+263.79)/60)/60;      %[hrs] time required to climb
Vcl = ((38.33+43.33)/2) * 1.944;    %projected horizontal speed of the airplane during the climb, in kts
Vde = 43.33 * 1.944;                %projected horizontal speed of the airplane during the descent, in kts
t_de = (263.79/60)/60;              %[hrs] time required to descent
t_cr = (953.50/60)/60;              %[hrs] time required to cruise
Vc = 66.67 * 1.944;                 % design cruise speed in kts
Vman = 250;                         %below 10000ft 250ks, maneuvering due ATC constraints

Nyr = 5;                            %number of years during which the airplane is operated

n_captain = 1;
n_copilot = 0;
n_flightengineer = 0;

Kj = 0.26;                          %vacation pay, cost of training, crew premium, crew insurance, and payroll tax
CEF = 2.53;                         %[-] Cost Escalation Factor for 2024 (pag 47 Carugo 4)
SAL_1989 = 75000;
%CEF_1989 = 3.01;                   %Cost Escalation Factor 1989
AH_j = 1500;                        %number of flight hours per year for a crew member
%TEF_1990 = 7;                      %(?) Travel expense factor
%CEF_1990 = 3.05;

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

%PARTE FUEL DA CAMBIARE pag 111
% Wf_used
% Wf_bl = Wf_used;                                   %block fuel used in lbs - 0.5per­cent increase in fuel consumption per year
% FP                                                 %price of fuel in USD/US gallon
% FD                                                 %fuel density in lbs/gallon
% Wol_bl = Wf_used/70;                               %weight of oil and lubricants used, in lbs - VEDI PAG 114
% OLP                                                %price of oil and lubricants in USD/gallon
% OD                                                 %oil density in lbs/gallon
% Cpol = (Wf_bl/Rbl) * (FP/FD) + (Wol_bl/Rbl) * (OLP/OD);
% %Cpol = 1.05 * (Wf_bl/Rbl) * (FP/FD)

% Energy used
Energy_used = 200; % example in kWh
PE = 0.15; % energy price in USD/kWh (2024)

% operative cost per energy
C_en = Energy_used * PE / Rbl;

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

f_ins_hull = 0.015;                               %0.005 to 0.030 USD/USD/airplane/year - annual hull insurance rate in USD per USD airplane price
AMP = 2250000;                                    %airplane market price in USD (2024)

Rl_ap = 18.40*CEF;                                      %airplane maintenance labor rate per manhour in USD/hr - 15-20
MHR_meng_bl = 0.03;                               %number of engine maintenance hours needed per block hour per engine
Rl_eng = 5*CEF;                                      %engine maintenance labor rate per man hour in USD/hr - 5-15
%Cmat_apblhr = 183.72;                            %airframe and systems maintenan­ce materials cost per airplane block
Cmat_apblhr = 7*CEF;                              %airframe and systems maintenan­ce materials cost per airplane block (2024)

MHR_map_flt = 0.04;                               %maintenance man hours per flight hour
%MHR_meng_flt = 0.07; 
C_mat_engblhr = 10*CEF;                               %engine maintenance materials cost per engine per airplane block hour in USD/hr (2024)
famb_lab = 1;                                     %overhead distribution factors for labor cost PAG 128
Ne = 8;                                           %(flie excel input data)
famb_mat = 0.35;                                  %overhead distribution factors for material cost PAG 128
MHR_eng_bl = 0.03; 
%MHR_mat_engblhr = 0.07; 
 
Fdap = 0.25;                                     %airframe depreciation factor
AEP = 2000000;                                   %airplane estimated price (pag 63 carugo) (2024)
EP = 150000;                                     %engine price (per engine) in USD - Battery Included (2024)
PP = 9000;                                       %price per propeller in USD 5,000-10,000 (2024)
Np = 8;                                          %(flie excel input data)
ASP = 150000;                                    %avionics systems price per airplane in USD (2024)
DP_ap = 5;                                       %airplane depreciation period 10ys
Fdeng = 0.25;                                    %engine depreciation factor
DP_eng = 3;                                      %engine depreciation period
Fdav = 0.15;                                     %avionics systems depreciation fac­tor
DP_av = 5;                                       %avionics systems depreciation period
Fdprp = 0.03;                                    %propeller depreciation factor - 2-5%
DP_prp = 3;                                      %propeller depreciation period
Fdapsp = 0.25;                                   %airplane spare parts depreciation factor
Fapsp = 0.10;                                    %airplane spare parts factor
DP_apsp = 2;                                     %airplane spare parts depreciation period
Fdengsp = 0.25;                                  %engine spare parts depreciation factor
Fengsp = 0.25;                                   %engine spare parts factor
ESPPF = 1.50;                                    %engine spare parts price factor
DP_engsp = 2;                                    %depreciation period for engine spare parts

Capnf = 0;                                        %navigation feecharged per airplane per flight in USD/flight
%Crt = f_rt * DOC;

Nacq = 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%% INDIRECT COST %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Cmls                            %cost for meal service
% Cpax_ins                        %cost for passenger insurance
% Ccat                            %cost for cabin attendants
% Cpax_gen                        %cost associated with the follow­ ing items: Passenger handling; 
%                                 %Passenger baggage handling; Sales and reservations; Security; Miscellaneous passenger costs
% 
% Cmgef                           %cost of maintaining ground equipment and facilities
% Cdgef                           %cost of depreciation of ground equipment and facilities
% 
% Caps                            %cost of airplane service
% Capctrl                         %cost of airplane control
% Cfrt                            %cost associated with handling freight
% IOCpse                          %indirect operating cost for promo­tion, sales and entertainment
%                                 %Commissions to travel agencies; Publicity and advertising campaigns; Entertainment
% 
% IOCgaa                          %indirect operating cost for general administrative expenses

%Programoperating_cost (??)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Operating cost

t_gm = 0.51*(10)^(-6)*(Wto) + 0.125;                                            %[hrs]time spent in ground maneuvers
Rcl=Vcl*t_cl;                                                                   %distance to climb and to accelerate to cruise speed, in nm
Rde=Vde*t_de;                                                                   %distance covered during the des­ cent, in nm
tman = 0.25*(10)^(-6)*Wto+0.0625;                                               %time spent in ATC maneuvering
Rman=Vman*tman;                                                                 %distance covered while maneu­vering due to ATC constraints (in nm)
Vcr = Vc;                                                                       %design cruise speed of an airplane in nautical miles per hour (kts)
%t_cr = (1.06*Rbl - Rcl - Rde + Rman)/Vcr;                                      %[hrs] time spent in cruise
t_bl = t_gm + t_cl + t_cr + t_de;                                               %[hrs] block time in hours
Vbl=Rbl/t_bl;                                                                   %block speed in nm/hr

U_ann_bl = (10^3)*(3.4546*(t_bl) + 2.994 - (12.289*((t_bl)^2) - 5.6626*(t_bl) + 8.964)^(1/2));      %annual utilization in block hours
                                                                                                    %U_ann_bl also expressed in flight hours U_ann_fli pg102
R_bl_ann = Vbl*U_ann_bl;                                                        %total annual block miles flown 

n_cj = n_captain + n_copilot + n_flightengineer;                                %number of crew members
%SAL_j = SAL_1989 * (CEF/CEF_1989);
SAL_j = SAL_1989 * (CEF);
%TEF_j = TEF_1990 * (CEF/CEF_1990);                                              %[USD] travel expense factor associated with each type of crew member
%Ccrew = n_cj*((1 + Kj)/Vbl) * (SAL_j/AH_j) + (TEF_j/Vbl);                       %[USD/nm] 
Ccrew = n_cj*((1 + Kj)/Vbl) * (SAL_j/AH_j);                                      %[USD/nm] (2024)


Cins = (f_ins_hull * AMP) / (U_ann_bl * Vbl);                                   %0.02 * DOC
DOCflt = Ccrew + C_en + Cins;                                                   %direct operating cost of flying in USD/nm (2024)

t_flt = t_cl + t_cr + t_de;
MHR_map_bl = (MHR_map_flt) * (t_flt / t_bl);                                    %number of airframe and systems maintenance man hours needed per block hour
%6.7 mhrs/blhr

Clab_ap = 1.03 * (MHR_map_bl * Rl_ap) / (Vbl);                                  %labor cost of airframe and systems in USD/nm (2024)
Clab_eng = 1.03 * (1.3) * (Ne) * (MHR_meng_bl) * (Rl_eng) / (Vbl);              %labor cost of engine(s) maintenan­ce in USD/nm (2024)
%MHR_flthr = MHR_map_flt + MHR_meng_flt; %0.3 pag 120
MHR_flthr = MHR_map_bl * t_bl / t_flt;
Cmat_ap = 1.03*(Cmat_apblhr)/Vbl;                                               %cost of maintenance materials for the airframe and systems in USD/nm (2024)
Cmat_eng = 1.03 * (Ne) * (C_mat_engblhr) / (Vbl);                       %cost of maintenance materials for the engines in USD/nm (2024)
Camb = 1.03 * (famb_lab * MHR_map_bl * Rl_ap + Ne * MHR_eng_bl * Rl_eng + famb_mat * (Cmat_apblhr + Ne * C_mat_engblhr )) / Vbl;   %applied maintenance burden in USD/nm
DOCmaint = Clab_ap + Clab_eng +	Cmat_ap + Cmat_eng + Camb;                      %direct operating cost of maintenance in USD/nm (2024)


Cdap = (Fdap * (AEP - Ne * EP - Np * PP - ASP)) / (DP_ap * U_ann_bl * Vbl);     %cost of airplane depreciation without engines and without propellers in USD/nm (2024)
Cdeng = (Fdeng * Ne * EP) / (DP_eng * U_ann_bl * Vbl);                          %cost of engine depreciation in USD/nm (2024)
Cdprp = (Fdprp * Np * PP) / (DP_prp * U_ann_bl * Vbl);                          %cost of depreciation of propellers in USD/nm (2024)
Cdav = (Fdav * ASP) / (DP_av * U_ann_bl * Vbl);                                 %cost of depreciation of avionics sys­tems in USD/nm (2024)
Cdapsp = (Fdapsp * Fapsp * (AEP - Ne * EP)) / (DP_apsp * U_ann_bl * Vbl);       %cost of depreciation of airplane spare parts in USD/nm (2024)
Cdengsp = (Fdengsp * Fengsp * Ne * EP * ESPPF) / (DP_engsp * U_ann_bl * Vbl);   %cost of depreciation of engine spare parts in USD/nm (2024)
Cbattery = (Fdeng * ASP) / (DP_eng * U_ann_bl * Vbl);
DOCdepr = Cdap + Cdeng + Cdprp + Cdav + Cdapsp + Cdengsp + Cbattery;                       %direct operating cost of depre­ciation in USD/nm (2024)


Caplf = 0.002 * Wto;                    %airplane landing fees in USD/lbs (2024)
Clf = Caplf/(Vbl*t_bl);                 %direct operating cost due to landing fees in USD/nm (2024)
Cnf  = Capnf/(Vbl*t_bl);                %cost of navigation fees in USD/nm (2024)
f_rt = 0.001 + (10^(-8)) * Wto;         % factor depending on airplane size
Crt = 0.01;                             %direct cost of registry taxes expres­ sed in USD/nm (2024)
DOClnr = Clf + Cnf + Crt;               %direct operating cost of landing fees, navigation fees and registry taxes in USD/nm

%DOCfin = 0.007 * DOC;                   %direct operating cost of financing in USD/nm

%DOC = DOCflt + DOCmaint + DOCdepr + DOClnr + DOCfin;
DOC = 1.007*(DOCflt + DOCmaint + DOCdepr + DOClnr);% / (1 - 0.7 - f_rt);             %USD/nm

DOCfin = 0.007 * DOC;                      

%%%%%%%%%%%%%%%%%%%%%%%%%%%% INDIRECT COST %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% IOCpax = Cmls + Cpax_ins + Ccat + Cpax_gen;	       %indirect operating cost for passen­ ger services. expressed in USD/nm
% 
% IOCsta = Cmgef + Cdgef;                              %indirect cost for maintaining and depreciating ground equipment and ground facilities
% 
% IOCascf = Caps + Capctrl + Cfrt;                     %indirect operating cost for air­plane and traffic servicing
% 
% IOC = IOCpax + IOCsta + IOCascf + IOCpse + IOCgaa;

IOC = 0.3 * DOC;

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

C_ops_dir=DOC*R_bl_ann*Nyr;
C_ops_ind=IOC*R_bl_ann*Nyr;
C_ops = C_ops_dir * Nacq + C_ops_ind * Nacq

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

% pie graph
% Data
data = [DOCflt, DOCmaint, DOCdepr, DOClnr, DOCfin, IOC];

% Labels
labels = {'DOCflt', 'DOCmaint', 'DOCdepr', 'DOClnr', 'DOCfin', 'IOC'};

figure;
h = pie(data);


delete(findall(gcf, 'Type', 'text'));
colors = get(findobj(h, 'Type', 'Patch'), 'FaceColor');
total = sum(data);
percentValues = data / total * 100;

% Add legend
legendEntries = cell(size(labels));
for i = 1:length(labels)
    legendEntries{i} = sprintf('%s: %.2f%%', labels{i}, percentValues(i));
end

legend(legendEntries, 'Location', 'eastoutside');

title('Operative cost - data percentage');

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

C_annual = C_ops/5;

C_mn = C_annual/R_bl_ann;

C_mission = C_mn * Rbl;

C_mission_time = C_mission/t_bl

C_per_pass = C_mission / 4

C_per_nm_per_pass = C_per_pass / Rbl

C_bl = C_annual / U_ann_bl;

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

% pie graph
% Data
data = [DOC, IOC];

% Labels
labels = {'DOC', 'IOC'};

figure;
h = pie(data);


delete(findall(gcf, 'Type', 'text'));
colors = get(findobj(h, 'Type', 'Patch'), 'FaceColor');
total = sum(data);
percentValues = data / total * 100;

% Add legend
legendEntries = cell(size(labels));
for i = 1:length(labels)
    legendEntries{i} = sprintf('%s: %.2f%%', labels{i}, percentValues(i));
end

legend(legendEntries, 'Location', 'eastoutside');

title('Operative cost - data percentage');

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