Need help debugging given matlab code for school assignment. Recieving weird error need help moving past

Here is the original code which I have half debugged, the professor left areas with the words "INSERT EXPRESSION" to areas I must fill in. I have ran my code and almost debugged all the errors, but reciving an error message as such:

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate.

RCOND = 4.590499e-21.

> In BEM (line 137)

Error using horzcat

Dimensions of arrays being concatenated are not consistent.

Error in BEM (line 140)

err = norm([a-an ap-apn]);

^

Here is the code which i have partially debugged (I guess up until line 140). Also "BEM" is the name of the matlab code while I have screengrabs of the loaded files, MATLAB won't let me share them in so I am hoping I attache them within this post. This error is very weird and I'm sure it's a matrix or vector/scalar issue just not sure how to go about it but it seems to be originating from the variable "a":

%Description:  This script will estimate the performance coefficient and
%thrust coefficient over a range of tip-speed ratios using the blade
%element momentum theory formulation in the 'AeroDyn Theory Manual'
%authored by P.H. Moriarty and A.C. Hansen.  This simple code requires a
%table of relevant blade information in addition to three airfoil tables.
%A description of the contents of these tables, which are loaded from text
%files, can be found below.
%Clear workspace
clear all
%Set inputs (These inputs are for the NREL 5 MW reference wind turbine)
%Load airfoil data
R = 63; %Blade radius (m)
RH = 1.5; %Hub radius (m)
Uinf = 11.4; %Wind speed (m/s)
B = 3; %Number of blades (-)
rho = 1.225; %Density of air (kg/m^3)
RPMList = [2:1:16]; %Rotor speed (rpm)
%Load blade node data
%Columns: 1- radial position (m), 2 - twist (deg), 3 - dr (m)
%3 - chord (m), %4 - %airfoil # (-)
b = load('blade.txt'); %Blade specification is modified NREL 5 MW
%Columns: 1 - beta in deg, 2 - lift coefficient, 3 - drag coefficient,
%4 - pitching moment coefficient (not used).
a1 = load('airfoil1.txt'); %Airfoil 1 is a cylinder with CD = 0.5
a2 = load('airfoil2.txt'); %Airfoil 2 is a Delft University DU30
a3 = load('airfoil3.txt'); %Airfoil 3 is a NACA 64(3)-618
%Separate out blade information
r = b(:,1); %Node radial position
beta = b(:,2)*pi/180; %Node twist (converted to radians)
dr = b(:,3); %Element width dr
c = b(:,4); %Node chord
at = b(:,5); %Airfoil #
%Determine number of elements
ne = length(r);
%Compute local solidity
sp = (1/(2*pi))*B.*c./r;
%Set BEM solver solution criteria
tol = 0.005;
maxiter = 1000;
%Determine number of rotor speed (TSR) values to compute Cp and Ct values
NTSR = length(RPMList);
%Initialize TSR, Cp and Ct
TSR = zeros(NTSR,1); %Tip-speed ratio
Cp = zeros(NTSR,1); %Performance coefficient
Ct = zeros(NTSR,1); %Thrust coefficient
%Loop through rotor speed values
for j = 1:NTSR;
    
    %Set RPM
    RPM = RPMList(j);
    
    %Compute local tip speed ratio for all elements
    tsr = (RPM*2*pi/60)*(1/Uinf)*r;
    
    %Initialize total thrust and total torque
    Ttotal = 0; %Thrust (N)
    Qtotal = 0; %Torque (N-m)
    Ftotal = 0; %Rotor plane blade shear force (N)
    
    %Loop through blade elements to compute torque and thrust increment dT
    %and dQ
    for i = ne:-1:1;
        
        %Set initial guess for induction factors
        a = 0.25*(2+pi*tsr(i)*sp(i)...
            -sqrt(4-4*pi*tsr(i)*sp(i)+pi*(tsr(i)^2)*sp(i)...
            *(8*beta(i)+pi*sp(i))));  %Axial induction factor
        ap = 0;  %Tangential induction factor
        
        %Set iteration counter
        ct = 0;
        
        %Set initial error
        err = 100;
        
        %Iterate to find final BEM induction factors
        while ((err > tol) && (ct < maxiter));
            
            %Update counter
            ct = ct+1;
            
            %Compute inflow angle (rad)
            phi = atan(Uinf*(1-a)./((RPM*2*pi/60)*r(i)*(1+ap)));
            
            %Estimate angle of attack in degrees
            alpha = (180/pi)*(phi-beta(i));
            
            %Obtain lift and drag coefficients
            if (at(i) == 1);
                CL = interp1(a1(:,1),a1(:,2),alpha);
                CD = interp1(a1(:,1),a1(:,3),alpha);
            elseif (at(i) == 2);
                CL = interp1(a2(:,1),a2(:,2),alpha);
                CD = interp1(a2(:,1),a2(:,3),alpha);
            else; 
                CL = interp1(a3(:,1),a3(:,2),alpha);
                CD = interp1(a3(:,1),a3(:,3),alpha);                  
            end;
            
            %Estimate the thrust coefficient
            CT = (sp*(1-a)^2*(CL.*cos(phi)+CD.*sin(phi)))/sin(phi).^2;
            
            %Compute tip loss
            Ft = (2/pi)*cos(exp(-((B*(R-r))/2.*r*sin(phi)))).^-1;
            
            %Compute hub loss
            Fh = (2/pi)*acos(exp(-(B*(r(i)-RH)./(2*RH*sin(phi)))));
            
            %Compute total loss
            F = Ft.*Fh;
            
            %Compute updated axial induction factor
            if (CT <= 0.96*F);  %Blade is not highly loaded
                an = (1+(4*F.*(sin(phi).^2))/(sp(i)*(CL.*cos(phi)+CD.*sin(phi)))).^(-1);
            else;  %Blade is highly loaded (Glauert Correction)
                an = .4;
            end;
            
            %Compute updated tangential induction factor
            apn = (-1+(4*sin(phi).*cos(phi))/(sp*(CL.*sin(phi)-CD.*cos(phi)))).^-1;
            
            %Compute error
            err = norm([a-an ap-apn]);
            
            %Update induction factors
            a = an;
            ap = apn;
        end;
        
        %Compute Ustar
        Ustar = sqrt((Uinf*(1-a))^2+((RPM*2*pi/60)*r(i)*(1+ap))^2);
        
        %Update thrust
        dT = B*.5*rho*Ustar^2*c*(CL*cos(phi)+CD*sin(phi)); %Thrust increment
        Ttotal = Ttotal+dT; 
        
        %Update torque
        dQ = B*0.5*rho*(Ustar^2)*...
            (CL*sin(phi)-CD*cos(phi))*c(i)*r(i)*dr(i); %Torque increment
        Qtotal = Qtotal+dQ;
        
        %Obtain values need to compute blade edge moment
        dF = dQ/r(i);
        Ftotal = Ftotal+dF;
        
        %Store out-of-rotor-plane moment My
        if (i == ne);
            My(j,i) = (dT/B)*dr(i)/2;
        else;
            My(j,i) = My(j,i+1)+(dT/B)*dr(i)/2+((Ttotal-dT)/B)*dr(i);
        end;
        
        %Store in-rotor-plane moment Mx
        if (i == ne);
            Mx(j,i) = (dF/B)*dr(i)/2;
        else;
            Mx(j,i) = Mx(j,i+1)+(dF/B)*dr(i)/2+((Ftotal-dF)/B)*dr(i);
        end;    
        
    end;
    
    %Compute total power
    Ptotal = Qtotal*RPM*2*pi/60;
    
    %Compute TSR for rotor speed j
    TSR(j) = (RPM*2*pi/60)*R/Uinf;
    
    %Compute Cp for rotor speed j
    Cp(j) = 4*a*(1-a)^2;
    
    %Compute Ct for rotor speed j
    Ct(j) = Ttotal/(0.5*rho*pi*(R^2)*(Uinf^2));
    
end;

ANY HELP IS GREATLY APPRECIATED AND THANK YOU

Answer 1

[Walter Roberson]

b = load('blade.txt'); %Blade specification is modified NREL 5 MW

b will generally be a matrix

r = b(:,1); %Node radial position
c = b(:,4); %Node chord

r and c will generally be vectors.

sp = (1/(2*pi))*B.*c./r;

because r and c are generally vectors, sp will generally be a vector.

            apn = (-1+(4*sin(phi).*cos(phi))/(sp*(CL.*sin(phi)-CD.*cos(phi)))).^-1;

sp is generally a vector. You have a / operation in which sp appears in the denominator, so the / operator is matrix-right-divide rather than element-by-element division. It happens that the matrix right divide operator here is calculating on a singular matrix, so you get the RCOND error.