% Script to estimate the amount of fuel burned in the cruise segment
% Programmer: A. Trani
% date: February 2023
% Version 1.5
clear
clc

% Given:
%
%  Temperature off ISA 
% Aircraft model to use
%
% Functions used by this program: drag03.m and densityOffISA.m
% 
% % ______________ Cruise problem ___________

    % Load the aircraft model and delta temperature form ISA
    deltaTemp = 0;                      % 0 = ISA
    Boeing_787Class

    % Conversion factors, cruise altitude and Mach number

    kmToSM = 1.609;              % km to statute miles
    kmToNM = 1.852;              % km to nautical miles
    hcruise = 36000/3.28;        % Altitude in meters 
    mach_cruise = 0.82;          % Desired cruise Mach number 
    
    % This new procedure divives the cruise phase into a numerical
    % integration procedure with noSteps as the number of steps to be used
    % in the integration process
    
    initialCruiseMass = 220000;                         % initial mass at cruise (kg)
    currentCruiseMass(1) = initialCruiseMass;  % Define the first cell of a new array with cruise mass
    currentDistanceTraveled(1) = 0;
    cruiseDistance = 3500;                                 % nm from TOC to TOD
    noSteps = 1000;                                              % Number of steps to estimate fuel burn
    
    cruiseDistanceInterval = cruiseDistance*    kmToNM*1000 /noSteps;      % meters 
    phi = 0.0;  % bank angle 
    
    % Preallocate the variables drag, LD, cl and cd (speed up the
    % computations)

    drag = zeros(1,noSteps);
    LD = zeros(1,noSteps);
    cl = zeros(1,noSteps);
    cd = zeros(1,noSteps);
    
    % Start a loop to calculate the atmospheric and aerodynamic parameters
    
    for i=1:noSteps
    
            % Calculate drag at cruise speed using the drag03.m function
            % Calculate the atmospheric parameters: speed of sound,
            % density, temperature and density altitude using
            % densityOffISA.m function

            [speedOfSound,density,temp,dAlt] = densityOffISA(deltaTemp,hcruise);
            [drag(i),LD(i),cl(i),cd(i)]=drag_03(mach_cruise,density,speedOfSound,currentCruiseMass(i),macht,Cdoct,phi);   % find drag 
            
            % Calculate the product of Mach * L/D ratio
            MLD(i) = LD(i) .* mach_cruise;
            
            % Calculate fuel flow in cruise
            fuelFlowCruise(i) = drag(i) * tsfc;          % fuel flow in cruise (N/s)

            % Calculation of time in cruise 
            vtas = mach_cruise * speedOfSound;                         % true airspeed (m/s)
            cruiseTimeInterval = cruiseDistanceInterval / vtas ;     % time interval per cruise segement in seconds
            cruiseMassConsumed(i) = fuelFlowCruise(i) * cruiseTimeInterval / g;                   % units are kilograms

            currentCruiseMass(i+1) = currentCruiseMass(i) - cruiseMassConsumed(i);           % initial mass at cruise (kg)
            currentDistanceTraveled(i+1) = currentDistanceTraveled(i) + cruiseDistanceInterval; 

    end
    
    % Collect terms and report to Comand window
    
    cruiseFuel = sum(cruiseMassConsumed);
    cruiseDistance = cruiseDistanceInterval*noSteps;
    cruiseTime = cruiseTimeInterval*noSteps;
    finalMass = initialCruiseMass - cruiseFuel;
   
%     clc
%     
%     disp(' ')
%     disp([' Cruise Fuel (kg)  ',num2str(cruiseFuel)])
%     disp([' Cruise Time (minutes)  ',num2str((cruiseTimeInterval*noSteps)/60)])
%     disp([' Cruise Distance (nm)  ',num2str((cruiseDistanceInterval*noSteps)/(1000*kmToNM))])
%     
    clc
    
    disp(' ')
    disp([' Computational Steps   ',num2str(noSteps)])
    disp(' ')
    disp([' Cruise Distance (nm)  ',num2str((cruiseDistanceInterval*noSteps)/(1000*kmToNM))])
    disp([' Cruise Fuel (kg)  ',num2str(cruiseFuel)])
    disp([' Cruise Time (minutes)  ',num2str((cruiseTimeInterval*noSteps)/60)])
    disp([' Initial Mass (kg)  ',num2str(initialCruiseMass)])
    disp([' Final Mass (kg)  ',num2str(finalMass)])


    % Make some plots to showcase how the aerodynamic parameters change as
    % a function of distance traveled
    
    % Plot of drag vs distance traveled (km)
    figure
    plot(currentDistanceTraveled(2:end)/1000,drag,'o--')
    xlabel('Distance Traveled (km)','fontsize',24)
    ylabel('Drag (N)','fontsize',24)
    grid
    
    figure
    plot(currentDistanceTraveled(2:end)/1000,LD,'^r--')
    xlabel('Distance Traveled (km)','fontsize',24)
    ylabel('L/D Ratio','fontsize',24)
    grid
    
    figure
    plot(currentDistanceTraveled(2:end)/1000,MLD,'ob--')
    xlabel('Distance Traveled (km)','fontsize',24)
    ylabel('M*L/D Ratio','fontsize',24)
    grid
    