function []=vorlat()
%Vortex lattice method.

	%This takes user input on the number of nodes to use on the
	%upper and lower surface and the 4 digit NACA number.
   clf
   writetofile=0;
   cr=input('Please input root chord as fraction of the span.\n?');
   ct=input('Please input tip chord as fraction of the span.\n?');
   L=input('Please input the leading edge sweepback angle in degrees.\n?');
   r=input('Please input the number panels as a matrix where the second number is even.\nEx. [2 6]\n?');
   numpanel=r(1)*r(2);
   
	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
   
   %Calculate endpoints of panels for visual representation
   XEndPoints=zeros(r(2)+1,r(1)+1);
   YEndPoints=zeros(r(2)+1,r(1)+1);
   for i=1:r(2)+1
      for j=1:r(1)+1
         XEndPoints(i,j)=0+(j-1)/r(1);
         YEndPoints(i,j)=-.5+(i-1)/r(2);
      end
   end
   XEndPoints=Transformation(XEndPoints,YEndPoints,ct,cr,L*pi/180);
   
   %Calculate vortex end points and control points
   XVor=(3*XEndPoints(:,1:r(1))+XEndPoints(:,2:r(1)+1))/4;
   YVor=YEndPoints(:,1:r(1));
   XContP=(XEndPoints(:,1:r(1))+3*XEndPoints(:,2:r(1)+1))/4;
   XContP=(XContP(1:r(2),:)+XContP(2:r(2)+1,:))/2;
   YContP=(YEndPoints(:,1:r(1))+YEndPoints(:,2:r(1)+1))/2;
   YContP=(YContP(1:r(2),:)+YContP(2:r(2)+1,:))/2;
   
   hold on
   plot(YVor,-XVor,'g-');
   plot(YContP,-XContP,'ro');
   plot(YEndPoints,-XEndPoints,'b-');
   plot(YEndPoints',-XEndPoints','b-');
   
   Xv=zeros(0,1);
   Yv=zeros(0,1);
   Xc=zeros(0,1);
   Yc=zeros(0,1);
   for j=1:2
	   for i=1:r(1)
         Xv=[Xv;XVor(r(2)/2*(j-1)+1:r(2)/2*j+1,i)];
         Yv=[Yv;YVor(r(2)/2*(j-1)+1:r(2)/2*j+1,i)];
         Xc=[Xc;XContP(r(2)/2*(j-1)+1:r(2)/2*j,i)];
         Yc=[Yc;YContP(r(2)/2*(j-1)+1:r(2)/2*j,i)];
      end
   end
   
   disp('Panel numbers and coordinates of the bound vortex end points.');
   for i=1:numpanel
   	disp(sprintf('%2d    %7.4f   %7.4f   %7.4f   %7.4f',i,Xv(i+floor(2*i/(r(2)+1))),Yv(i+floor(2*i/(r(2)+1))),Xv(1+i+floor(2*i/(r(2)+1))),Yv(1+i+floor(2*i/(r(2)+1)))));
   end
   disp(sprintf('\n\nControl point coordinates.'));
   for i=1:numpanel
   	disp(sprintf('%2d    %7.4f   %7.4f',i,Xc(i),Yc(i)));
   end
   
   if writetofile
	   fprintf(fid,'Panel numbers and coordinates of the bound vortex end points.\n');
   	for i=1:numpanel
   		fprintf(fid,'%2d    %7.4f   %7.4f   %7.4f   %7.4f\n',i,Xv(i+floor(2*i/(r(2)+1))),Yv(i+floor(2*i/(r(2)+1))),Xv(1+i+floor(2*i/(r(2)+1))),Yv(1+i+floor(2*i/(r(2)+1))));
	   end
   	fprintf(fid,'\n\nControl point coordinates.\n');
   	for i=1:numpanel
   		fprintf(fid,'%2d    %7.4f   %7.4f\n',i,Xc(i),Yc(i));
   	end
   end

   A=zeros(numpanel/2);
   for i=1:numpanel/2
      for j=1:numpanel
			d1=sqrt((Xv(j+floor(2*(j-1)/(r(2))))-Xc(i))^2+(Yv(j+floor(2*(j-1)/r(2)))-Yc(i))^ 2);
			d2=sqrt((Xv(1+j+floor(2*(j-1)/r(2)))-Xc(i))^2+(Yv(1+j+floor(2*(j-1)/r(2)))-Yc(i))^2);
			d3=(Xv(1+j+floor(2*(j-1)/r(2)))-Xc(i))*(Yv(j+floor(2*(j-1)/r(2)))-Yc(i))-(Xv(j+floor(2*(j-1)/r(2)))-Xc(i))*(Yv(1+j+floor(2*(j-1)/r(2)))-Yc(i));
			d4=Yv(j+floor(2*(j-1)/r(2)))-Yc(i);
			d5=Yv(1+j+floor(2*(j-1)/r(2)))-Yc(i);
			vbc=((Xv(1+j+floor(2*(j-1)/r(2)))-Xv(j+floor(2*(j-1)/r(2))))*(Xv(j+floor(2*(j-1)/r(2)))-Xc(i))+(Yv(1+j+floor(2*(j-1)/r(2)))-Yv(j+floor(2*(j-1)/r(2))))*(Yv(j+floor(2*(j-1)/r(2)))-Yc(i)))/d1;
			vbc=vbc-((Xv(1+j+floor(2*(j-1)/r(2)))-Xv(j+floor(2*(j-1)/r(2))))*(Xv(1+j+floor(2*(j-1)/r(2)))-Xc(i))+(Yv(1+j+floor(2*(j-1)/r(2)))-Yv(j+floor(2*(j-1)/r(2))))*(Yv(1+j+floor(2*(j-1)/r(2)))-Yc(i)))/d2;
			vbc=vbc/d3;
			vab=1-(Xv(j+floor(2*(j-1)/r(2)))-Xc(i))/d1;
			vab=vab/d4;
			vcd=1-(Xv(1+j+floor(2*(j-1)/r(2)))-Xc(i))/d2;
			vcd=-vcd/d5;
         w=vab+vbc+vcd;
			if j>numpanel/2
				m=fix((j-numpanel/2-.5)/(r(2)/2))+1;
            n=(2*m-1)*(r(2)/2)+1-(j-numpanel/2);
				A(i,n)=A(i,n)+w;
			else
				A(i,j)=w;
         end
      end
   end
   disp('Calculations from this point on are for the left wing only.');
   disp('Values for the right wing may be obtained using the wing summetry.');
   
   VorStr=-inv(A)*ones(numpanel/2,1);
   disp(sprintf('\nVortex Strengths:'));
   for i=1:numpanel/2
      disp(sprintf('G(%1d)/(4*pi*b*U*alpha)   %7.4f',i,VorStr(i)));
   end
   disp(sprintf('\nWing Lift Curve Slope              %7.4f',32*pi/((cr+ct)*r(2))*ones(1,numpanel/2)*VorStr));
   disp(sprintf('Wing Area                          %7.4f',(cr+ct)/2));
   disp(sprintf('Aspect Ratio                       %7.4f',2/(cr+ct)));
   disp(sprintf('Taper Ratio                        %7.4f',ct/cr));
   disp(sprintf('Leading Edge Sweepback Angle       %7.4f',L));
   disp(sprintf('Trailing Edge Sweepback Angle      %7.4f',atan((.5*tan(L*pi/180)+ct-cr)/.5)*180/pi));
   disp(sprintf('Quater Chord Line Sweepback Angle  %7.4f',atan((.5*tan(L*pi/180)+(ct-cr)/4)/.5)*180/pi));
   
   if writetofile 	
		fprintf(fid,'\nVortex Strengths:\n');
   	for i=1:numpanel/2
      	fprintf(fid,'G(%1d)/(4*pi*b*U*alpha)   %7.4f\n',i,VorStr(i));
   	end
   	fprintf(fid,'\nWing Lift Curve Slope              %7.4f\n',32*pi/((cr+ct)*r(2))*ones(1,numpanel/2)*VorStr);
   	fprintf(fid,'Wing Area                          %7.4f\n',(cr+ct)/2);
   	fprintf(fid,'Aspect Ratio                       %7.4f\n',2/(cr+ct));
   	fprintf(fid,'Taper Ratio                        %7.4f\n',ct/cr);
   	fprintf(fid,'Leading Edge Sweepback Angle       %7.4f\n',L);
   	fprintf(fid,'Trailing Edge Sweepback Angle      %7.4f\n',atan((.5*tan(L*pi/180)+ct-cr)/.5)*180/pi);
   	fprintf(fid,'Quater Chord Line Sweepback Angle  %7.4f\n',atan((.5*tan(L*pi/180)+(ct-cr)/4)/.5)*180/pi);
	   fclose(fid);
		disp(['Created filename',blanks(1),filename]);
   end

function [xnew]=Transformation(x,y,ct,cr,L);
	tan_u=(.5.*tan(L)+ct.*x-cr.*x)/.5;
	xnew=cr.*x+abs(y).*tan_u;