The Matlab function listed below returns the coordinates of an airfoil with a specified maximum thickness and turning angles. The thickness distribution used is from the NACA 4-series airfoils. The camber line is computed as a cubic polynomial before it is rotated to match the input turning angles.

This function takes 5 input arguments. $$t$$ is the approximate maximum thickness as a percentage of chord, $$x_{\text{turn}}$$ is the approximate location of the turning point as a percentage of chord, $$\theta_1$$ and $$\theta_2$$ are turning angles defined next, and $$N$$ is the number of chord-wise points to return. The definition of the angles $$\theta_1$$ and $$\theta_2$$ is demonstrated in this sketch: function [ return_struct ] = ...
airfoil_gen( t, x_turn, theta1, theta2, N )
% Generates an coordinates of turbine airfoil
% blades with specified thickness and turning
% angles

phi = (theta1 + theta2)/2;

b = [tan(phi);
tan(-phi);
0;
0];

A = [0 1 0 0;
0 1 2 3;
1 0 0 0;
0 1 2*x_turn 3*x_turn^2];

a = A\b;

ii = linspace(0, 1, N);
xc = (cos(pi - ii*pi) + 1)/2;

yc = a(1) + a(2)*xc + a(3)*xc.^2 + a(4)*xc.^3;
dyc = a(2) + 2*a(3)*xc + 3*a(4)*xc.^2;

% Use NACA 4-series thickness distribution
yt = 5*t*(0.2969*sqrt(xc) - 0.126*xc ...
- 0.3516*xc.^2 + 0.2843*xc.^3 - 0.1015*xc.^4);

xu = xc - yt.*sin(dyc);
yu = yc + yt.*cos(dyc);

xl = xc + yt.*sin(dyc);
yl = yc - yt.*cos(dyc);

rotation_angle = phi - theta1;

% Rotate to match input angles:
C = [cos(rotation_angle) sin(rotation_angle);
-sin(rotation_angle) cos(rotation_angle)];

for ii = 1:N
coords = [xc(ii); yc(ii)];
coords = C*coords;
xc(ii) = coords(1);
yc(ii) = coords(2);

coords = [xu(ii); yu(ii)];
coords = C*coords;
xu(ii) = coords(1);
yu(ii) = coords(2);

coords = [xl(ii); yl(ii)];
coords = C*coords;
xl(ii) = coords(1);
yl(ii) = coords(2);
end

return_struct.xl = xl;
return_struct.yl = yl;
return_struct.xu = xu;
return_struct.yu = yu;
return_struct.xc = xc;
return_struct.yc = yc;

end


These airfoils are definitely not optimal but they may be good to use as initial guesses in aerodynamic shape optimization. Below are several plots of airfoil shapes using the script above. The baseline parameters are $$t = 0.12$$, $$\theta_1 = 15^\circ$$, $$\theta_2 = 30^\circ$$, and $$x_{\text{turn}} = 0.35$$. Each of the plots vary these parameters one at a time.

Maximum thickness variation: $$\theta_1$$ variation: $$\theta_2$$ variation: $$x_{\text{turn}}$$ variation: If you have any comments or suggestions, let me know if the comments below.