% Script to use calculate the drag and thrust of the aircraft for a given
% altitude and Mach number condition
% conditions

% Code assumes the drag polar is of the form:

% Cd = Cdo + Cdi

% Programmer: Toni Trani
% date: august 26, 2021

% calls:  drag03.m (to calculate the drag produced by the aircraft)
% calls:  thrust_calculationNoLoss.m (to calculate thrust)
% calls: isan.m (to calculate the atmospheric conditions)

clear all
global macht Cdoct phi neng

% Load the aircraft model
% Boeing738_class
% B787_class
B777_class

% Define the aircraft parameters 
% mass = 70000;          % kilograms for Boeing 737-800 class vehicle
mass = 290000;

conversion_ft2meters = 1/3.28;

% Define the aircraft altitude and speed conditions

altitude = 36000 * conversion_ft2meters;          % altitude (meters)
machNumber = 0.82;                                           % Mach number (dimensionless)
    
% Calculate the atmospheric conditions 

[speedOfSound,density,temperature] = isan(altitude);               % calculates speed of sound in m/s
        vtas = machNumber .* speedOfSound;                                                    % true airspeed (m/s)
   
% Estimate the drag coefficient, lift coefficient , drag and L/D ratio
        phi = 0;                % bank angle (rad)
        [drag,LD,lift_coefficient,drag_coefficient]=drag_03(machNumber,density,speedOfSound,mass,macht,Cdoct,phi);

% Calculate the thrust produced by each engine at the given altitude and Mach number

[thrust] = thrust_calculationNoLoss(thrust_table,mach_table,lapse_rate_factor,machNumber,density);
totalThrust = neng * thrust;

% Produce output numbers for Drag and Lift

disp(' ')
disp(['Altitude (m)  ',num2str(altitude)]);
disp(['True Airspeed (m/s)  ',num2str(vtas)]);
disp(['Lift Coeff. (dim)  ',num2str(lift_coefficient)]);
disp(['Drag Coeff. (dim)  ',num2str(drag_coefficient)]);
disp(['Drag  (N)  ',num2str(drag)]);

disp(' ')
disp(['Thrust per engine at  ', num2str(altitude), ' meters  ', '   in Newtons   ', num2str(thrust)]);
disp(['Total thrust in Newtons  ', num2str(totalThrust)]);
   