Loop to create matrix with 2 inputs and 2 outputs

Hello,

I recently had an assignment to create a code and got it right, but afterwards it was recommended to me that it would have been easier to do it a different way. The way the code below is set up was to be more like a calculator. You input the 2 values and it displays 2 outputs. This was very tedious to do and a loop would have been much quicker.

For my situation I have 2 inputs and 2 outputs that are a function of the 2 inputs. I want to have the inputs change through a loop and save the output variables to a matrix.

The value of TINT needs to be set to range from 1100 to 1800 at intervals of 100 while the OPR needs to be 4 to 30 increments of 2. I have a basic understanding of how to do this if it was only one variable changing. I am not sure how to do 2 variables.

%------------------------------------------------------
%Given Values                                         |
Ra = 0.287;   %kJ/kg/k                                |
ga = 1.4;     %Gamma Air                              |
gg = 1.3333;   %Gamma gas                              |
Nd = 0.93;    %Efficiency of Diffuser                 |
Nc = 0.87;    %Efficiency of Compressor               |
Nb = 0.98;    %Efficiency of Combustion Chamber       |
Nm = 0.99;    %Efficiency of Mechanical               |
Nt = 0.90;    %Efficiency of Turbines                 |
Nn = 0.95;    %Efficiency of Nozzle                   |
Ma = 0.84;    %Mach                                   |
Alt = 5;      %km                                     |
%OPR = 9;      %Overal Pressure Ratio                 |   
prompt = 'Overall Pressure Ratio? ';  %               |
OPR = input(prompt);                  %               |
PRClp = sqrt(OPR); %Pressure Ratio over LP Compressor |
PRChp = sqrt(OPR); %Pressure Ratio over HP Compressor |
deltap0 = 0.04;   % 4% of delievery a loss            |
%-----------------------------------------------------|

%From Handout 5 based on Altitude of 5km--------------|
Ta = 255.7;    %K                                     |
Pa = 54.05;    %kPa                                   |
rho = 0.736;  %kg/m^3                                 |
u = 1.63;     %N*s/m^2                                |
%-----------------------------------------------------|

Va = Ma * sqrt(ga*Ra*Ta*1000);   %in m/s
Cpoa = (ga*Ra)/(ga-1);          %in kJ/kg/K

T01 = Ta*(1 + ((ga-1)/2)*Ma^2);
P01 = Pa*(1 + Nd*((ga-1)/2)*(Ma^2))^(ga/(ga-1));
P02 = P01*PRClp;
T02 = T01 + ((T01*((PRClp)^((ga-1)/ga)-1)/Nc));
Wclp = Cpoa*(T01-T02);
P03 = P02*PRChp;
T03 = T02 + ((T02*((PRChp)^((ga-1)/ga)-1)/Nc));
Wchp = -Cpoa*(T03-T02);
Tr = T03;
   % Tp
    prompt = 'Temperature of Products? ';
    Tp = input(prompt);
    T04 = Tp;

      if (Tp <= 1600 && Tp>=400)
              COAp = 299180; COBp = 37.85; COCp = -4571.9; CO2Ap = 56835; 
              CO2Bp = 66.27; CO2Cp = -11634; H2OAp = 88923; H2OBp = 49.36; 
              H2OCp = -7940.8; N2Ap = 31317; N2Bp = 37.46; N2Cp = -4559.3;
              O2Ap = 43388; O2Bp = 42.27; O2Cp = -6635.4;

          %Tp > 1600K
          elseif (Tp > 1600);
              COAp = 309070; COBp = 39.29; COCp = -6201; CO2Ap = 93048; 
              CO2Bp = 68.58; CO2Cp = -16979; H2OAp = 154670; H2OBp = 60.43; 
              H2OCp = -19212; N2Ap = 44639; N2Bp = 39.32; N2Cp = -6753.4; 
              O2Ap = 127010; O2Bp = 46.25; O2Cp = -18798;
      end

MC = 14.4;
MH = 24.9;
MO = 0;
Hrp = -8561991.6;
Ycc = MC + (MH/4) + (MO/2);
Ymin = Ycc - (MC/2); 
Mfuel = MC*12 + MH + MO*16;
Mair = 28.97;

%Find f
    %based on Tr < 1600K
    COAr = 299180; COBr = 37.85; COCr = -4571.9; CO2Ar = 56835;
    CO2Br = 66.27; CO2Cr = -11634; H2OAr = 88923; H2OBr = 49.36;
    H2OCr = -7940.8; N2Ar = 31317; N2Br = 37.46; N2Cr = -4559.3;
    O2Ar = 43388; O2Br = 42.27; O2Cr = -6635.4;

      %based on Excess Air
      N1 = 0;
      N2 = MC;
      N3 = (MH/2); 

      HbarCO = 0;
      HbarCO2 = (CO2Ap + CO2Bp*Tp + CO2Cp*log(Tp)) - (CO2Ar + CO2Br*Tr + CO2Cr*log(Tr));
      HbarH2O = (H2OAp + H2OBp*Tp + H2OCp*log(Tp)) - (H2OAr + H2OBr*Tr + H2OCr*log(Tr));
      HbarO2 = (O2Ap + O2Bp*Tp + O2Cp*log(Tp)) - (O2Ar + O2Br*Tr + O2Cr*log(Tr));
      HbarN2 = (N2Ap + N2Bp*Tp + N2Cp*log(Tp)) - (N2Ar + N2Br*Tr + N2Cr*log(Tr));

      syms Y
      eqn = Hrp + (N1 * 282800) + (N2 * HbarCO2) + (N3 * HbarH2O) + ((Y - Ycc)* HbarO2) + ((3.76 * Y)* HbarN2)==0;

      Y = solve(eqn,Y);

      %Finding f   
      fideal= ((1/(4.76*Y)))*((Mfuel)/(Mair));
      f = fideal/0.98;  %Getting f actual from ideal and efficiency

     P04 = (1- deltap0)*P03;
     Cp0g = (gg*Ra)/(gg-1);      %in kJ/kg/K
     T05 = T04 +((Wchp)/(Nm*(1+0.0168170)*Cp0g));  %K
     P05 = (P04)*(1-((1-(T05/T04))/Nt))^(gg/(gg-1));  %kPa
     T06 = T05 + (Wclp/(Nm*(1+f)*Cp0g));
     %double(T06);
     P06 = P05*(1-((1-(T06/T05))/Nt))^(gg/(gg-1));  %Page 26!!!!
    double(P06);
    %Return to the particular converging nozzle of the turbojet problem
    P0i = P06;
    T0i = T06;
    T0e = T06;
    PstarovP06 = (1-(1/Nn)*(1-(2/(gg+1))))^(gg/(gg-1));
    PaovP06 = Pa/P06;

    %Choke Test
    if(PstarovP06 >= PaovP06)
        choked = 1;
    elseif (PstarovP06 > PaovP06)
        choked =2;
    end

    %If Choked
    if(choked == 1)
        Me = 1;
        Pe = P06*(PstarovP06);
        Te = T06*(2/(gg+1));
        rhoe = Pe/(Ra*Te);   %kg/m^3
        Ve = sqrt(gg*Ra*1000*Te);
        double(Ve);   

    %If not choked
    elseif (choked == 2)
        Pe = Pa;  
        Te = T06*(1-(Nn(1-(Pe/P06)^((gg-1)/gg))));
        rhoe = Pe/(Ra*Te);   %kg/m^3
        Me = sqrt(((T06/Te)-1)*(2/(gg-1)));   %Me<1 if not choked
        Ve = Me*sqrt(gg*Ra*1000*Te);
    end

%Compute Fs (Specific Thrust)
Fs = ((1+f)*Ve)-Va +(Pe*1000-Pa*1000)*((1+f)/(rhoe*Ve));  %N/(kg/sec)
double(Fs)
%Compute tsfc
tsfc = f/Fs*3600;
double(tsfc)             %kgfuel/hour/N

Below is the way that I tried to edit the code to accomplish this and it has been running for about 5 minutes with no results.

clc clear format long g

%------------------------------------------------------
%Given Values                                         |
Ra = 0.287;   %kJ/kg/k                                |
ga = 1.4;     %Gamma Air                              |
gg = 1.3333;   %Gamma gas                              |
Nd = 0.93;    %Efficiency of Diffuser                 |
Nc = 0.87;    %Efficiency of Compressor               |
Nb = 0.98;    %Efficiency of Combustion Chamber       |
Nm = 0.99;    %Efficiency of Mechanical               |
Nt = 0.90;    %Efficiency of Turbines                 |
Nn = 0.95;    %Efficiency of Nozzle                   |
Ma = 0.84;    %Mach                                   |
Alt = 5;      %km                                     |
%OPR = 9;      %Overal Pressure Ratio                 |   
%prompt = 'Overall Pressure Ratio? ';  %               |
%OPR = input(prompt);                  %               |
%PRClp = sqrt(OPR); %Pressure Ratio over LP Compressor |
%PRChp = sqrt(OPR); %Pressure Ratio over HP Compressor |
deltap0 = 0.04;   % 4% of delievery a loss            |
%-----------------------------------------------------|

%From Handout 5 based on Altitude of 5km--------------|
Ta = 255.7;    %K                                     |
Pa = 54.05;    %kPa                                   |
rho = 0.736;  %kg/m^3                                 |
u = 1.63;     %N*s/m^2                                |
%-----------------------------------------------------|

     % Tp
      %prompt = 'Temperature of Products? ';
      %Tp = input(prompt);
      %T04 = Tp;
     for TINT = 1000:1700
         TINT = TINT + 100;
         for OPR = 2:28
             OPR = OPR + 2;
             PRClp = sqrt(OPR); %Pressure Ratio over LP Compressor 
             PRChp = sqrt(OPR); %Pressure Ratio over HP Compressor
              Va = Ma * sqrt(ga*Ra*Ta*1000);   %in m/s
              Cpoa = (ga*Ra)/(ga-1);          %in kJ/kg/K

              T01 = Ta*(1 + ((ga-1)/2)*Ma^2);
              P01 = Pa*(1 + Nd*((ga-1)/2)*(Ma^2))^(ga/(ga-1));
              P02 = P01*PRClp;
              T02 = T01 + ((T01*((PRClp)^((ga-1)/ga)-1)/Nc));
              Wclp = Cpoa*(T01-T02);
              P03 = P02*PRChp;
              T03 = T02 + ((T02*((PRChp)^((ga-1)/ga)-1)/Nc));
              Wchp = -Cpoa*(T03-T02);
              Tr = T03;
              Tp = TINT;
              T04 = Tp;

      if (Tp <= 1600 && Tp>=400)
              COAp = 299180; COBp = 37.85; COCp = -4571.9; CO2Ap = 56835; 
              CO2Bp = 66.27; CO2Cp = -11634; H2OAp = 88923; H2OBp = 49.36; 
              H2OCp = -7940.8; N2Ap = 31317; N2Bp = 37.46; N2Cp = -4559.3;
              O2Ap = 43388; O2Bp = 42.27; O2Cp = -6635.4;

          %Tp > 1600K
          elseif (Tp > 1600);
              COAp = 309070; COBp = 39.29; COCp = -6201; CO2Ap = 93048; 
              CO2Bp = 68.58; CO2Cp = -16979; H2OAp = 154670; H2OBp = 60.43; 
              H2OCp = -19212; N2Ap = 44639; N2Bp = 39.32; N2Cp = -6753.4; 
              O2Ap = 127010; O2Bp = 46.25; O2Cp = -18798;
      end

MC = 14.4;
MH = 24.9;
MO = 0;
Hrp = -8561991.6;
Ycc = MC + (MH/4) + (MO/2);
Ymin = Ycc - (MC/2); 
Mfuel = MC*12 + MH + MO*16;
Mair = 28.97;

%Find f
    %based on Tr < 1600K
    COAr = 299180; COBr = 37.85; COCr = -4571.9; CO2Ar = 56835;
    CO2Br = 66.27; CO2Cr = -11634; H2OAr = 88923; H2OBr = 49.36;
    H2OCr = -7940.8; N2Ar = 31317; N2Br = 37.46; N2Cr = -4559.3;
    O2Ar = 43388; O2Br = 42.27; O2Cr = -6635.4;

      %based on Excess Air
      N1 = 0;
      N2 = MC;
      N3 = (MH/2); 

      HbarCO = 0;
      HbarCO2 = (CO2Ap + CO2Bp*Tp + CO2Cp*log(Tp)) - (CO2Ar + CO2Br*Tr + CO2Cr*log(Tr));
      HbarH2O = (H2OAp + H2OBp*Tp + H2OCp*log(Tp)) - (H2OAr + H2OBr*Tr + H2OCr*log(Tr));
      HbarO2 = (O2Ap + O2Bp*Tp + O2Cp*log(Tp)) - (O2Ar + O2Br*Tr + O2Cr*log(Tr));
      HbarN2 = (N2Ap + N2Bp*Tp + N2Cp*log(Tp)) - (N2Ar + N2Br*Tr + N2Cr*log(Tr));

      syms Y
      eqn = Hrp + (N1 * 282800) + (N2 * HbarCO2) + (N3 * HbarH2O) + ((Y - Ycc)* HbarO2) + ((3.76 * Y)* HbarN2)==0;

      Y = solve(eqn,Y);

      %Finding f   
      fideal= ((1/(4.76*Y)))*((Mfuel)/(Mair));
      f = fideal/0.98;  %Getting f actual from ideal and efficiency

     P04 = (1- deltap0)*P03;
     Cp0g = (gg*Ra)/(gg-1);      %in kJ/kg/K
     T05 = T04 +((Wchp)/(Nm*(1+0.0168170)*Cp0g));  %K
     P05 = (P04)*(1-((1-(T05/T04))/Nt))^(gg/(gg-1));  %kPa
     T06 = T05 + (Wclp/(Nm*(1+f)*Cp0g));
     %double(T06);
     P06 = P05*(1-((1-(T06/T05))/Nt))^(gg/(gg-1));  %Page 26!!!!
    double(P06);
    %Return to the particular converging nozzle of the turbojet problem
    P0i = P06;
    T0i = T06;
    T0e = T06;
    PstarovP06 = (1-(1/Nn)*(1-(2/(gg+1))))^(gg/(gg-1));
    PaovP06 = Pa/P06;

    %Choke Test
    if(PstarovP06 >= PaovP06)
        choked = 1;
    elseif (PstarovP06 > PaovP06)
        choked =2;
    end

    %If Choked
    if(choked == 1)
        Me = 1;
        Pe = P06*(PstarovP06);
        Te = T06*(2/(gg+1));
        rhoe = Pe/(Ra*Te);   %kg/m^3
        Ve = sqrt(gg*Ra*1000*Te);
        double(Ve);   

    %If not choked
    elseif (choked == 2)
        Pe = Pa;  
        Te = T06*(1-(Nn(1-(Pe/P06)^((gg-1)/gg))));
        rhoe = Pe/(Ra*Te);   %kg/m^3
        Me = sqrt(((T06/Te)-1)*(2/(gg-1)));   %Me<1 if not choked
        Ve = Me*sqrt(gg*Ra*1000*Te);
    end

%Compute Fs (Specific Thrust)
Fs = ((1+f)*Ve)-Va +(Pe*1000-Pa*1000)*((1+f)/(rhoe*Ve));  %N/(kg/sec)
double(Fs);

%Compute tsfc
tsfc = f/Fs*3600;
double(tsfc);            %kgfuel/hour/N
Fs(TINT) = Fs;
       end
   end
   Fs

Thank you for any help you may be able to provide,

Colby

Answer 1

[Geoff Hayes]

Colby - you almost have it. Your code for the two for loops is as follows

     for TINT = 1000:1700
         TINT = TINT + 100;
         for OPR = 2:28
             OPR = OPR + 2;

Since TINT and OPR are your indexing variables, you don't need to increment either through code. This will happen automatically so long as you set your step value correctly in the for loop. The above can be rewritten as

     for TINT = 1100:100:1800
         for OPR = 4:2:30
             % etc.
         end
     end

So on the first iteration of the outer for loop, TINT will be assigned a value on 1100, and subsequent iterations of this loop will assign values of 1200, 1300, 1400,...,1800 to TINT.

The inner for loop will initialize *OPR to 4 with subsequent iterations assigning it values of 6,8,10,...,30.

[colbyjh]

Thank you very much for the quick answer! My only other question is how to save it to a vector or matrix so that I can see all of the values for Fs and tsfc. If it was just one variable changing and 1 loop, i think what is posted below would work:

Fs(TINT) = Fs;

If I do it this way, then it would get over written once the 2nd variable changes since there will be sever iterations with the same TINT, but they will have different OPR

[Geoff Hayes]

You probably don't want to use

 Fs(TINT)

since TINT is at least 1100 and so your output matrix Fs will be larger than it should be. Also, by doing the above, you have created a matrix named Fs which is the same name as your local variable (scalar) Fs. So on each iteration of your inner for loop, you will be overwriting the matrix with the scalar. (In the end, once all iterating has completed, you would have a matrix but it would only have one non-zero value which would be the last Fs).

You will want to create an mxn array where m is the number of TINT values and n is the number of OPR values (so an 8x14 matrix) as

 FsArray   = zeros(8,14);
 TsfcArray = zeros(8,14);

On each iteration of your inner for loop, you will need to update the appropriate element in each array. Something like

 FsArray(u,v)   = Fs;
 TsfcArray(u,v) = tsfc;

The trick then is to keep track of u and v as they range from 1 to 8 and 1 to 14 respectively.

[colbyjh]

Thank you so much!! I had been trying to figure this out for so long and now it is working as I had hoped!

Below is the finished code in case in benefits anyone else.

clc
clear
format long g

%------------------------------------------------------
%Given Values                                         |
Ra = 0.287;   %kJ/kg/k                                |
ga = 1.4;     %Gamma Air                              |
gg = 1.3333;   %Gamma gas                              |
Nd = 0.93;    %Efficiency of Diffuser                 |
Nc = 0.87;    %Efficiency of Compressor               |
Nb = 0.98;    %Efficiency of Combustion Chamber       |
Nm = 0.99;    %Efficiency of Mechanical               |
Nt = 0.90;    %Efficiency of Turbines                 |
Nn = 0.95;    %Efficiency of Nozzle                   |
Ma = 0.84;    %Mach                                   |
Alt = 5;      %km                                     |
deltap0 = 0.04;   % 4% of delievery a loss            |
%-----------------------------------------------------|

%From Handout 5 based on Altitude of 5km--------------|
Ta = 255.7;    %K                                     |
Pa = 54.05;    %kPa                                   |
rho = 0.736;  %kg/m^3                                 |
mu = 1.63;     %N*s/m^2                                |
%-----------------------------------------------------|

%Create array to accept all answers for Fs and TINT---|
FsArray = zeros(14,8);                               %|
TsfcArray = zeros(14,8);                             %|
%-----------------------------------------------------|

       u = 1;     %location of first OPR
       v = 1;  %locatino of first TINT
     for TINT = 1100:100:1800
         u = 1;
         for OPR = 4:2:30
             PRClp = sqrt(OPR); %Pressure Ratio over LP Compressor 
             PRChp = sqrt(OPR); %Pressure Ratio over HP Compressor
              Va = Ma * sqrt(ga*Ra*Ta*1000);   %in m/s
              Cpoa = (ga*Ra)/(ga-1);          %in kJ/kg/K

              T01 = Ta*(1 + ((ga-1)/2)*Ma^2);
              P01 = Pa*(1 + Nd*((ga-1)/2)*(Ma^2))^(ga/(ga-1));
              P02 = P01*PRClp;
              T02 = T01 + ((T01*((PRClp)^((ga-1)/ga)-1)/Nc));
              Wclp = Cpoa*(T01-T02);
              P03 = P02*PRChp;
              T03 = T02 + ((T02*((PRChp)^((ga-1)/ga)-1)/Nc));
              Wchp = -Cpoa*(T03-T02);
              Tr = T03;
              Tp = TINT;
              T04 = Tp;

      if (Tp <= 1600 && Tp>=400)
              COAp = 299180; COBp = 37.85; COCp = -4571.9; CO2Ap = 56835; 
              CO2Bp = 66.27; CO2Cp = -11634; H2OAp = 88923; H2OBp = 49.36; 
              H2OCp = -7940.8; N2Ap = 31317; N2Bp = 37.46; N2Cp = -4559.3;
              O2Ap = 43388; O2Bp = 42.27; O2Cp = -6635.4;

          %Tp > 1600K
          elseif (Tp > 1600);
              COAp = 309070; COBp = 39.29; COCp = -6201; CO2Ap = 93048; 
              CO2Bp = 68.58; CO2Cp = -16979; H2OAp = 154670; H2OBp = 60.43; 
              H2OCp = -19212; N2Ap = 44639; N2Bp = 39.32; N2Cp = -6753.4; 
              O2Ap = 127010; O2Bp = 46.25; O2Cp = -18798;
      end

MC = 14.4;
MH = 24.9;
MO = 0;
Hrp = -8561991.6;
Ycc = MC + (MH/4) + (MO/2);
Ymin = Ycc - (MC/2); 
Mfuel = MC*12 + MH + MO*16;
Mair = 28.97;

%Find f
    %based on Tr < 1600K
    COAr = 299180; COBr = 37.85; COCr = -4571.9; CO2Ar = 56835;
    CO2Br = 66.27; CO2Cr = -11634; H2OAr = 88923; H2OBr = 49.36;
    H2OCr = -7940.8; N2Ar = 31317; N2Br = 37.46; N2Cr = -4559.3;
    O2Ar = 43388; O2Br = 42.27; O2Cr = -6635.4;

      %based on Excess Air
      N1 = 0;
      N2 = MC;
      N3 = (MH/2); 

      HbarCO = 0;
      HbarCO2 = (CO2Ap + CO2Bp*Tp + CO2Cp*log(Tp)) - (CO2Ar + CO2Br*Tr + CO2Cr*log(Tr));
      HbarH2O = (H2OAp + H2OBp*Tp + H2OCp*log(Tp)) - (H2OAr + H2OBr*Tr + H2OCr*log(Tr));
      HbarO2 = (O2Ap + O2Bp*Tp + O2Cp*log(Tp)) - (O2Ar + O2Br*Tr + O2Cr*log(Tr));
      HbarN2 = (N2Ap + N2Bp*Tp + N2Cp*log(Tp)) - (N2Ar + N2Br*Tr + N2Cr*log(Tr));

      syms Y
      eqn = Hrp + (N1 * 282800) + (N2 * HbarCO2) + (N3 * HbarH2O) + ((Y - Ycc)* HbarO2) + ((3.76 * Y)* HbarN2)==0;

      Y = solve(eqn,Y);

      %Finding f   
      fideal= ((1/(4.76*Y)))*((Mfuel)/(Mair));
      f = fideal/0.98;  %Getting f actual from ideal and efficiency

     P04 = (1- deltap0)*P03;
     Cp0g = (gg*Ra)/(gg-1);      %in kJ/kg/K
     T05 = T04 +((Wchp)/(Nm*(1+0.0168170)*Cp0g));  %K
     P05 = (P04)*(1-((1-(T05/T04))/Nt))^(gg/(gg-1));  %kPa
     T06 = T05 + (Wclp/(Nm*(1+f)*Cp0g));
     %double(T06);
     P06 = P05*(1-((1-(T06/T05))/Nt))^(gg/(gg-1));  %Page 26!!!!
    double(P06);
    %Return to the particular converging nozzle of the turbojet problem
    P0i = P06;
    T0i = T06;
    T0e = T06;
    PstarovP06 = (1-(1/Nn)*(1-(2/(gg+1))))^(gg/(gg-1));
    PaovP06 = Pa/P06;

    %Choke Test
    if(PstarovP06 >= PaovP06)
        choked = 1;
    elseif (PstarovP06 > PaovP06)
        choked =2;
    end

    %If Choked
    if(choked == 1)
        Me = 1;
        Pe = P06*(PstarovP06);
        Te = T06*(2/(gg+1));
        rhoe = Pe/(Ra*Te);   %kg/m^3
        Ve = sqrt(gg*Ra*1000*Te);
        double(Ve);   

    %If not choked
    elseif (choked == 2)
        Pe = Pa;  
        Te = T06*(1-(Nn(1-(Pe/P06)^((gg-1)/gg))));
        rhoe = Pe/(Ra*Te);   %kg/m^3
        Me = sqrt(((T06/Te)-1)*(2/(gg-1)));   %Me<1 if not choked
        Ve = Me*sqrt(gg*Ra*1000*Te);
    end

%Compute Fs (Specific Thrust)
Fs = ((1+f)*Ve)-Va +(Pe*1000-Pa*1000)*((1+f)/(rhoe*Ve));  %N/(kg/sec)
double(Fs);

%Compute tsfc
tsfc = f/Fs*3600;
double(tsfc);            %kgfuel/hour/N

%Add values to the array------|
FsArray(u,v) = Fs;           %|
TsfcArray(u,v) = tsfc;       %|
%-----------------------------|

u = u+1;  %update u array element counter                  
       end
   v = v + 1;  % update v array element counter
   end

FsArray
TsfcArray

Thank you,

Colby