%
% Generate Gambit Vertex map for NACA airfoils
% Author Kyle Crandall
%

clear all;
clc;

m = 0;
p = 0;
t = 15;
c = 1;
Xthe = 0;
Ythe = 0;
Zthe = 0;
points = 25;
alpha_i = 0;

name = sprintf('NACA%d%d%d', m, p, t);

Xthe = Xthe*pi()/180;
Ythe = Ythe*pi()/180;
Zthe = Zthe*pi()/180;
m = m/100;
p = p/10;
t = t/100;
alpha = alpha_i * -pi/180;
R = [cos(alpha), -sin(alpha);sin(alpha), cos(alpha)];

x = zeros(1, points);
y_c = zeros(1, points);
y_t = zeros(1, points);
U = zeros(3, points);
L = zeros(3, points);
U_tilt = U;
L_tilt = L;
out = zeros(points*2 - 2, 3);

Rx = [1 0 0; 0 cos(Xthe) -sin(Xthe); 0 sin(Xthe) cos(Xthe)];
Ry = [cos(Ythe) 0 sin(Ythe); 0 1 0; -sin(Ythe) 0 cos(Ythe)];
Rz = [cos(Zthe) -sin(Zthe) 0; sin(Zthe) cos(Zthe) 0; 0 0 1];

for i = 2:points
	x(i) = (1-cos(((i - 1) * (90/(points - 1))) * pi/180))*c;
	%x(i) = x(i - 1) + c/points;
	y_t(i) = (t*c/0.2)*((0.2969*((x(i)/c)^0.5))-(0.126*(x(i)/c))-(0.3516*((x(i)/c)^2))+(0.2843*((x(i)/c)^3))-(0.1015*((x(i)/c)^4)));
	
	if x(i) < p
		y_c(i) = (m * (x(i)/(p^2)))*(2*p - x(i)/c); 
	else
		y_c(i) = m*(c-x(i))/((1-p)^2) * (1 + x(i)/c - 2*p); 
	end
	
	theta = (y_c(i) - y_c(i - 1))/(x(i) - x(i-1));
	
	U(1, i) = x(i) - y_t(i)*sin(theta);
	U(2, i) = y_c(i) + y_t(i)*cos(theta);
	
	L(1, i) = x(i) + y_t(i)*sin(theta);
	L(2, i) = y_c(i) - y_t(i)*cos(theta);
	
	U_tilt(1:2, i) = R * U(1:2, i);
	L_tilt(1:2, i) = R * L(1:2, i);
	
	out(i, :) = U_tilt(:, i)';
	
	if i ~= points
		out(length(out) - (i-2), :) = L_tilt(:, i)';
	end
end

% out(length(out) - 3, :) = [-3*c,  2*c, 0];
% out(length(out) - 2, :) = [-3*c, -2*c, 0];
% out(length(out) - 1, :) = [ 3*c, -2*c, 0];
% out(length(out) - 0, :) = [ 3*c,  2*c, 0];

figure(1);
clf;
hold on;
plot(x, y_c);
plot(U(1,:), U(2, :), 'r');
plot(L(1,:), L(2, :), 'g');
plot([U(1, 1), U(1, length(U))], [U(2, 1), U(2, length(U))], 'k');
hold off;
title(name);
axis([0, c, -c, c]);

figure(2);
clf;
hold on
plot(U_tilt(1,:), U_tilt(2, :), 'r');
plot(L_tilt(1,:), L_tilt(2, :), 'g');
plot([U_tilt(1, 1), U_tilt(1, length(U))], [U_tilt(2, 1), U_tilt(2, length(U))], 'k');
hold off;
title(sprintf('%s: with attack %d deg', name, alpha));
axis([0, c, -c, c]);
out = Rx*Ry*Rz*out';
out = out';
dlmwrite(sprintf('%s_%d_%d.txt', name, alpha_i, c), out, ',');