Load Recorded Flight Data for Analysis

The recorded data contains the following flight parameters:

load('astflight.mat');

Extract Flight Parameters from Loaded Data

MATLAB® variables are created for angle of attack (alpha), sideslip angle (beta), body angular rates (omega), and altitude (alt) from recorded data. The convangvel function is used to convert body angular rates from radians per second (rad/s) to degrees per second (deg/s).

alpha = fltdata(:,2); 
beta  = fltdata(:,3);
omega = convangvel( fltdata(:,5:7), 'rad/s', 'deg/s' );
alt   = fltdata(:,10);

Compute True Airspeed from Indicated Airspeed

In this set of flight data, indicated airspeed (IAS) was recorded. Indicated airspeed (IAS) is displayed in the cockpit instrumentation. To perform calculations, true airspeed (TAS), the airspeed without measurement errors, is typically used.

Measurement errors are introduced through the pilot-static airspeed indicators used to determine airspeed. These measurement errors are density error, compressibility error and calibration error. Applying these errors to true airspeed results in indicated airspeed.

The following data is the airspeed calibration table for the airspeed indicator of the aircraft with zero flap deflection. The airspeed calibration table converts indicated airspeed (IAS) to calibrated airspeed (CAS) by removing the calibration error.

flaps0IAS = 40:10:140;
flaps0CAS = [43 51 59 68 77 87 98 108 118 129 140];

The indicated airspeed (IAS) from the flight and airspeed calibration table are used to determine the calibrated airspeed (CAS) for the flight.

CAS = interp1( flaps0IAS, flaps0CAS, fltdata(:,4) );

The atmospheric properties, temperature (T), speed of sound (a), pressure (P), and density (rho), are determined at altitude for standard day using the atmoscoesa function.

[T, a, P, rho]= atmoscoesa( alt );

Once the calibrated airspeed (CAS) and the atmospheric properties are determined, the true airspeed (Vt) can be calculated using the correctairspeed function.

Vt = correctairspeed( CAS, a, P, 'CAS', 'TAS' );

Import Digital DATCOM Data for Aircraft

Use the datcomimport function to bring the Digital DATCOM data into MATLAB. The units for this aerodynamic information are feet and degrees.

data = datcomimport( 'astflight.out', true, 0 );

It can be seen in the Digital DATCOM output file and examining the imported data that

$$C_{Y\beta },$$ $$C_{n\beta },$$ $$C_{lq},$$ and $$C_{mq}$$

have data only in the first alpha value. By default, missing data points are set to 99999 and filled using the "previous" method of fillmissing. The missing data points are filled with the values for the first alpha, since these data points are meant to be used for all alpha values.

data{1}.cyb = fillmissing(data{1}.cyb, "previous", "MissingLocations", data{1}.cyb == 99999);
data{1}.cnb = fillmissing(data{1}.cnb, "previous", "MissingLocations", data{1}.cnb == 99999);
data{1}.clq = fillmissing(data{1}.clq, "previous", "MissingLocations", data{1}.clq == 99999);
data{1}.cmq = fillmissing(data{1}.cmq, "previous", "MissingLocations", data{1}.cmq == 99999);

Interpolate Stability and Dynamic Derivatives at Flight Conditions

The stability and dynamic derivatives in the digital DATCOM structure are 3-D tables that are functions of Mach number, angle of attack in degrees, and altitude in feet. To perform 3-D linear interpolation (griddedInterpolant), the indices for the derivative tables are required to be monotonic full grid of points. Indices of this form are generated by the ndgrid function.

[mnum, alp, h] = ndgrid( data{1}.mach, data{1}.alpha, data{1}.alt );

Since the angular units of the derivatives are in degrees, the units of angle of attack (alpha) are converted from radians to degrees by the function convang.

alphadeg = convang( alpha, 'rad', 'deg' );

The Mach numbers for the flight are calculated by the function machnumber using speed of sound (a) and airspeed (Vt).

Mach = machnumber( convvel( [Vt zeros(size(Vt,1),2)], 'kts', 'm/s' ), a );

GriddedInterpolant can be used to linearly interpolate derivative tables to find static and dynamic derivatives at the flight conditions.

F = griddedInterpolant( mnum, alp, h, pagetranspose(data{1}.cd), 'linear');
cd = F(Mach, alphadeg, alt);
F = griddedInterpolant( mnum, alp, h, pagetranspose(data{1}.cyb), 'linear');
cyb = F(Mach, alphadeg, alt);
F = griddedInterpolant( mnum, alp, h, pagetranspose(data{1}.cl), 'linear');
cl = F(Mach, alphadeg, alt);
F = griddedInterpolant( mnum, alp, h, pagetranspose(data{1}.cyp), 'linear');
cyp = F(Mach, alphadeg, alt);
F = griddedInterpolant( mnum, alp, h, pagetranspose(data{1}.clad), 'linear');
clad = F(Mach, alphadeg, alt);

Compute Aerodynamic Coefficients

Once the derivatives are found for the flight conditions, aerodynamic coefficients can be calculated.

Reference lengths and areas used in the aerodynamic coefficient computation are extracted from the digital DATCOM structure.

cbar = data{1}.cbar;
Sref = data{1}.sref;
bref = data{1}.blref;

The angular units of the derivatives are in degrees, so the units of sideslip angle (beta) are converted from radians to degrees by the function convang.

betadeg  = convang( beta, 'rad', 'deg' );

To calculate the aerodynamic coefficients, the body angular rates (omega) need to be given in the stability axes, like the derivatives. The function dcmbody2stability generates the direction cosine matrix for body axes to stability axes (Tsb) when the sideslip angle (beta) is set to zero.

Tsb = dcmbody2stability( alpha );

The rate of change in angle of attack (alpha_dot) is also needed to find angular rates in stability axis (omega_stab). The function diff is used on alpha in degrees divided by data sample time (0.50 seconds) to approximate the rate of change in angle of attack (alpha_dot).

alpha_dot = diff( alphadeg/0.50 );

The last value of alpha_dot is held to keep the length of alpha_dot consistent with other arrays in this calculation. This is needed because the diff function returns an array that is one value shorter than the input.

alpha_dot = [alpha_dot; alpha_dot(end)];

The angular rates in stability axis (omega_stab) are computed for the flight data. The angular rates are reshaped into a 3-D matrix to be multiplied with the 3-D matrix for the direction cosine matrix for body axes to stability axes (Tsb).

omega_temp = reshape((omega - [zeros(size(alpha)) alpha_dot zeros(size(alpha))])',3,1,length(omega));

for k = length(omega):-1:1
    omega_stab(k,:) = (Tsb(:,:,k)*omega_temp(:,:,k))'; 
end

Compute the drag coefficient (CD), the side force coefficient (CY), and the lift coefficient (CL). The convvel function is used to get the units of airspeed (Vt) consistent with those of the derivatives.

CD = cd;
CY = ( cyb.*betadeg ) + ((( cyp.*omega_stab(:,1) )*bref )/( 2./convvel(Vt, 'kts', 'ft/s') ));
CL = cl + ((( clad.*alpha_dot )*cbar )/( 2./convvel(Vt, 'kts', 'ft/s') ));

Compute Forces

Aerodynamic coefficients for drag, side force and lift are used to compute aerodynamic forces.

Dynamic pressure is needed to calculate the aerodynamic forces. The function dpressure compute dynamic pressure from the airspeed (Vt) and density (rho). The convvel function is used to get the units of airspeed (Vt) consistent with those of density (rho).

qbar = dpressure( convvel( [Vt zeros(size(Vt,1),2)], 'kts', 'm/s' ), rho );

To find forces in body axes, the direction cosine matrix for stability axes to body axes (Tbs) is needed. Direction cosine matrix for stability axes to body axes (Tbs) is the transpose of direction cosine matrix for body axes to stability axes (Tsb). To take the transpose of a 3-D array, the pagetranspose function is used.

Tbs = pagetranspose(Tsb);

Looping through the flight data points, aerodynamic forces are computed and converted from stability to body axes. The convpres function is used to get the units of dynamic pressure (qbar) consistent with those of the reference area (Sref).

for k = length(qbar):-1:1
    forces_lbs(k,:) = Tbs(:,:,k)*(convpres(qbar(k),'Pa','psf')*Sref*[-CD(k); CY(k); -CL(k)]); 
end

A constant thrust is estimated in the body axes.

thrust = ones(length(forces_lbs),1)*[200 0 0];

The constant thrust estimate is added to aerodynamic forces and units are converted to metric.

forces = convforce((forces_lbs + thrust),'lbf','N');

Estimate G Forces

Using the calculated forces, estimate G forces during the flight.

Accelerations are estimated using the calculated forces and mass converted into kilograms using convmass. Accelerations are converted to G forces using convacc.

N = convacc(( forces/convmass(84.2,'slug','kg') ),'m/s^2',"G's");
N = N .* [1,1,-1];

G forces are plotted over the flight.

h1 = figure;
plot(fltdata(:,1), N);
xlabel('Time (sec)')
ylabel('G Force')
title('G Forces over the Flight')
legend('Nx','Ny','Nz','Location','Best')

close(h1);