function [AR,xLE]=ARfit(P,t)

global tmax
% This function calculates the AR of a rectangle fitted inside an airfoil
% described by the cubic spline function in P. The value t is the thickness
% if the required rectangle (<t/c max).
% OUTPUT:
%                               AR=Aspect ratio of fitted rectangle
%                               xLE=Forward x-position of rectangle.

% Penalty functions for going out of bounds

flag=0;
x1=tmax/1000;
if (t>tmax-x1)
   tsave=t;
   t=tmax-x1;
   flag=2;
  
elseif (t<x1)
        tsave=t;
   t=x1;
   flag=1;
end;
   
N=50;
theta=(pi-logspace(-1,pi,N))/2;
theta(1)=[]; 
theta=[fliplr(theta) pi-theta];
f1=ppval(P,theta).*sin(theta)-t/2;
f2=[f1(2:N*2-2) f1(N*2-2)];
FF=abs(sign(f1)-sign(f2));
roots=find(FF>0);
thetaTE=reversefit(P,t/2,theta(roots(1)));
rTE=ppval(P,thetaTE);
xTE=rTE*cos(thetaTE);
thetaLE=reversefit(P,t/2,theta(roots(2)));
rLE=ppval(P,thetaLE);
xLE=rLE*cos(thetaLE);
l=abs(xTE-xLE);
AR=l/t;
   
if flag==1
   x3=2*x1;
   y1=AR;
   y2=1.1*AR;
   a = (y1-y2)/(-2*x3*x1+x1^2);
   b=(y1-a*x1^2-y2)/x1;
   AR=a*tsave^2+b*tsave+y2;
elseif flag==2
   x3=tmax-x1*.001;
   x2=tmax;
   x1=tmax-2*x1;
   a=-AR/((x2^2-x3^2)-2*x1*(x2-x3));
   b=-2*a*x1;
   c=-a*x2^2-b*x2;
   AR=a*tsave^2+b*tsave+c;
end;