% Script to estimate the life cycle cost of a low boom aircraft
% Method is adapted from L. Nicolai textbook
% Aircraft and Airship Design: Volume 1, AIAA 2012

% Nicolai's method employs equations originally developed at the
% RAND Corporation

% Script to estimate the development of cost of an aircaft
% Using Nicolai's method and modified DAPCA IV equations

% Function calculateHourlyRates estimates hourly rates for engineering and
% production given a year

% Programmer: A. Trani

clear
close all
clc

    % Enter the key variables (data shown is similar to a large twin engine commercial aircraft likem the Boeing 777-300ER):
    Qproduction             = 250:10:700;                      % unites produced over the life cycle
    MaxSpeed_knots     = 500;                      % knots - maximum speed at cruise altitude
    operating_empty_weight_lb = 370000;    % Pounds
    thrust_surrogate_lb  =  115000;                 % pounds per engine

% Global variables

profitMargin      = 1.15;                               % Assumes 15% profit margin in the aircraft development program
yearOfAnalysis   = 1998;                             % Year with cost values in Nicolai's equations
engineModel      = 2;                                   % 1 = Nicolai, 2 = Our model
noEnginesProduction = 2;                        
noEnginesDevelopment = 3;
amortizationMultiplier   = 2;
avionicsInstallationMultiplier = 1.2;           % 20% of the avionics cost to install during production
avionicsCost = 6e6;                                  % $ in 1998 based on a survey done by Trani using B/C Aviation Avionics Data

% Input variables:

% Wempty_lb                       = Empty weight (lb)
% MaxSpeed_knots              = Maximum speed in knots at best altitude
% Qproduction                     = Total number of aircraft produced

Qdevelopment                                       = 7;                      % Development aircraft 
yearOfAnalysis                                      = 1998;
QC_to_Manufacturing_Hours_Ratio      = 0.13;                 % Quality control to manufacturing hours ratio

% Aircraft data

% inletTemperature = 2660;                    % Degrees Rankine
noEngines = 2;
noEnginesLCC = 2;                                   % ratio of no engines used in the life cycle
    
% Start the aircraft development cost equations

Qtotal  = Qdevelopment + Qproduction;

% Airframe engineering cost

% Calculate the engineering hours

EHours_DTE            = 4.86 * operating_empty_weight_lb .^0.777 .* MaxSpeed_knots .^ 0.894 .* Qdevelopment .^ 0.163;
EHours_Total          = 4.86 * operating_empty_weight_lb .^0.777 .* MaxSpeed_knots .^ 0.894 .* Qtotal.^ 0.163; 
EHours_Production  = EHours_Total  - EHours_DTE ;

% Estimate the hourly rates for all four activities (cost)

[hourlyRateTooling,hourlyRateEngineering,hourlyRateManufacturing,hourlyRateQC] = calculateHourlyRates(yearOfAnalysis);

EngineeringDTE_Cost             = EHours_DTE * hourlyRateEngineering;
EngineeringProduction_Cost   = EHours_Total  * hourlyRateEngineering;
Engineering_DTE_Cost            = EHours_Production  + EngineeringProduction_Cost;

% Development support cost

DevelopmentCost = 66 * operating_empty_weight_lb .^0.63 .* MaxSpeed_knots .^ 1.3;           % in 1998 dollars

% Flight test and operations

FlightTestCost = 1852 * operating_empty_weight_lb .^0.325 .* MaxSpeed_knots .^ 0.822 .* Qdevelopment .^ 1.21;

% Tooling cost

ToolingHours_DTE     = 5.99 * operating_empty_weight_lb .^0.777 .* MaxSpeed_knots .^ 0.696 .* Qdevelopment .^ 0.263;
ToolingHours_Total   = 5.99 * operating_empty_weight_lb .^0.777 .* MaxSpeed_knots .^ 0.696 .* Qtotal .^ 0.263;
ToolingHours_Production  = ToolingHours_Total  - ToolingHours_DTE;

ToolingDTE_Cost             = ToolingHours_DTE * hourlyRateTooling;
ToolingProduction_Cost   = ToolingHours_Production  * hourlyRateTooling;
ToolingTotal_Cost           = ToolingHours_Total * hourlyRateTooling;

% Manufacturing labor cost

ManLaborHours_DTE            = 7.37 * operating_empty_weight_lb .^0.82 .* MaxSpeed_knots .^ 0.484 .* Qdevelopment .^ 0.641;
ManLaborHours_Total          = 7.37 * operating_empty_weight_lb .^0.82 .* MaxSpeed_knots .^ 0.484 .* Qtotal  .^ 0.641;
ManLaborHours_Production  = ManLaborHours_Total  - ManLaborHours_DTE;

ManLabor_DTE_Cost             = ManLaborHours_DTE  * hourlyRateManufacturing;
ManLabor_Production_Cost   = ManLaborHours_Production  * hourlyRateManufacturing;
ManLabor_Total_Cost           = ManLaborHours_Total * hourlyRateManufacturing;

% Quality control cost

QualityControl_Hours_DTE            = QC_to_Manufacturing_Hours_Ratio * ManLaborHours_DTE  ;
QualityControl_Hours_Total          = QC_to_Manufacturing_Hours_Ratio * ManLaborHours_Total;
QualityControl_Hours_Production  = QC_to_Manufacturing_Hours_Ratio * ManLaborHours_Production;

QualityControl_DTE_Cost             = QualityControl_Hours_DTE  * hourlyRateQC;
QualityControl_Production_Cost   = QualityControl_Hours_Total    * hourlyRateQC;
QualityControl_Total_Cost           = QualityControl_Hours_Production * hourlyRateQC;


% Manufacturing material cost (in 1998 dollars)

ManMaterial_DTE_Cost             = 16.39 * operating_empty_weight_lb .^0.921 .* MaxSpeed_knots .^ 0.621 .* Qdevelopment .^ 0.799;
ManMaterial_Production_Cost   = 16.39 * operating_empty_weight_lb .^0.921 .* MaxSpeed_knots .^ 0.621 .* Qproduction .^ 0.799;
ManMaterial_Total_Cost          =  ManMaterial_DTE_Cost +  ManMaterial_Production_Cost;
  

% Engine cost

% Original equation to find the engine cost using the RAND equations
% EngineCost = 2306 * (0.043 .* thrust_surrogate_lb + 243.3 * MaxMach + 0.969 .* inletTemperature -2228);

if engineModel == 1
    % Simplified expression in Nicolai
    EngineCost = 520 * thrust_surrogate_lb .^0.8356;
elseif engineModel == 2
    % Our own regression to find the engine cost based on a regression
    % Equation using modern engines (Trani 2020)
    EngineCost = 1.0428e-08*thrust_surrogate_lb.^3 - 0.0028093*thrust_surrogate_lb.^2 + 329.55*thrust_surrogate_lb - 3.8218e+5;        % dollars 1998
end
    
    DTE_EngineCost                   = EngineCost * noEnginesDevelopment *  Qdevelopment;
    ProductionEngineCost          = EngineCost * noEnginesProduction * Qproduction;
    LCC_EngineCost                   = ProductionEngineCost  + DTE_EngineCost;

% Avionics cost

AvionicsCost = avionicsCost * avionicsInstallationMultiplier;          % dollars in 1998
LCC_AvionicsCost                = AvionicsCost * Qtotal;
DTE_AvionicsCost                = AvionicsCost * Qdevelopment;
Production_AvionicsCost      = AvionicsCost * Qproduction;

% Summary

totalProgramCost = LCC_AvionicsCost  + LCC_EngineCost + ManMaterial_Total_Cost + ...
    QualityControl_Total_Cost + ManLabor_Total_Cost + ToolingTotal_Cost + FlightTestCost + ...
    DevelopmentCost + Engineering_DTE_Cost;

totalDTECost = DTE_AvionicsCost  + DTE_EngineCost  + ManMaterial_DTE_Cost + ...
    QualityControl_DTE_Cost  + ManLabor_DTE_Cost  + ToolingDTE_Cost + FlightTestCost + ...
    DevelopmentCost + Engineering_DTE_Cost;

UnitCost_DTE = totalDTECost ./ (Qproduction/amortizationMultiplier);

totalProductionCost = Production_AvionicsCost  + ProductionEngineCost  + ManMaterial_Production_Cost + ...
    QualityControl_Production_Cost  + ManLabor_Production_Cost  + ToolingProduction_Cost + ...
    EngineeringProduction_Cost;

UnitCost_Production = totalProductionCost/ Qproduction;

costPerAirframe_1998 = UnitCost_DTE + UnitCost_Production;

% Use CPI index to convert 1998 costs into 2020 as needed

yearVector = [1998 2000 2010 2020];
CPI_index = [1.000 1.044 1.341 1.596];

cpi_value = interp1(yearVector,CPI_index,2020);

% costPerAirframe_1998                 = totalProgramCost ./ Qtotal;
costPerAirframe_wProfit                 = costPerAirframe_1998 * profitMargin;          % profit margin
adjustedCostPer_Airframe_2020    = costPerAirframe_wProfit * cpi_value;

% Make a plot of aircraft unit cost versus units produced in the life cycle

plot(Qproduction,adjustedCostPer_Airframe_2020,'o--r')
xlabel('Aircraft Production','fontsize',24)
ylabel('Cost per Airframe ($2020)','fontsize',24)
grid




