Finding the corresponding combination of values for a given outcome

I have created this code to calculate the thickness and camber of an airfoil by measuring two points (y1 and y2) in the upper and lower surface (z_airfoil is the surface). My desired output is a combination of values of "a" and "b" that give me the desired thickness(trel1,trel2,trel3) and camber(hmax1,hmax2). I need 6 pairs of values of "a" and "b" as I want to compare the thickness_rel and camber values I measured with: trel1 and hmax1, trel2 and hmax1, trel3 and hmax1, and the same combination with hmax2. I don't know how to get this as I get confused when it comes to the indices inside the for loops. Could someone help me figuring out this?

% Input data (assumed)
a = 0:0.1:1;
b = 0:0.1:1;

% Desired relative thickness and camber
trel1 = 0.08;   % 8%
trel2 = 0.14;   % 14%
trel3 = 0.25;   % 25%

hmax1 = 0;     % 0%
hmax2 = 0.08;  % 8%

%Flow along contour of Circle and Joukowsky Airfoil
theta1 = pi:-0.01:0;    %radians
theta2 = pi:0.01:2*pi;

for j=2:length(a)
    for k=2:length(b)
        
    c(j,k) = -a(j) + sqrt(d^2-b(k)^2);
    
%Up surface
for i=1:length(theta1)
R_up(i,j,k) = -(a(j)*cos(theta1(i))-b(k)*sin(theta1(i))) + sqrt((a(j)*cos(theta1(i)))^2-(2*(a(j)*cos(theta1(i)))*(b(k)*sin(theta1(i))))+(b(k)*sin(theta1(i)))^2+c(j,k)^2+(2*c(j,k)*a(j)));
z_circle_up(i) = R_up(i,j,k)*(cos(theta1(i))+1i*sin(theta1(i)));
z_airfoil_up(i) = z_circle_up(i)+c(j,k)^2/z_circle_up(i);
x1(i)=real(z_airfoil_up(i));
y1(i)=imag(z_airfoil_up(i));
end

%Down surface
for i=1:length(theta2)
R_down(i,j,k) = -(a(j)*cos(theta2(i))-b(k)*sin(theta2(i))) + sqrt((a(j)*cos(theta2(i)))^2-(2*(a(j)*cos(theta2(i)))*(b(k)*sin(theta2(i))))+(b(k)*sin(theta2(i)))^2+c(j,k)^2+(2*c(j,k)*a(j)));
z_circle_down(i) = R_down(i,j,k)*(cos(theta2(i))+1i*sin(theta2(i)));
z_airfoil_down(i) = z_circle_down(i)+c(j,k)^2/z_circle_down(i);
x2(i)=real(z_airfoil_down(i));
y2(i)=imag(z_airfoil_down(i));
end

%Measuring the distance

l(j,k) = 4*(a(j)+c(j,k))^2/(2*a(j)+c(j,k));    %chord

thickness_rel = max(abs(y2-y1))./l;
camber = max(abs(0.5*(y2+y1)))./l;
    end
end

[Bob Thompson]

What is d?

Also, I'm not familiar with the shape of a Joukowsky Airfoil. How do you know that the code isn't working? I was able to get it to run with substituting d = 1, so I'm guessing the shape is wrong some how, but more information would be helpful.

[Maria Sarcos]

Yes sorry d=1!, I'm not sure how to find the values of a and b that I want.. I assume I need several if statements to compare the thickness and camber calculated to the desired ones and then interpolate to find a and b, just don't know how to implement it..

[Bob Thompson]

Hmm, ok. How are you actually comparing thickness and camber to trel and hmax? The calculated results are coming out as 11 x 11 arrays, and the inputs are single values.

[Maria Sarcos]

I would need to interpolate through the whole array the value for trel I want and then find the index and value corresponding to a and b, don't really know how to code that

[Bob Thompson]

Hmm, ok. One of the first problems you're going to run into is that the calculated thickness values do not bound the trel values. A quick way to code a check if trel is within the thickness results could be the following:

if min(thickness_re,[],'all')<trel1 & max(~isinf(thickness_rel)==1,[],'all')>trel3
    % If you made it here, then the thickness values can be interpolated.
end