function [aircraft, state] = astSkyHogg(source)
%ASTSKYHOGG Fixed-wing object representation of the Sky Hogg aircraft
%
% [AIRCRAFT, STATE] = astSkyHogg()
%
% returns the fixed-wing object, AIRCRAFT, and fixed-wing cruise state, STATE,
% of a Sky Hogg aircraft.
% The Sky Hogg aircraft is a single engine, standard configuration aircraft
% modeled using 3 degrees of freedom and a thrust lookup table.
%
% Limitations:
%
% astSkyHogg requires Simulink(R).
%
% Examples:
%   Create a fixed-wing aircraft object and fixed-wing state object
%   using the Sky Hogg aircraft model
%     [aircraft, state] = astSkyHogg()
%
%   Create a fixed-wing aircraft object and fixed-wing state object
%   using the Sky Hogg aircraft model imported from DATCOM
%     [aircraft, state] = astSkyHogg("DATCOM")

% References
%
%[1] Cannon, M, Gabbard, M, Meyer, T, Morrison, S, Skocik, M, Woods, D. "Swineworks D-200 Sky Hogg Design Proposal." AIAA(R)/General Dynamics Corporation Team Aircraft Design Competition, 1991-1992.
%[2] Turvesky, A., Gage, S., and Buhr, C., "Accelerating Flight Vehicle Design", MATLAB(R)? Digest, January 2007.
%[3] Turvesky, A., Gage, S., and Buhr, C., "Model-based Design of a New Lightweight Aircraft", AIAA paper 2007-6371, AIAA Modeling and Simulation Technologies Conference and Exhibit, Hilton Head, South Carolina, Aug. 20-23, 2007.

% Copyright 2020-2021 The MathWorks, Inc.

arguments
    source (1,1) string {mustBeMember(source, ["DATCOM", "Demo"])} = "Demo"
end

switch source
    case "DATCOM"
        % values taken from the datcom file
        
        allData = datcomimport('astSkyHoggDatcom.out', false, 0);
        skyHoggData = allData{1};
        skyHoggData.cyb = fillmissing(skyHoggData.cyb, "previous", "MissingLocations", skyHoggData.cyb == 99999);
        skyHoggData.cnb = fillmissing(skyHoggData.cnb, "previous", "MissingLocations", skyHoggData.cnb == 99999);
        skyHoggData.clq = fillmissing(skyHoggData.clq, "previous", "MissingLocations", skyHoggData.clq == 99999);
        skyHoggData.cmq = fillmissing(skyHoggData.cmq, "previous", "MissingLocations", skyHoggData.cmq == 99999); 
        
        skyHogg = Aero.FixedWing();
        skyHogg.Properties.Name = "SkyHogg";
        skyHogg.DegreesOfFreedom = "3DOF";
        [aircraft, state] = datcomToFixedWing(skyHogg, skyHoggData);
        
        aircraft.Properties.Type = "General Aviation";
        
        aircraft.Thrusts(1).Properties.Name = "Propeller";
        aircraft.Thrusts(1).Coefficients.MultiplyStateVariables = "off";
        aircraft.Thrusts(1).Coefficients.NonDimensional = "off";
        aircraft.Thrusts(1).Coefficients.ReferenceFrame = "Body";
        
        aircraft.Surfaces.Properties.Name = "Delta";
        aircraft = aircraft.update();
        
        aircraft.Surfaces.MaximumValue = 20;
        aircraft.Surfaces.MinimumValue = -20;
        aircraft.DegreesOfFreedom = "3DOF";
        
        % thrust
        % thrust values taken from the asbSkyHogg simulink model
        Thrust = convforce(2512.574, 'N', 'lbf');
        prop = Simulink.LookupTable;
        prop.Table.Value = Thrust*(1+tanh(-3:3));
        prop.Breakpoints(1).Value = (1/6)*(3+(-3:3));
        prop.Breakpoints(1).FieldName = "Propeller";
        
        aircraft = aircraft.setCoefficient("CX", "Zero", {prop}, "Component", "Propeller");

        state = setupControlStates(state, aircraft);
        state = state.setState(["Propeller", "Delta"], [0.54, 8.6517e-04]);
        state.Mass = 1299.214;
        state.Inertia{:,:} = [5787.969 0 117.64;0 6928.93 0;-117.64 0 11578.329];
        
    case "Demo"
        % values taken from the 3DOF Sky Hogg Longitudinal Airframe in the
        % flightcontrol support package
        
        aircraft = Aero.FixedWing;
        aircraft.AngleSystem = "Degrees";
        aircraft.Properties.Name = "SkyHogg";
        aircraft.Properties.Type = "General Aviation";
        aircraft.ReferenceArea = 20.9775;
        aircraft.ReferenceLength= 1.7526;
        aircraft.ReferenceSpan = 12.54;
        aircraft.DegreesOfFreedom = "3DOF";
        
        alphaLT = Simulink.LookupTable;
        alphaLT.Breakpoints(1).FieldName = "Alpha";
        alphaLT.Breakpoints(1).Value = [-16     -12     -8      -4      -2      0       2       4       8       12];
        
        alphaLTCD = copy(alphaLT);
        alphaLTCL = copy(alphaLT);
        alphaLTCM = copy(alphaLT);
        
        alphaLTCD.Table.Value = [0.1150    0.0790    0.0470    0.0310    0.0270    0.0270    0.0290    0.0340    0.0540    0.0890];
        alphaLTCL.Table.Value = [-1.4210   -1.0920   -0.6950   -0.3120   -0.1320    0.0410    0.2180    0.4020    0.7860    1.1860];
        alphaLTCM.Table.Value = [0.0775    0.0663    0.0530    0.0337    0.0217    0.0073   -0.0090   -0.0263   -0.0632   -0.1235];
        aircraft.Coefficients = Aero.FixedWing.Coefficient(1, "Values", ...
            {alphaLTCD; 0; alphaLTCL; 0; alphaLTCM; 0});
        
        aircraft.Surfaces(1).Properties.Name = "HorizontalStabilizer";
        aircraft.Surfaces(1).Surfaces(1).Properties.Name = "Elevator";
        aircraft.Surfaces(1).Surfaces(1).Coefficients.StateVariables = "Elevator";
        aircraft.Surfaces(1).Surfaces(1).Coefficients.MultiplyStateVariables = "off";
        aircraft.Surfaces(1).Surfaces(1).Controllable = "on";
        
        deltaLT = Simulink.LookupTable;
        deltaLT.Breakpoints(1).FieldName = "Elevator";
        deltaLT.Breakpoints(1).Value = [-20     -10     0 10 20];
        
        alphaDeltaLT = Simulink.LookupTable;
        alphaDeltaLT.Breakpoints(1).FieldName = "Alpha";
        alphaDeltaLT.Breakpoints(1).Value = [-16        -12     -8      -4      -2      0       2       4       8       12];
        alphaDeltaLT.Breakpoints(2).FieldName = "Elevator";
        alphaDeltaLT.Breakpoints(2).Value = [-20        -10     0 10 20];
        
        deltaLTCD = copy(alphaDeltaLT);
        deltaLTCL = copy(deltaLT);
        deltaLTCM = copy(deltaLT);
        
        deltaLTCD.Table.Value = [
            0.0046  0.0029  0 -0.0009 -0.0008
            0.0039  0.0024  0 -0.0005 -0.0001
            0.0034  0.002   0  0        0.0005
            0.0028  0.0016  0  0.0004   0.0011
            0.0025  0.0014  0  0.0006   0.0014
            0.0022  0.0012  0  0.0008   0.0016
            0.0019  0.001   0  0.001    0.0019
            0.0016  0.0008  0  0.0012   0.0022
            0.001   0.0003  0  0.00160  0.0028
            0.0004 -0.0001 0  0.00210  0.0034
            ];
        deltaLTCL.Table.Value = [-0.0510   -0.0380         0    0.0380    0.0520];
        deltaLTCM.Table.Value = [0.0842    0.0601   -0.0001   -0.0601   -0.0843];
        
        aircraft.Surfaces(1).Surfaces(1).Coefficients.Values = {deltaLTCD; 0; deltaLTCL; 0; deltaLTCM; 0};
        aircraft.Surfaces(1).Surfaces.MaximumValue = 20;
        aircraft.Surfaces(1).Surfaces.MinimumValue = -20;
        
        aircraft.Thrusts(1).Properties.Name = "Propeller";
        aircraft.Thrusts(1).Coefficients.StateVariables = "Propeller";
        aircraft.Thrusts(1).Coefficients.MultiplyStateVariables = "off";
        aircraft.Thrusts(1).Coefficients.NonDimensional = "off";
        aircraft.Thrusts(1).Coefficients.ReferenceFrame = "Body";
        
        % thrust
        Thrust = 2512.574;
        prop = Simulink.LookupTable;
        prop.Table.Value = Thrust*(1+tanh(-3:3));
        prop.Breakpoints(1).Value = (1/6)*(3+(-3:3));
        prop.Breakpoints(1).FieldName = "Propeller";
        
        aircraft = aircraft.setCoefficient("CX", "Propeller", {prop}, "Component", "Propeller");
        
        aircraft = aircraft.update();

        % same values used in the demo, converted to degrees
        state = Aero.FixedWing.State(1, "AngleSystem", "Degrees");
        h = 2000;
        [T, a, p, rho] = atmoscoesa(h);
        state.AltitudeMSL = h;
        state.Environment.Temperature = T;
        state.Environment.SpeedOfSound = a;
        state.Environment.Pressure = p;
        state.Environment.Density = rho;
        state.U = 98.5;
        state.W = 1.52;
        state.Theta = rad2deg(0.0154);
        state.Mass = 1299.214;
        state.Inertia{:,:} = [5787.969 0 117.64;0 6928.93 0;-117.64 0 11578.329];
        state.CenterOfGravity = [0 0 0];
        
        state = setupControlStates(state, aircraft);
        state = state.setState(["Propeller", "Elevator"], [0.54, 8.6517e-04]);
end
end

