Anyone able to find the error in this aerofoil code?

if true
  %943639  %uvwxyz
%NACAwxyz %chord length z %u angle of attack
clc
%
NACA = string([4 4 7 4]);
L = length(NACA);
c = 3;
if L == 4
M = str2double(NACA(1))/100; % maximun chamber
P = str2double(NACA(2))/10; % maximun chamber along the cord
t1 = str2double(NACA(3));
t2 = str2double(NACA(4)); 
tmax =((10*t1+t2)/100)*(c/10); % max thickness
end
%%Mean Camber Line
x = linspace(0,c,1000);
N = P*c;
if 0<=x<=(P*c)
  dZdX = ((M*2)/P.^2)*(P-(x/c));
  Z = ((M*c)/P.^2)*(2*P*(x/c)-(x/c).^2);
elseif (P*c)<=x<=c
  dZdX = ((M*2)/(1-P).^2)*(P-(x/c));
  Z = ((M*c)/(1-P).^2)*((1-2*P)+2*P*(x/c)-(x/c).^2);
end
%%Equation for Upper and Lower Surface
%
Y = (tmax/0.2)*(0.2969*sqrt(x/c)-0.1281*(x/c)-0.3516*((x/c).^2)+0.2843*((x/c).^3)-0.1015*((x/c).^4));
%
theta = atan(dZdX); 
%
xup = x-Y*sin(theta);
yup = Z+Y*cos(theta);
xlow = x+Y*sin(theta);
ylow = Z-Y*cos(theta);

 fprintf("Printing Graph\n");
plot(xup,yup,'g'); hold on
plot(xlow,ylow,'g'); hold on
fprintf("\nPrinting Mean Camber Line");
plot(xup,Z,'b-.'); hold on
plot(0:3,0); hold on
axis([-0.5 3.5 -0.4 0.7])

 end

The Error is:

Error using * Inner matrix dimensions must agree.

Error in Aerofoil (line 33) xlow = x+Y*sin(theta);

Thank you

[Jan]

All we see is the code. We cannot guess why you assume, that it contains an error. Do you get an error message or does the result differ from your expectations? Please be so kind and add more information by editing the question (not as comment or answer).

Why do you create a string at first only to convert it back to a double afterwards by str2double?

[David Oshidero]

I was trying to find a way so I can easily change the NACA values as I wanted to test it for 3 different values. However when i test it for the one which I need to compare the graph produced to the graph I already have its wrong.

[Jan]

I assume this is more direct:

NACA = [4 4 7 4];
M  = NACA(1) / 100; % maximun chamber
P  = NACA(2) / 10;  % maximun chamber along the cord
t1 = NACA(3);
t2 = NACA(4);

While "the graph is wrong" has a meaning for you, the readers in the forum cannot imagine, what you expect. Please read: https://www.mathworks.com/matlabcentral/answers/6200-tutorial-how-to-ask-a-question-on-answers-and-get-a-fast-answer.

Answer 1

[Jan]

if 0<=x<=(P*c)

This will most likely not do what you expect. Unfortunately there is no comment, which defines clearly, what the intention is. But I guess you want:

if 0<=x & x<=(P*c)

Note that "0<=x<=(P*c)" is evaluated from left to right: 0<=x is either TRUE or FALSE, which is treated as 1 or 0. In the next step you get "1<=P*c" or "0<=P*c".

The next problem is that x is a vector, but the condition of an if must be scalar. Therefore Matlab inserts this automatically:

if all(0<=x<=(P*c))

So you need either a for loop over all elements of x, or logical indexing:

index    = (0<=x & x<=P*c);
dZdX     = ((M*2)/P.^2)*(P-(x(index)/c));
Z(index) = ((M*c)/P.^2)*(2*P*(x(index)/c)-(x(index)/c).^2);

The same for "elseif (P*c)<=x<=c".

[David Oshidero]

I will also try to change this. Thank you

[David Oshidero]

For logical indexing would i use the if statements aswell?

[Alec Britton]

I am also stuck with what he's asking could you elaborate

[Jan]

Z    = zeros(size(x));  % Or NaN(size(x)) ?
dZdX = zeros(size(x));  % Or NaN(size(x)) ?

index       = (0<=x & x<=P*c);
dZdX(index) = ((M*2)/P.^2)*(P-(x(index)/c));
Z(index)    = ((M*c)/P.^2)*(2*P*(x(index)/c)-(x(index)/c).^2);

index    = (P*c<=x & x<=c)
dZdX     = ((M*2)/(1-P).^2)*(P-(x(index)/c));
Z(index) = ((M*c)/(1-P).^2)*((1-2*P)+2*P*(x(index)/c)-(x(index)/c).^2);

The TRUE and FALSE values of the logical indexing replace the need for a loop and the IF branches.