function [cp,cl]=dpan(n,alpha,x,y)

% Doublet panel method (c) N.D.Sandham 2006
% inputs: n=number of points, alpha=incidence, xb,yb=panel end point co-ordinates
% output: Cp=pressure at xb,yb, Cl=lift coefficient
%
twopiinv=1.0/(2.0*pi);
for i=1:n
  xc(i)=0.5*(x(i)+x(i+1));
  yc(i)=0.5*(y(i)+y(i+1));
  dx=x(i+1)-x(i);
  dy=y(i+1)-y(i);
  beta(i)=atan2(dy,dx);
  panlen(i)=sqrt(dx^2+dy^2);
end
%
for i=1:n
  for j=1:n 
    xt=xc(i)-x(j);
    yt=yc(i)-y(j);
    x1=xt*cos(beta(j))+yt*sin(beta(j));
    y1=-xt*sin(beta(j))+yt*cos(beta(j));
    a(i,j)=twopiinv*(atan2(y1,x1)-atan2(y1,x1-panlen(j)));
  end
% self-influence coeff  
  a(i,i)=0.5;
% wake panel (i.e. effect of vortex at trailing edge)
  a(i,n+1)=twopiinv*atan2((yc(i)-y(n+1)),(x(n+1)-xc(i)));
  rhs(i)=-(xc(i)*cos(alpha)+yc(i)*sin(alpha));
end

% Kutta condition
a(n+1,1)=1.0;
a(n+1,n)=-1.0;
a(n+1,n+1)=-1.0;
rhs(n+1)=0.0;

% solve for doublet strengths
sol=a\rhs';

% find lift coefficient
cl=2.0*sol(n+1);

% generate pressure coefficients
for i=2:n
    qt(i)=-(sol(i)-sol(i-1))*2.0/(panlen(i-1)+panlen(i));
    cp(i)=1.0-qt(i)^2;
end    
cp(1)=1.0;
cp(n+1)=1.0;