Array indices must be positive integers or logical values.

My code worked fine before trying to incorporate the for loops. Then it keeps giving me this error and I don't know where it went wrong. I have attached my code in a txt file. Any help is appreciated!

clear; clc; close

sweepRange = 5:5:35; % My specified sweep angle in degrees to begin
ARRange = 6:0.1:12; % My specified aspect ratio to begin

%% For Loops for Sweep and AR
for i = 1:length(sweepRange)
    sweep = sweepRange(i);
    for j = 1:length(ARRange)
        AR = ARRange(j);

        % 2023 Design Specifications
        P = 210; % Number of passengers
        cargo = 4000; % Total cargo weight in lbs
        range = 3500; % Max scheduled range in nmi
        h_alt = 35000; % Initial cruise altitude in ft
        M_alt = 0.80; % Initial cruise Mach number
        TOFL = 6000; % Takeoff field length in ft at sea-level on a hot day (84 F)
        Approach = 135; % Approach speed in kts

        % Chosen Specifications
        lambda = 0.35; % Wing taper ratio
        airfoil = 2; % (1 for conventional and 2 for supercritical)
        EngineType = 2; % (1 for JT8D and 2 for JT9D)
        NumEngines = 2; % Specifying the number of engines
        EngineMount = 1; % (1 for wing-mounted engines and 2 for fuselage mounted engines)
        abreast = 6; % Number of seats abreast
        aisle = 1; % Number of aisles
        aluminum = 1; % (1 for aluminum and 2 for no aluminum)
        gamma = 1.4; % Ratio of specific heats for air
        R = 1718; % Gas constant for air in English units
        T = 394.08; % Temperature at 35000 ft in degrees Rankine
        sound = (sqrt(gamma*R*T))*(1/1.68781); % Speed of sound at 35000 ft in kts
        V_cruise = M_alt*sound; % Cruise speed in kts
        R_ao = range + 200 + (V_cruise)*(0.75); % Scheduled all out range in nmi

        % Range loop values
        range_total = 100; % Setting the value of all out range to trigger the loop
        FuelAdjust = 0; % An adjustment factor for fuel carried

        % Thrust loop values
        ThrustAdjust = 0; % An adjustment factor for thrust loading
        T_R_JT_IC = 11e4; % Setting an initial value for thrust required at cruise for JT8D/9D engine
        T_a_JT_cruise = 1;

        %% Outer Loop for Thrust
        while T_R_JT_IC > T_a_JT_cruise
            %% Loop for Range
            while abs(R_ao - range_total) > 10
                C_Lguess = 0.60; % Initial value chosen for lift coefficient
                C_LIC = 0.1; % Creating a system variable for initial cruise lift coefficient
                %% Section 1: Initial Cruise Lift Coefficient Calculations
                while abs(C_LIC - C_Lguess) >= 0.001
                    C_L = C_Lguess;
                    % Change in divergent Mach number
                    if airfoil == 1
                        DeltaM_Div = -0.2156511*(C_L)^2 - 0.1007080*(C_L) + 0.1209273; % Figure 2 (Conventional)
                    else
                        DeltaM_Div = 0.8579659*(C_L)^3 - 1.8221479*(C_L)^2 + 1.0693644*C_L - 0.1794305; % Figure 2 (Supercritical)
                    end
                    M_Div = M_alt + 0.004 - DeltaM_Div; % Divergent Mach number
                    % Thickness to chord ratio for different airfoil types and sweep angles
                    if airfoil == 1
                        if (sweep >= 0 && sweep < 10)
                            tc_ave = (-0.6356774*(M_Div) + 0.5735735)*((10 - sweep)/(10 - 0)) + (-0.6188624*(M_Div) + 0.5654024)*(1 - ((10 - sweep)/(10 - 0)));
                        elseif (sweep >= 10 && sweep < 15)
                            tc_ave = (-0.6188624*(M_Div) + 0.5654024)*((15 - sweep)/(15 - 0)) + (-0.5951630*(M_Div) + 0.5527027)*(1 - ((15 - sweep)/(15 - 10)));
                        elseif (sweep >= 15 && sweep < 20)
                            tc_ave = (-0.5951630*(M_Div) + 0.5527027)*((20 - sweep)/(20 - 15)) + (-0.5689636*(M_Div) + 0.5391453)*(1 - ((20 - sweep)/(20 - 15)));
                        elseif (sweep >= 20 && sweep < 25)
                            tc_ave = (-0.5689636*(M_Div) + 0.5391453)*((25 - sweep)/(25 - 20)) + (-0.5348202*(M_Div) + 0.5206442)*(1 - ((25 - sweep)/(25 - 20)));
                        elseif (sweep >= 25 && sweep < 30)
                            tc_ave = (-0.5348202*(M_Div) + 0.5206442)*((30 - sweep)/(30 - 25)) + (-0.5056037*(M_Div) + 0.5060521)*(1 - ((30 - sweep)/(30 - 25)));
                        elseif (sweep >= 30 && sweep < 35)
                            tc_ave = (-0.5056037*(M_Div) + 0.5060521)*((35 - sweep)/(35 - 30)) + (-0.4698165*(M_Div) + 0.4869880)*(1 - ((35 - sweep)/(35 - 30)));
                        elseif (sweep >= 35 && sweep < 40)
                            tc_ave = (-0.4698165*(M_Div) + 0.4869880)*((40 - sweep)/(40 - 35)) + (-0.4305368*(M_Div) + 0.4656397)*(1 - ((40 - sweep)/(40 - 35)));
                        elseif sweep == 40
                            tc_ave = -0.4305368*(M_Div) + 0.4656397;
                        end
                    end
                    if airfoil ~= 1
                        if (sweep >= 0 && sweep < 5)
                            tc_ave = (-38.2571833*(M_Div)^3 + 89.0718974*(M_Div)^2 - 69.7890623*(M_Div) + 18.4846319)*((5 - sweep)/(5 - 0)) + (-27.8560713*(M_Div)^3 + 66.1635208*(M_Div)^2 - 53.0102696*(M_Div) + 14.4005209)*(1 - ((5 - sweep)/(5 - 0)));
                        elseif (sweep >= 5 && sweep < 10)
                            tc_ave = (-27.8560713*(M_Div)^3 + 66.1635208*(M_Div)^2 - 53.0102696*(M_Div) + 14.4005209)*((10 - sweep)/(10 - 5)) + (-19.1627368*(M_Div)^3 + 47.0249724*(M_Div)^2 - 39.0047747*(M_Div) + 11.0001550)*(1 - ((10 - sweep)/(10 - 5)));
                        elseif (sweep >= 10 && sweep < 15)
                            tc_ave = (-19.1627368*(M_Div)^3 + 47.0249724*(M_Div)^2 - 39.0047747*(M_Div) + 11.0001550)*((15 - sweep)/(15 - 10)) + (-21.8271021*(M_Div)^3 + 53.5679708*(M_Div)^2 - 44.3673843*(M_Div) + 12.4756828)*(1 - ((15 - sweep)/(15 - 10)));
                        elseif (sweep >= 15 && sweep < 20)
                            tc_ave = (-21.8271021*(M_Div)^3 + 53.5679708*(M_Div)^2 - 44.3673843*(M_Div) + 12.4756828)*((20 - sweep)/(20 - 15)) +(-15.2177170*(M_Div)^3 + 38.8450823*(M_Div)^2 - 33.5369749*(M_Div) + 9.8605802)*(1 - ((20 - sweep)/(20 -15)));
                        elseif (sweep >= 20 && sweep < 25)
                            tc_ave = (-15.2177170*(M_Div)^3 + 38.8450823*(M_Div)^2 - 33.5369749*(M_Div) + 9.8605802)*((25 - sweep)/(25 - 20)) + (-18.4322440*(M_Div)^3 + 47.5880517*(M_Div)^2 - 41.4576687*(M_Div) + 12.2611522)*(1 - ((25 - sweep)/(25 - 20)));
                        elseif (sweep >= 25 && sweep < 30)
                            tc_ave = (-18.4322440*(M_Div)^3 + 47.5880517*(M_Div)^2 - 41.4576687*(M_Div) + 12.2611522)*((30 - sweep)/(30 - 25)) + (-16.7368110*(M_Div)^3 + 45.0535774*(M_Div)^2 - 40.9221296*(M_Div) + 12.6083969)*(1 - ((30 - sweep)/(30 - 25)));
                        elseif (sweep >= 30 && sweep < 35)
                            tc_ave = (-16.7368110*(M_Div)^3 + 45.0535774*(M_Div)^2 - 40.9221296*(M_Div) + 12.6083969)*((35 - sweep)/(35 - 30)) +(-37.2527592*(M_Div)^3 + 99.9053005*(M_Div)^2 - 90.0209481*(M_Div) + 27.3326317)*(1 - ((35 - sweep)/(35 - 30)));
                        elseif (sweep >= 35 && sweep < 40)
                            tc_ave = (-37.2527592*(M_Div)^3 + 99.9053005*(M_Div)^2 - 90.0209481*(M_Div) + 27.3326317)*((40 - sweep)/(40 - 35)) + (-403.3066991*(M_Div)^3 + 1060.1835747*(M_Div)^2 - 930.6486246*(M_Div) + 272.9330195)*(1 - ((40 - sweep)/(40 - 35)));
                        elseif sweep == 40
                            tc_ave = -403.3066991*(M_Div)^3 + 1060.1835747*(M_Div)^2 - 930.6486246*(M_Div) + 272.9330195;
                        end
                    end
                    var = (cosd(sweep))^2 *(tc_ave)^2 *AR; % Figure 3 variable
                    C_LmaxLD = 103.2508180*(var)^3 - 65.1539504*(var)^2 + 16.2173873*(var) + 2.0142218; % Max C_L during landing
                    C_LmaxTO = 81.6463897*(var)^3 - 63.5364168*(var)^2 + 16.9175815*(var) + 1.0394459; % Max C_L during takeoff
                    rho_SL84 = 2.265677968e-3; % Sea-level density at 84 degrees Fahrenheit
                    rho_SL = 2.3769e-3; % Standard sea-level density
                    sigma = rho_SL84/rho_SL; % Density ratio at sea-level
                    WS_LD = (Approach/1.3)^2 *((sigma*C_LmaxLD)/296); % Wing loading during landing in lbs/ft^2
                    if EngineType == 1
                        Wf_WTO_JT8D = (5.1856130e-13)*(R_ao)^3 - (1.1830472e-8)*(R_ao)^2 + (1.2480300e-4)*(R_ao) + 6.9173913e-3 + FuelAdjust; % Figure 4 equation
                        Wf_WTO = Wf_WTO_JT8D;
                        x = 0; % Coefficient in step 10
                        WS_TO = (WS_LD)/(1-(x*Wf_WTO)); % Wing loading at takeoff in lbs/ft^2
                        WS_IC = 0.965*WS_TO; % Wing loading at initial cruise in lbs/ft^2
                        p_alt = 499.34; % Pressure at 35000 ft in lbs/ft^2
                        p_SL = 2116.2; % Pressure at sea-level in lbs/ft^2
                        delta = p_alt/p_SL; % Pressure ratio
                        C_LIC = (WS_IC)/(1481*delta*(M_alt)^2); % Initial cruise C_L
                        C_Lfinal = C_L; % Final converging C_L
                    elseif EngineType ~= 1
                        Wf_WTO_JT8D = (5.1856130e-13)*(R_ao)^3 - (1.1830472e-8)*(R_ao)^2 + (1.2480300e-4)*(R_ao) + 6.9173913e-3 + FuelAdjust; % Figure 4 equation
                        SFC_JT9D = 0.61; % Specific fuel consumption for JT9D engine at 85% max continuous cruise thrust
                        SFC_JT8D = 0.78; % Specific fuel consumption for JT8D engine at 85% max continuous cruise thrust
                        Wf_WTO = Wf_WTO_JT8D *(SFC_JT9D/SFC_JT8D); % Fuel weight over takeoff weight for JT9D engine
                        x = 0; % Coefficient in step 10
                        WS_TO = (WS_LD)/(1-(x*Wf_WTO)); % Wing loading at takeoff in lbs/ft^2
                        WS_IC = 0.965*WS_TO; % Wing loading at initial cruise in lbs/ft^2
                        p_alt = 499.34; % Pressure at 35000 ft in lbs/ft^2
                        p_SL = 2116.2; % Pressure at sea-level in lbs/ft^2
                        delta = p_alt/p_SL; % Pressure ratio
                        C_LIC = (WS_IC)/(1481*delta*(M_alt)^2); % Initial cruise C_L
                        C_Lfinal = C_L; % Final converging C_L
                    end
                    if(C_LIC > C_Lfinal)
                        C_Lguess = C_L + 0.001;
                    else
                        C_Lguess = C_L - 0.001;
                    end
                end
                disp(['Lift coefficient converges at  ' num2str(C_Lfinal)]);

                %% Section 2: Thrust Loading Calculations
                if NumEngines == 2
                    K = 28.1629049*(TOFL/1000) - 8.0256314; % Figure 5 value
                elseif NumEngines == 3
                    K = 31.3478489*(TOFL/1000) - 6.6448483;
                elseif NumEngines == 4
                    K = 32.3531006*(TOFL/1000) + 2.2111345;
                end
                WT_07VLO = (K*sigma*C_LmaxTO)/(WS_TO); % Thrust loading at 70% L/O velocity
                V_LO = 1.2*sqrt((296*WS_TO)/(sigma*C_LmaxTO)); % Liftoff velocity in kts
                M_LO = V_LO/(661*sqrt(sigma)); % Liftoff Mach number
                M_07LO = 0.7*M_LO; % 70 percent liftoff Mach number
                if EngineType == 1
                    T_SLS = 14500; % Sea-level static thrust for JT8D engine in lbs
                    T_M07LO = 7482.6479132*(M_07LO)^2 - 8676.7374528*(M_07LO) + 14819.6019504; % Sea-level dry takeoff thrust at 70% liftoff Mach number in lbs
                elseif EngineType == 2
                    T_SLS = 45500; % Sea-level static thrust for JT9D in lbs
                    T_M07LO = 1000*(37.4386807*(M_07LO)^2 - 47.3982379*(M_07LO) + 45.3366105); % Sea-level dry takeoff thrust at 70% liftoff Mach number in lbs
                end
                WT = WT_07VLO*(T_M07LO/T_SLS) + ThrustAdjust; % Thrust loading
                %                 disp(['Thrust loading is ' num2str(WT)]);

                %% Section 3: Takeoff Weight Calculations
                eta = 3.75; % Ultimate load factor
                tc_bar = tc_ave + 0.03; % Adjusted t/c
                % Constant k_w for calculating the wing weight
                if EngineMount == 1
                    k_w = 1.00;
                elseif EngineMount ~= 1
                    k_w = 1.03;
                end
                % Wing weight
                W_w = (0.00945 *(AR)^0.8 *(1 + lambda)^0.25 *k_w *(eta)^0.5 )/((tc_bar)^0.4 *cosd(sweep) *(WS_TO)^0.695 );
                % Fuselage weight
                k_f = 11.5; % Constant for calculating fuselage weight
                length = (3.76*(P/abreast) + 33.2); % Fuselage length in ft
                diameter = (1.75*abreast + 1.58*aisle + 1.0); % Fuselage diameter in ft
                W_fus = 0.6727 *k_f *(length)^0.6 *(diameter)^0.72 *(eta)^0.3;
                % Landing gear weight
                W_LG = 0.040;
                % Nacelles and pylons weight
                W_NP = (0.0555)/(WT);
                % Tail surface weight
                if EngineMount == 1
                    k_ts = 0.17;
                elseif EngineMount == 2
                    k_ts = 0.25;
                end
                W_TS = k_ts*W_w;
                % Power plant weight
                W_PP = 1/(3.58*WT);
                % Fuel weight
                W_f = 1.0275*(Wf_WTO);
                TankWeight = 0.0175*W_f; % Fuel tank weight
                UnuseableFuel = 0.01*W_f; % Wasted fuel weight
                % Payload weight
                PassengerWeight = 165*P; % Weight per passenger in lbs
                BaggageWeight = 50*P; % Baggage weight per passenger in lbs
                W_PL = PassengerWeight + BaggageWeight + cargo;
                % Fixed equipment weight
                crew = 2; % Number of flight crew
                attendant = ceil(P/50); % Number of flight attendants
                W_fe = 132*P + 300*NumEngines + 0.035 + 260*crew + 170*attendant;
                % Combining like terms
                a = W_w + W_TS; % (W_TO)^1.195 terms
                b = W_fus; % (W_TO)^0.235 terms
                c = W_LG + W_NP + W_PP + W_f + 0.035 -1; % (W_TO) terms
                d = W_PL + (W_fe - 0.035); % Constant terms
                W_TO = 1e6; % Guessing a value for takeoff weight in lbs
                W = a*(W_TO)^1.195 + b*(W_TO)^0.235 + c*(W_TO) + d;
                % While loop to iterate W_TO value until convergance
                while abs(W) > 10
                    W = a*(W_TO)^1.195 + b*(W_TO)^0.235 + c*(W_TO) + d;
                    if W < 10
                        W_TO = W_TO - 10;
                    else
                        W_TO = W_TO + 10;
                    end
                end
                W_TO_Al = W_TO - 0.06*W_w - 0.06*W_TS - 0.06*W_fus; % Advanced technology factors for aluminum
                if aluminum == 1
                    TakeoffWeight = W_TO_Al;
                elseif aluminum == 2
                    TakeoffWeight = W_TO;
                end
                disp(['Takeoff weight is ' num2str(TakeoffWeight) ' lbs'])
                S = TakeoffWeight/WS_TO; % Wing planform area in ft^2
                b = sqrt(AR*S); % Wing span in ft
                mac = S/b; % Approximate mean aerodynamic chord in ft
                T = TakeoffWeight/WT; % Total sea-level static thrust produced by the engine in lbs
                T_e = T/NumEngines; % Sea-level static thrust per engine in lbs
                %                 disp(['Sea-level static thrust per engine is ' num2str(T_e) ' lbs']);
                %                 disp(['Wing planform area is ' num2str(S) ' ft^2'])

                %% Section 4: Drag Calculations
                nu_30000 = 3.4884e-4; % Kinematic viscosity in ft^2/s at 30000 ft
                sound_30000 = 995.2917522; % Speed of sound in ft/s at 30000 ft
                M_30000 = 0.5; % Mach number at 30000 ft
                RN_30000 = (M_30000*sound_30000)/nu_30000; % Reynolds number per unit length at 30000 ft
                RN_w = RN_30000*mac; % Wing Reynolds number
                RN_fus = RN_30000*length; % Fuselage Reynolds number
                c_f_w = (81.7340305*(RN_w)^(-0.1951868))/1e3; % Wing skin friction coefficient
                c_f_fus = (81.7340305*(RN_fus)^(-0.1951868))/1e3; % Fuselage skin friction coefficient
                S_wet_w = 2*1.02*(S-diameter*30); % Wing wetted area in ft^2
                S_wet_fus = 0.9*pi*diameter*length; % Fuselage wetted are in ft^2
                z_w = ((2 - M_30000^2)*cosd(sweep))/sqrt(1 - M_30000^2 *(cosd(sweep))^2); % z from Shevell Fig 11.3
                k_wing = 1 + z_w*tc_ave + 100*(tc_ave)^4; % Wing form factor
                f_w = k_wing*c_f_w*S_wet_w; % Wing flat plate drag area in ft^2
                fineness = length/diameter; % Fuselage fineness ratio
                k_fus = -0.0016263*(fineness)^3 + 0.0408285*(fineness)^2 - 0.3672062*(fineness) + 2.3068476; % Fuselage form factor
                f_fus = k_fus*c_f_fus*S_wet_fus; % Fuselage flat plate drag area in ft^2
                f_ts = 0.35*f_w; % Tail surfaces flat plate drag area
                k_n = 1.25; % Nacelle form factor
                c_f_n = c_f_w; % Nacelle skin friction coefficient
                S_wet_n = 2.1*sqrt(T_e)*NumEngines; % Nacelle wetted area in ft^2
                f_n = k_n*c_f_n*S_wet_n; % Nacelle flat plate drag area in ft^2
                f_p = 0.2*f_n; % Pylons flat plate drag area in ft^2
                f_total = (f_w + f_fus + f_ts + f_n + f_p)*1.06; % Total flat plate drag area
                C_D0 = f_total/S; % Profile drag coefficient
                e = 1/(1.035 + 0.38*C_D0*pi*AR); % Oswald efficiency factor
                %                 disp(['Flat plate drag area is ' num2str(f_total) ' ft^2']);
                %                 disp(['Profile drag coefficient is ' num2str(C_D0)])
                %                 disp(['Oswald efficiency factor is ' num2str(e)]);

                %% Section 5: Range During Climb Calculations
                rho_20000 = 1.2673e-3; % Density at 20000 ft
                sigma_climb = rho_20000/rho_SL; % Climb density ratio
                W_climb = (1.965/2)*TakeoffWeight; % Average climb weight in lbs
                V_climb = (1.3*12.9)*((f_total*e)^(-1/4)) *sqrt(W_climb/(sigma_climb*b)); % Climb airspeed in kts
                a_20000 = 1037.380661; % Speed of sound at 20000 ft in ft/s
                M_climb = (V_climb*1.688)/a_20000; % Climb Mach number
                T_R_climb = ((sigma_climb*f_total*V_climb^2)/296) + (94.1/(sigma_climb*e*V_climb^2))*(W_climb/b)^2; % Thrust required during climb in lbs
                if EngineType == 1
                    SFC_20000 = (0.2211680)*(M_climb)^3 - (0.3020208)*(M_climb)^2 + (0.4227834)*(M_climb) + 0.5723053; % An estimate for SFC at 20000 ft
                    T_a_20000 = -1819.9312960*(M_climb)^3 + 5223.9814259*(M_climb)^2 - 4461.5269824*(M_climb) + 7509.1235238; % JT8D Max Climb Thrust at 20000 ft
                    T_a_climb = (T_e/T_SLS)*T_a_20000; % Available thrust per engine during climb
                elseif EngineType == 2
                    SFC_15000 = 0.0759755*(M_climb)^2 + 0.3476766*(M_climb) + 0.3342225; % An estimate for SFC at 15000 ft
                    SFC_25000 = 0.3702403*(M_climb) + 0.3417182; % An estimate for SFC at 25000 ft
                    SFC_20000 = (SFC_15000 + SFC_25000)/2; % An estimate for SFC at 20000 ft
                    T_a_15000 = 1e3*(-9.7850170*(M_climb)^3 + 23.5203048*(M_climb)^2 -25.1318476*(M_climb) + 27.0910853); % JT9D Max Climb Thrust at 15000 ft
                    T_a_25000 = 1e3*(-5.0119266*(M_climb)^3 + 11.6581282*(M_climb)^2 -9.9497828*(M_climb) + 17.3467089); % JT9D Max Climb Thrust at 25000 ft
                    T_a_20000 = (T_a_15000 + T_a_25000)/2; % JT9D Max Climb thrust at 20000 ft
                    T_a_climb = (T_e/T_SLS)*T_a_20000; % Available thrust per engine during climb
                end
                RoC = (101*(T_a_climb*NumEngines - T_R_climb)*V_climb)/W_climb; % R/C in ft/min
                t_climb = h_alt/RoC; % Time to climb in min
                range_climb = V_climb*(t_climb/60); % Climb range in nmi
                W_f_climb = T_a_climb*SFC_20000*(t_climb/60); % Fuel used during climb in lbs
                %                 disp(['Available thrust per engine during climb is ' num2str(T_a_climb) ' lbs'])
                %                 disp(['Rate of climb is ' num2str(RoC) ' ft/min'])
                %                 disp(['Time to climb is ' num2str(t_climb) ' min'])

                %% Section 6: Range During Cruise Calculations
                W_0 = W_TO - W_f_climb; % Initial cruise weight in lbs
                W_1 = (1 - Wf_WTO)*W_TO; % Final cruise weight in lbs
                C_LavgCruise = ((W_0 + W_1)/(2*S))/(1481*delta*M_alt^2); % Average C_L during cruise
                C_Di = (C_LavgCruise^2)/(pi*AR*e); % Induced drag coefficient
                DeltaC_DC = 0.001; % Compressibility drag coefficient
                C_D = C_D0 + C_Di + DeltaC_DC; % Total drag Coefficient
                LiftToDrag = C_LavgCruise/(C_D); % Lift to drag ratio
                T_R_cruise = ((W_0 + W_1)/2)/LiftToDrag; % Total thrust required at cruise
                T_R_JT = (T_R_cruise/NumEngines)*(T_SLS/T_e); % Thrust required for JT8D/9D engine
                % SFC value at given M and cruise altitude
                if EngineType == 1
                    SFC_35000 = (7.8188283e-13)*(T_R_JT)^4 - (4.3009761e-9)*(T_R_JT)^3 + (8.7558319e-6)*(T_R_JT)^2 - (8.0052903e-3)*(T_R_JT) + 3.7025121;
                elseif EngineType == 2
                    SFC_35000 = (-4.5023006e-12)*(T_R_JT)^3 + (7.3349955e-8)*(T_R_JT)^2 - (4.1068768e-4)*(T_R_JT) + 1.4133773;
                end
                range_cruise = ((M_alt*sound)/SFC_35000)*LiftToDrag*log(W_0/W_1); % Cruise range in nmi
                range_total = range_climb + range_cruise; % Total range in nmi
                Wf_WTO = (W_0 - W_1)/W_TO; % Recalculating the fuel fraction
                if range_total < R_ao
                    FuelAdjust = FuelAdjust + 0.001;
                else
                    FuelAdjust = FuelAdjust - 0.001;
                end
            end
            %             disp(['Average cruise lift coefficient is ' num2str(C_LavgCruise)])
            %             disp(['Induced drag coefficient during cruise is ' num2str(C_Di)])
            %             disp(['L/D during cruise is ' num2str(LiftToDrag)])
            %             disp(['SFC during cruise is ' num2str(SFC_35000)])
            %             disp(['Scheduled and calculated all out range differ by ' num2str(R_ao - range_total) ' nmi'])

            %% Section 7: Thrust Check
            C_LIC_check = (W_0/S)/(1481*delta*M_alt^2); % Recalculating initial cruise lift coefficient
            C_Di_IC = (C_LIC_check^2)/(pi*AR*e); % Recalculating initial cruise induced drag coefficient
            C_D_IC = C_D0 + C_Di_IC + 0.001; % Total drag coefficient at initial cruise
            LiftToDrag_IC = C_LIC_check/C_D_IC; % Initial cruise L/D ratio
            T_R_IC = (W_0/LiftToDrag_IC)/NumEngines; % Thrust required at initial cruise per engine
            T_R_JT_IC = T_R_IC*(T_SLS/T_e); % Thrust required per JT8D/9D engine at initial cruise
            if EngineType == 1
                % Thrust available at cruise for JT8D
                T_a_JT_cruise = -1568.2384575*(M_alt)^3 + 4892.0566740*(M_alt^2) - 3218.3519572*(M_alt) + 4032.2542516;
                % Thrust available at cruise for JT9D
            elseif EngineType == 2
                T_a_JT_cruise = -24987.9909806*(M_alt)^3 + 51566.7213979*(M_alt)^2 - 32354.7223933*(M_alt) + 15650.2395719;
            end
            if T_R_JT_IC > T_a_JT_cruise
                fprintf('NOT ENOUGH THRUST AT TOP OF CLIMB\n');
                ThrustAdjust = ThrustAdjust + 0.01;
            else
                fprintf('ENOUGH THRUST AT TOP OF CLIMB\n');
                %                 disp(['Thrust available at top of climb is ' num2str(T_a_JT_cruise) ' lbs'])
            end
        end

        %% Section 8: Gradients
        % Gradients - 1st T/O Segment
        C_LTO1 = C_LmaxTO/(1.2)^2; % Lift coefficient at 1st takeoff segment
        C_LTO1ratio = (C_LTO1/C_LmaxTO); % Ratio of lift coefficients at 1st takeoff segment
        % Incremental profile drag coefficient during 1st takeoff segment
        DeltaC_D0_TO1 = 0.1086760*(C_LTO1ratio)^3 - 0.0407543*(C_LTO1ratio)^2 - 0.0493280*(C_LTO1ratio) + 0.0321092;
        DeltaC_D_Gear = C_D0; % Incremental drag coefficient from landing gear
        C_D_TO1 = C_D0 + DeltaC_D0_TO1 + DeltaC_D_Gear + ((C_LTO1^2)/(pi*AR*e)); % Drag coefficient during 1st takeoff segment
        LifttoDrag_TO1 = C_LTO1/C_D_TO1; % L/D ratio during 1st takeoff segment
        T_R_TO1 = TakeoffWeight/LifttoDrag_TO1; % Thrust required during 1st takeoff segment
        % Sea-level Max Thrust/Dry Takeoff Thrust (JT8D/JT9D)
        if EngineType == 1
            T_a_TO1_JT = 7482.6479132*(M_LO)^2 - 8676.7374528*(M_LO) + 14819.6019504;
        elseif EngineType == 2
            T_a_TO1_JT = 1000*(37.4386807*(M_LO)^2 - 47.3982379*(M_LO) + 45.3366105);
        end
        T_a_TO1 = (T_e/T_SLS)*T_a_TO1_JT; % Thrust available per engine during 1st takeoff segment
        Gradient1 = ((((NumEngines - 1) * T_a_TO1) - T_R_TO1)/(TakeoffWeight)) * 100; % 1st takeoff segment gradient
        % 1st takeoff segment gradient check
        if Gradient1 > 0
            fprintf('Gradient 1 requirement met\n')
            %             disp(['Gradient 1 = ' num2str(Gradient1) '%'])
        else
            fprintf('Gradient 1 requirement not met\n');
        end

        % Gradients - 2nd T/O Segment
        C_D_TO2 = C_D_TO1 - DeltaC_D_Gear; % Drag coefficient during 2nd takeoff segment
        C_LTO2 = C_LTO1; % Lift coefficient during 2nd takeoff segment
        LifttoDrag_TO2 = C_LTO2/C_D_TO2; % L/D ratio during 2nd takeoff segment
        T_R_TO2 = TakeoffWeight/LifttoDrag_TO2; % Total thrust required during 2nd takeoff segment
        T_a_TO2 = T_a_TO1; % Thrust available per engine during 2nd takeoff segment
        Gradient2 = (((NumEngines - 1)*T_a_TO2 - T_R_TO2)/TakeoffWeight)*100; % 2nd takeoff segment gradient
        % 2nd takeoff segment gradient check
        if Gradient2 > 2.4
            fprintf('Gradient 2 requirement met\n')
            %             disp(['Gradient 2 = ' num2str(Gradient2) '%'])
        else
            fprintf('Gradient 2 requirement not met\n')
        end

        % Gradients - 3rd T/O Segment
        % Clean Wing Configuration Graph Selection
        % At sweep between 0 and 15
        if (sweep >= 0 && sweep < 15)
            C_LMaxClean = (-301.3947698*(tc_ave)^3 + 77.2839301*(tc_ave)^2 - 1.9494582*(tc_ave) + 0.9322080)*((15 - sweep)/(15 - 0)) + (-338.4238696*(tc_ave)^3 + 89.1438114*(tc_ave)^2 - 3.0248621*(tc_ave) + 0.9078337)*(1 - ((15 - sweep)/(15 - 0)));
            % At sweep between 15 and 35
        elseif (sweep >= 15 && sweep < 35)
            C_LMaxClean = (-338.4238696*(tc_ave)^3 + 89.1438114*(tc_ave)^2 - 3.0248621*(tc_ave) + 0.9078337)*((35 - sweep)/(35 - 15)) + (-319.7572509*(tc_ave)^3 + 84.3328349*(tc_ave)^2 - 2.7221484*(tc_ave) + 0.8542277)*(1 - ((35 - sweep)/(35 - 15)));
            % At sweep of 35
        elseif sweep == 35
            C_LMaxClean = -319.7572509*(tc_ave)^3 + 84.3328349*(tc_ave)^2 - 2.7221484*(tc_ave) + 0.8542277;
        end
        rho_1000 = 2.3081e-3; % Density at 1000 ft
        sigma_1000 = rho_1000/rho_SL84; % Density ratio at 1000 ft on a hot day
        V_TO3 = 1.2*sqrt((296*WS_TO)/(sigma_1000*C_LMaxClean)); % 3rd takeoff segment velocity in kts
        sound_1000 = sqrt(gamma*R*515.12)/1.688; % Speed of sound in kts
        M_TO3 = V_TO3/sound_1000; % Mach number during 3rd takeoff segment
        C_LTO3 = C_LMaxClean/(1.2^2); % 3rd takeoff segment lift coefficient
        C_D_TO3 = C_D0 + (C_LTO3^2/(pi*AR*e)); % 3rd takeoff segment drag coefficient
        LifttoDrag_TO3 = C_LTO3/C_D_TO3; % 3rd takeoff segment L/D ratio
        T_R_TO3 = TakeoffWeight/LifttoDrag_TO3; % Total thrust required during 3rd takeoff segment
        % Max Continuous/Max Climb Thrust (JT8D/JT9D)
        if EngineType == 1
            T_a_JT_TO3 = -4416.7762017*(M_TO3)^3 + 11231.9430229*(M_TO3)^2 - 9754.6007222*(M_TO3) + 12663.6807136;
        elseif EngineType == 2
            T_a_JT_TO3 = -20500.7401839*(M_TO3)^3 + 42302.4847473*(M_TO3)^2 - 43323.5787736*(M_TO3) + 37961.2002238;
        end
        T_a_TO3 = (T_e/T_SLS)*T_a_JT_TO3; % Thrust available per engine during 3rd takeoff segment
        Gradient3 = (((NumEngines - 1)*T_a_TO3 - T_R_TO3)/TakeoffWeight)*100; % 3rd segment gradient
        % 3rd takeoff segment gradient check
        if Gradient3 > 1.2
            fprintf('Gradient 3 requirement met\n')
            %             disp(['Gradient 3 = ' num2str(Gradient3) '%'])
        else
            fprintf('Gradient 3 requirement not met\n')
        end

        % Approach
        C_LAp = C_LmaxTO/(1.3)^2; % Approach lift coefficient
        C_LApRatio = C_LAp/C_LmaxTO; % Approach lift coefficient ratio
        % Incremental profile drag during approach
        DeltaC_D0Ap = 0.1086760*(C_LApRatio)^3 - 0.0407543*(C_LApRatio)^2 - 0.0493280*(C_LApRatio) + 0.0321092;
        C_DAp = C_D0 + DeltaC_D0Ap + (C_LAp^2 /(pi*AR*e)); % Approach drag coefficient
        LifttoDrag_Ap = C_LAp/C_DAp; % Apprach L/D ratio
        LandingWeight = WS_LD*S; % Landing weight
        T_R_Ap = LandingWeight/LifttoDrag_Ap; % Total thrust required during approach
        V_Ap = sqrt((296*WS_LD)/(sigma*C_LAp)); % Approach speed in kts
        M_Ap = V_Ap/sound; % Approach Mach number
        % Max Continuous/Max Climb Thrust (JT8D/JT9D)
        if EngineType == 1
            T_a_JT_Ap = -4416.7762017*(M_Ap)^3 + 11231.9430229*(M_Ap)^2 - 9754.6007222*(M_Ap) + 12663.6807136;
        elseif EngineType == 2
            T_a_JT_Ap = -20500.7401839*(M_Ap)^3 + 42302.4847473*(M_Ap)^2 - 43323.5787736*(M_Ap) + 37961.2002238;
        end
        T_a_Ap = (T_e/T_SLS)*T_a_JT_Ap; % Thrust available per engine during approach
        GradientAp = (((NumEngines - 1)*T_a_Ap - T_R_Ap)/LandingWeight)*100; % Approach gradient
        % Approach gradient check
        if GradientAp > 2.1
            fprintf('Approach gradient requirement met\n')
            %             disp(['Approach gradient = ' num2str(GradientAp) '%'])
        else
            fprintf('Approach gradient requirement not met\n')
        end

        % Landing
        C_L_LD = C_LmaxLD/(1.3^2);
        C_L_LDRatio = C_L_LD/C_LmaxLD; % Landing lift coefficient ratio
        % Incremental profile drag coefficient during landing
        DeltaC_D0LD = 0.0792529*(C_L_LDRatio)^3 + 0.0084234*(C_L_LDRatio)^2 - 0.0690409*(C_L_LDRatio) + 0.0413224;
        C_DLD = C_D0 + DeltaC_D0LD + DeltaC_D_Gear + (C_L_LD^2 /(pi*AR*e)); % Landing drag coefficient
        LifttoDrag_LD = C_L_LD/C_DLD; % Landing L/D ratio
        T_R_LD = LandingWeight/LifttoDrag_LD; % Total thrust required during landing
        V_LD = Approach; % Landing approach speed in kts
        sound_SL = (sqrt(gamma*R*543.67))/1.688; % Sea-level speed of sound in kts
        M_LD = V_LD/(sound_SL); % Landing Mach number
        % Sea-level Max Thrust/Dry Takeoff Thrust (JT8D/JT9D)
        if EngineType == 1
            T_a_JT_LD = 7482.6479132*(M_LD)^2 - 8676.7374528*(M_LD) + 14819.6019504;
        elseif EngineType == 2
            T_a_JT_LD = 1000*(37.4386807*(M_LD)^2 - 47.3982379*(M_LD) + 45.3366105);
        end
        T_a_LD = (T_e/T_SLS)*T_a_JT_LD; % Thrust available per engine during landing
        Gradient_LD = (((NumEngines*T_a_LD) - T_R_LD)/(LandingWeight))*100; % Landing gradient
        % Landing gradient check
        if Gradient_LD > 3.2
            fprintf('Landing gradient requirement met\n')
            %             disp(['Landing gradient = ' num2str(Gradient_LD) '%'])
        else
            fprintf('Landing gradient not met\n')
        end

        %% Section 9: Direct Operating Cost
        % Block speed parameters
        D = range*1.15; % CAB trip distance in statute miles
        % All times are in hours
        t_gm = 0.25; % Ground maneuver time
        t_cl = 0.18; % Time to climb
        t_d = 0; % Time to descend
        t_am = 0.1; % Time for air maneuver
        K_a = 0.02*D; % Airway distance increment
        D_cl = 83; % Distance to climb
        D_d = 0;% Distance to descend
        V_cr = V_cruise; % Average true airpseed in cruise in kts
        t_cr = ((D + K_a + 20) - (D_cl + D_d))/V_cr; % Time at cruise altitude
        % Block speed in mph
        V_b = D/(t_gm + t_cl + t_d + t_cr + t_am);
        t_b = t_gm + t_cl + t_d + t_cr + t_am; % Block time

        % Block fuel parameters
        % Fuel is measured in lbs
        F_gm = 0; % Ground maneuver fuel for taxi and takeoff
        F_cl = 9624; % Fuel to climb to cruise altitude
        F_d = 0; % Fuel required to descend
        F_crAndF_am = T_R_cruise*SFC_35000*(t_cr + t_am); % Combined summation fuel term
        % Block fuel in lbs
        F_b = F_gm + F_cl + F_d + F_crAndF_am;

        % Flying Operations Cost
        % a. Flight crew
        TotalPayload = W_PL/2000; % In tons
        V_c = V_cruise*1.15; % Cruise speed in statute mph
        CostPerBlockHr = 17.849*(V_c*(TakeoffWeight/1e5))^0.3 + 40.83; % [$/block hr]
        C_TMCrew = (CostPerBlockHr)/(V_b*TotalPayload); % [$/ton-mile] for crew

        % b. Fuel and Oil
        % All fuel costs are in $/lb
        C_fuel = 0.28*(1/6.4); % Fuel cost for domestic flights
        C_ot = 2.15; % Cost of oil turbines
        C_TMFuel = 1.02*((F_b*C_fuel + NumEngines*C_ot*t_b*0.135)/(D*TotalPayload)); % [$/ton-mile] fuel

        % c. Hull Insurance
        % All costs use 1975 dollars
        W_FC = 0; % Flight crew weight
        W_a = TakeoffWeight*(1 - W_f) - W_PL - W_FC - W_PP*TakeoffWeight; % Net airplane weight
        C_a = 2.4e6 + 87.5*W_a; % Airplane cost
        C_e = 590e3 + 16*T_e; % Engine cost
        C_t = NumEngines*C_e + C_a; % Total airplane cost
        U = 630 + (4000/(1 + (1/(t_b + 0.5)))); % Utilization [block hrs/year]
        IRA = 0.01; % Insurance rate per $ value
        C_TMIns = (IRA*C_t)/(U*V_b*TotalPayload); % [$/ton-mile] insurance

        % Direct Maintenance
        % a. Airframe labor
        K_FHa = 4.9169*log10(W_a/1e3) - 6.425; % Direct labor [MH/FH]
        K_FCa = 0.21256*(log10(W_a/1e3))^3.7375; % Direct labor [MH/cycle]
        t_f = t_b - t_gm; % Flight time
        R_L = 8.60; % Labor rate in 1975 [$/hr]
        C_TMAirLabor = ((K_FHa*t_f + K_FCa)/(V_b*t_b*TotalPayload))*R_L; % [$/ton-mile] labor

        % b. Airframe material
        C_FHa = 1.5994*(C_a/1e6) + 3.2463; % Material cost [$/FH]
        C_FCa = 1.9229*(C_a/1e6) + 2.2504; % Material cost [$/cycle]
        C_TMAirMat = (C_FHa*t_f + C_FCa)/(V_b*t_b*TotalPayload); % [$/ton-mile] airframe materials

        % c. Engine labor
        K_FHe = (NumEngines*(T_e/1e3))/(0.82715*(T_e/1e3) + 13.639); % Direct labor [MH/FH]
        K_FCe = 0.20*NumEngines; % Direct labor [MH/cycle]
        C_TMEngLab = ((K_FHe*t_f + K_FCe)/(V_b*t_b*TotalPayload))*R_L; % [$/ton-mile] engine labor

        % d. Engine material
        C_FHe = (28.2353*(C_e/1e6) - 6.5176)*NumEngines; % Material cost [$/FH]
        C_FCe = (3.6698*(C_e/1e6) + 1.3685)*NumEngines; % Material cost [$/cycle]
        C_TMEngMat = (C_FHe*t_f + C_FCe)/(V_b*t_b*TotalPayload); % [$/ton-mile] engine materials

        % e. Total Maintenance - Burdened
        C_TMBur = (C_TMAirLabor + C_TMAirMat + C_TMEngLab + C_TMEngMat)*2; % [$/ton-mile] total

        % Depreciation
        D_a = 14; % Depreciation period [yrs]
        C_TMDep = ((C_TMBur + 0.06*(C_TMBur -NumEngines*C_e) + 0.3*NumEngines*C_e)/(D_a*U*V_b*TotalPayload)); % [$/ton-mile] depreciated

        % Total Direct Operating Cost
        % [$/ton-mile]
        DOCTM = C_TMCrew + C_TMFuel + C_TMIns + C_TMAirLabor + C_TMAirMat + C_TMEngLab + C_TMEngMat + C_TMBur + C_TMDep;
        % [%/pax-mile]
        DOCPM = DOCTM*(TotalPayload/P);
        disp(['Direct Operating Cost per ton-mile = $' num2str(DOCTM) '/ton-mile'])
        disp(['Direct Operating Cost per pax-mile = $' num2str(DOCPM) '/PAX-mile'])
        % Storing values in an array for the inner for loop
        DOC_Final(i,j) = DOCTM;
        TakeoffWeightFinal(i,j) = TakeoffWeight;
    end
end

%% Section 10: Plots
figure(1)
plot(ARRange,DOC_Final)
hold on

figure(2)
plot(ARRange,TakeoffWeightFinal)
hold on

Answer 1

[Torsten]

You define a variable "length" as

length = (3.76*(P/abreast) + 33.2); % Fuselage length in ft

This conflicts with the MATLAB function "length" to determine array sizes.

Rename the variable in your code.

I don't know if this will make your code work, but it's at least a problem you should correct.

[Kenneth Collio]

Thank you so much! That fixed my code.