function []=thinairfoil()
	%Thin Airfoil Theory

	%Initialize variables
	A=zeros(1,3);
	x=zeros(1,0);
	y=zeros(1,0);
	p=0;
	coeffor=ones(1,4);
	coefback=ones(1,4);

	%Choose the format to enter data.
	method=input('What way would you like to enter the airfoil data?\n 1) 4 digit airfoil number\n 2) two cubic equations\n 3) data points\n?');

	switch method
		case 1
			%Four digit airfoil method
			FourDigit=input('\nType in NACA 4 digit airfoil number.\n?');
			alpha=input('\nType in an angle of attack in degrees.\n?');
			alpha=alpha*pi/180;
			m=floor(FourDigit/1000)*1000;
			p=floor(FourDigit/100)*100-m;
			m=m/100000;
			p=p/1000;

			[TempA0,TempA1,TempA2,Tempx,Tempy]=FindA(0,-m/p^2,2*m/p,0,0,p);
			A=A+[TempA0,TempA1,TempA2];
			x=[x,Tempx];
			y=[y,Tempy];

			[TempA0,TempA1,TempA2,Tempx,Tempy]=FindA(0,-m/(1-p)^2,2*m*p/(1-p)^2,(m-2*m*p)/(1-p)^2,p,1);
			A=A+[TempA0,TempA1,TempA2];
			x=[x,Tempx];
			y=[y,Tempy];
		case 2
			%Two cubic equations of an airfoil
			while (coeffor(1)*p^3+coeffor(2)*p^2+coeffor(3)*p)~=(coefback(1)*p^3+coefback(2)*p^2+coefback(3)*p+coefback(4)) | dot(coefback,ones(1,4))~=0
				coeffor=input('\nType in coefficients for the forward camber of the airfoil.\n  ax^3 + bx^2 + cx\n  ex.[0 1 2] = x^2 + 2x\n?');
				p=input('\nEnter the point fraction of the chord where the forward camber ends and the back camber begins.\nIf only one equation is used type 1.\n?');
				if p~=1
					coefback=input('\nType in coefficients for the back camber of the airfoil.\n  ax^3 + bx^2 + cx + d\n  ex.[0 1 2 3] = x^2 + 2x + 3\n?');
				end
				if (coeffor(1)*p^3+coeffor(2)*p^2+coeffor(3)*p)~=(coefback(1)*p^3+coefback(2)*p^2+coefback(3)*p+coefback(4)) | dot(coefback,ones(1,4))~=0
					disp('');
					disp('Equations must connect and the y coordinate is 0 at x/c=1.');
				end
			end
		
			[TempA0,TempA1,TempA2,Tempx,Tempy]=FindA(coeffor(1),coeffor(2),coeffor(3),0,0,p);
			A=A+[TempA0,TempA1,TempA2];
			x=[x,Tempx];
			y=[y,Tempy];

			[TempA0,TempA1,TempA2,Tempx,Tempy]=FindA(coefback(1),coefback(2),coefback(3),coefback(4),p,1);
			A=A+[TempA0,TempA1,TempA2];
			x=[x,Tempx];
			y=[y,Tempy];
		case 3
			x=[0 0];
			y=[0 0];
			while x(1)~=0 | y(1)~=0 | x(length(x))~=1 | y(length(y))~=0
				x=input('\nPlease enter the x coordinates of the camber line.\n  ex.[0 .25 .5 .75 1]\n?');
				y=input('\nPlease enter the y coordinates of the camber line.\n  ex.[0 .05 .1 .05 0]\n?');
				if x(1)~=0 | y(1)~=0 | x(length(x))~=1 | y(length(y))~=0
					disp('');
					disp('Coordinates must begin at (0,0) and end at (1,0) for this theory.');
				end
			end
			for i=1:length(x)-1
				m=(y(i+1)-y(i))/(x(i+1)-x(i));
				b=-x(i)*m+y(i);
				[TempA0,TempA1,TempA2]=FindA(0,0,m,b,x(i),x(i+1));
				A=A+[TempA0,TempA1,TempA2];
			end
	end

	%Optional printing information to a file.
	if input('\nWould you like to print the information to a file?(y,n)\n?','s')=='y'
		filename=input('Enter a file name.\n?','s');
		filename=strcat(filename,'.dat');
		fid=fopen(filename,'w');
		fprintf(fid,'\n\nName                        Symbol   Value\n');
		fprintf(fid,'--------------------------  ------   -------\n');
		fprintf(fid,'Lift Coefficient            Cl      %9.5f\n',Cl);
		fprintf(fid,'Zero Lift Angle Name        a0l     %9.5f(%9.5f degrees)\n',a0l,a0l*180/pi);
		fprintf(fid,'Moment Coefficient (x=0  )  Cm0     %9.5f\n',Cm0);
		fprintf(fid,'Moment Coefficient (x=c/4)  Cmac    %9.5f\n',Cmac);
		fprintf(fid,'Center of Pressure          xcp     %9.5f\n',xcp);
		fclose(fid);
		disp(['Created filename',blanks(1),filename]);
	end

	%Plot camber line and arrow to show free stream angle of attack
	plot(x,y);
	hold on;
	axis([-.5,1.5,-.1,.1]);
	e=pi/18;
	plot([-0.1,-0.1-0.3*cos(alpha)],[0,-0.15*sin(alpha)],'k')
	plot([-0.1,-0.1-0.05*cos(e-alpha)],[0,0.025*sin(e-alpha)],'k')
	plot([-0.1,-0.1-0.05*cos(e+alpha)],[0,-0.025*sin(e+alpha)],'k')
	hold off;

	%Calculate aerodynamic effects
	a0l=-(A(1,1)+A(1,2)/2);               %Zero angle of attack(rad)
	A(1,1)=A(1,1)+alpha;
	Cl=2*pi*(A(1,1)+A(1,2)/2);             %Lift Coefficient
	Cm0=-pi/2*(A(1,1)+A(1,2)-A(1,3)/2);    %Moment Coefficient at tip
	Cmac=pi*(A(1,3)-A(1,2))/4;             %Moment Coefficient around aerodynamic center.
	xcp=((2*A(1,1)+2*A(1,2)-A(1,3))/(2*A(1,1)+A(1,2)))/4;      %Center of pressure

	%Print Information of Thin Airfoil
	disp(sprintf('\n\nName                        Symbol   Value'));
	disp(sprintf('--------------------------  ------   -------'));
	disp(sprintf('Lift Coefficient            Cl      %9.5f',Cl));
	disp(sprintf('Zero Lift Angle Name        a0l     %9.5f(%9.5f degrees)',a0l,a0l*180/pi));
	disp(sprintf('Moment Coefficient (x=0  )  Cm0     %9.5f',Cm0));
	disp(sprintf('Moment Coefficient (x=c/4)  Cmac    %9.5f',Cmac));
	disp(sprintf('Center of Pressure          xcp     %9.5f',xcp));
 
 
 
%********************************************************** 
function [A0,A1,A2,x,y]=FindA(a,b,c,d,x1,x2)

	x=x1:(x2-x1)/20:x2;
	y=a*x.*x.*x+b*x.*x+c*x+d;
	%Conversion to theta coordinates
	theta1=acos(1-2*x1);
	theta2=acos(1-2*x2);

	A01=1.125*a*theta1 + b*theta1 + c*theta1 - 1.5*a*sin(theta1) - b*sin(theta1) + 0.1875*a*sin(2*theta1);
	A02=1.125*a*theta2 + b*theta2 + c*theta2 - 1.5*a*sin(theta2) - b*sin(theta2) + 0.1875*a*sin(2*theta2);
	A0=(A01-A02)/pi;

	A11=-0.75*a*theta1 - 0.5*b*theta1 + 1.3125*a*sin(theta1) + b*sin(theta1)+c*sin(theta1) - 0.375*a*sin(2*theta1) - 0.25*b*sin(2*theta1) + 0.0625*a*sin(3*theta1);
	A12=-0.75*a*theta2 - 0.5*b*theta2 + 1.3125*a*sin(theta2) + b*sin(theta2)+c*sin(theta2) - 0.375*a*sin(2*theta2) - 0.25*b*sin(2*theta2) + 0.0625*a*sin(3*theta2);
	A1=2*(A12-A11)/pi;	

A21=0.1875*a*theta1 - 0.75*a*sin(theta1) - 0.5*b*sin(theta1) + 0.5625*a*sin(2*theta1) + 0.5*b*sin(2*theta1) + 0.5*c*sin(2*theta1) - 0.25*a*sin(3*theta1) - b*sin(3*theta1)/6+0.046875*a*sin(4*theta1);
A22=0.1875*a*theta2 - 0.75*a*sin(theta2) - 0.5*b*sin(theta2) + 0.5625*a*sin(2*theta2) + 0.5*b*sin(2*theta2) + 0.5*c*sin(2*theta2) - 0.25*a*sin(3*theta2) - b*sin(3*theta2)/6+0.046875*a*sin(4*theta2);
A2=2*(A22-A21)/pi;