function []=airfoil()
%This program computes velocity and pressure coefficient
%around an airfoil, whose contour is approximated by
%vortex panels of linearly varying strength.

	clf
	writetofile=0;
   hold off
	x=input('\nEnter the x values of the contour from tail to nose to tail.\nEx.[1 .5 0 .5 1]\n?');
	y=input('\nEnter the y values of the contour from tail to nose to tail.\nEx.[0 -.1 0 .1 0]\n?');
	%Initialize
   numpts=length(x);
   An=zeros(numpts);
   At=zeros(numpts);
  	cn1=zeros(numpts-1);
   cn2=zeros(numpts-1);
   ct1=zeros(numpts-1);
   ct2=zeros(numpts-1);
   v=zeros(numpts-1,1);
   Cp=zeros(numpts-1,1);
   x=x';
   y=y';
   
	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');
   	writetofile=1;
   end
   xmid=(x(2:numpts)+x(1:numpts-1))/2;
   ymid=(y(2:numpts)+y(1:numpts-1))/2;
   dx=x(2:numpts)-x(1:numpts-1);
   dy=y(2:numpts)-y(1:numpts-1);
   s=sqrt(dx.*dx+dy.*dy);
   theta=atan2(y(2:numpts)-y(1:numpts-1),x(2:numpts)-x(1:numpts-1));
   for i=1:numpts-1
      for j=1:numpts-1
         if i~=j
	         a=-(xmid(i)-x(j))*cos(theta(j))-(ymid(i)-y(j))*sin(theta(j));
   	      b=(xmid(i)-x(j))^2+(ymid(i)-y(j))^2;
      	   c=sin(theta(i)-theta(j));
	         d=cos(theta(i)-theta(j));
   	      e=(xmid(i)-x(j))*sin(theta(j))-(ymid(i)-y(j))*cos(theta(j));
      	   f=log(1+s(j)*(s(j)+2*a)/b);
         	g=atan2(e*s(j),b+a*s(j));
	         p=(xmid(i)-x(j))*sin(theta(i)-2*theta(j))+(ymid(i)-y(j))*cos(theta(i)-2*theta(j));
   	      q=(xmid(i)-x(j))*cos(theta(i)-2*theta(j))-(ymid(i)-y(j))*sin(theta(i)-2*theta(j));
      	   cn2(i,j)=d+.5*q*f/s(j)-(a*c+d*e)*g/s(j);
         	cn1(i,j)=.5*d*f+c*g-cn2(i,j);
	         ct2(i,j)=c+.5*p*f/s(j)+(a*d-c*e)*g/s(j);
   	      ct1(i,j)=.5*c*f-d*g-ct2(i,j);
         else
            cn1(i,j) = -1.0;
         	cn2(i,j) = 1.0;
         	ct1(i,j) = .5*pi;
         	ct2(i,j) = .5*pi;
         end
      end
   end
   
   %Create normal and tangent A matrices
   for i=1:numpts-1
      An(i,1)=cn1(i,1);
      An(i,numpts)=cn2(i,numpts-1);
      At(i,1)=ct1(i,1);
      At(i,numpts)=ct2(i,numpts-1);
      for j=2:numpts-1
         An(i,j)=cn1(i,j)+cn2(i,j-1);
         At(i,j)=ct1(i,j)+ct2(i,j-1);
      end
   end
   An(numpts,1)=1;
   An(numpts,numpts)=1;
   
   %Plot Airfoil Shape
   subplot(2,2,1),plot(x,y,'b-',x,y,'ro');
   subplot(2,2,1),axis([-.5,1.5,-.5,.5]);
   subplot(2,2,1),title('Shape of Airfoil');
   subplot(2,2,1),xlabel('x/c');
   subplot(2,2,1),ylabel('y/c');
   subplot(2,2,4),hold on
   if writetofile
      fprintf(fid,'    Shape of body\n');
		fprintf(fid,'     x            y\n');
		fprintf(fid,' ---------    --------\n');
      for i=1:length(x)
	   	fprintf(fid,'%9.5f    %9.5f\n',x(i),y(i));
      end
   end
   
   %Write information to screen
   alpha=input('\nEnter an angle of attack between -90 and 90 degrees.\n?');
   while alpha<=90 & alpha>=-90
	   alphar=alpha*pi/180;
      
      %Solve for gammas on panels
   	rhs=[sin(theta-alphar);0];
   	gamma=inv(An)*rhs;
   
   	for i=1:numpts-1
      	v(i)=cos(theta(i)-alphar);
      	for j=1:numpts
         	v(i)=v(i)+At(i,j)*gamma(j);
      	end
   	end
   	Cp=1-v.*v;
   
   	Cfx=Cp'*dy;
   	Cfy=-Cp'*dx;
   	Cm=Cp'*(dx.*xmid+dy.*ymid);
   
   	Cd=Cfx*cos(alphar)+Cfy*sin(alphar);
   	Cl=Cfy*cos(alphar)-Cfx*sin(alphar);
            
      disp('    Resulting Data around Airfoil');
   	disp('     i     x        y     theta       s     gamma       v       Cp');
      disp('    --- ------- -------- --------  ------- -------- -------- --------');
      for i=1:numpts
      	if i==numpts
         	String=sprintf('    %2d                            %7.4f',i,gamma(i));
         	disp(String);
      	else
         	String=sprintf('    %2d  %.4f  %7.4f  %7.4f  %7.4f  %7.4f  %7.4f  %7.4f',i,xmid(i),ymid(i),theta(i),s(i),gamma(i),v(i),Cp(i));
         	disp(String);
         end
      end
   	disp(['    Cl =',num2str(Cl)]);
   	disp(['    Cd =',num2str(Cd)]);
      disp(['    Cm =',num2str(Cm)]);
	   if writetofile
      for i=1:numpts
      	if i==numpts
         	String=sprintf('    %2d                            %7.4f',i,gamma(i));
         	disp(String);
      	else
         	String=sprintf('    %2d  %.4f  %7.4f  %7.4f  %7.4f  %7.4f  %7.4f  %7.4f',i,xmid(i),ymid(i),theta(i),s(i),gamma(i),v(i),Cp(i));
         	disp(String);
         end
      end
   	disp(['    Cl =',num2str(Cl)]);
   	disp(['    Cd =',num2str(Cd)]);
      disp(['    Cm =',num2str(Cm)]);
      	fprintf(fid,'\n\nAngle of Attack  =  %9.3f\n',alpha);
      	fprintf(fid,'\n    Resulting Data around Airfoil\n');
   		fprintf(fid,'     i     x        y     theta       s     gamma       v       Cp\n');
			fprintf(fid,'    --- ------- -------- --------  ------- -------- -------- --------\n');
      	for i=1:length(xmid)
	   		fprintf(fid,'    %2d  %.4f  %7.4f  %7.4f  %7.4f  %7.4f  %7.4f  %7.4f\n',i,xmid(i),ymid(i),theta(i),s(i),gamma(i),v(i),Cp(i));
      	end
   		fprintf(fid,'Cl = %9.5f\n',Cl);
   		fprintf(fid,'Cd = %9.5f\n',Cd);
   		fprintf(fid,'Cm = %9.5f\n',Cm);
   	end
      subplot(2,2,3),plot([xmid;xmid(1)],-[Cp;Cp(1)],'b-',[xmid;xmid(1)],-[Cp;Cp(1)],'ro');
      subplot(2,2,3),title('Pressures');
      subplot(2,2,3),xlabel('x/c');
      subplot(2,2,3),ylabel('-Cp');
      subplot(2,2,4),plot(alpha,Cl,'bd');
      subplot(2,2,4),plot(alpha,Cd,'rx');
      subplot(2,2,4),plot(alpha,Cm,'g*');
	   subplot(2,2,4),title('Force and Moment Coefficients');
	   subplot(2,2,4),xlabel('Angle of Attack(degrees)');
   	subplot(2,2,4),ylabel('Cl,Cd,Cm');
	   alpha=input('\nEnter another angle of attack.\nIf you are done type in 91.\n?');
   end
   if writetofile 	
	   fclose(fid);
		disp(['Created filename',blanks(1),filename]);
   end