%  This code solves Laplace's equation inside an ellipse using
%  a boundary integral equation.  The integral is approximated
%  using the trapezoidal rule, with points equally spaced in 
%  theta.  The dense linear system is solved directly for the
%  dipole density mu.  This is then used to evaluate the solution
%  at desired points inside the ellipse.
%
%  *NOTE* Since points are not equispaced in arc length, will not
%  see superalgebraic convergence (except for a circle).  In general
%  convergence should be second order.


%  Read input data.

npts = input(' Enter number of boundary points: ');
a = input(' Enter half major axis: ');
b = input(' Enter half minor axis: ');

x0 = 2*a;  y0 = 2*b;                          % any point outside ellipse

%  Set boundary points, normal derivatives, curvatures, and
%  arclengths.  Also set exact solution on the boundary.

x = zeros(npts,1);  y = zeros(npts,1);        % boundary points
nx = zeros(npts,1);  ny = zeros(npts,1);      % boundary normals
curv = zeros(npts,1);                         % curvatures
h = zeros(npts,1);                            % arclengths between points
soln = zeros(npts,1);                         % Dirichlet data

dtheta = (2*pi)/npts;
for i=1:npts,
  theta = i*dtheta;
  cost = cos(theta);  sint = sin(theta);
  x(i) = a*cost;  y(i) = b*sint;
  sig = sqrt(a^2*sint^2 + b^2*cost^2);
  nx(i) = b*cost/sig;  ny(i) = a*sint/sig;
  curv(i) = a*b/(sig^3);

%    Compute arclength by integrating over this piece.
  theta0 = theta - dtheta;
  nintpts = 10;
  dt = dtheta/nintpts;
  h(i) = 0;
  for ii=1:nintpts,
    t = theta0 + ii*dt;
    tmdt = t - dt;
    dxdtr = -a*sin(t);  dydtr = b*cos(t);
    dxdtl = -a*sin(tmdt);  dydtl = b*cos(tmdt);
    h(i) = h(i) + dt*.5*( sqrt(dxdtl^2+dydtl^2) + sqrt(dxdtr^2+dydtr^2) );
  end;

  soln(i) = 1 + 2*(x(i)-x0)/((x(i)-x0)^2+(y(i)-y0)^2);
end;

%  Form dense matrix.

A = zeros(npts,npts);  

for i=1:npts,
  for j=1:npts,
    if j==i,
      A(i,j) = .5*curv(j);
    else
      rx = x(j) - x(i);  ry = y(j) - y(i);
      rsq = rx^2 + ry^2;
      A(i,j) = ((rx*nx(j) + ry*ny(j))/rsq);
    end;
    jp1 = mod(j,npts) + 1;
    A(i,j) = A(i,j)*.5*(h(j)+h(jp1));
  end;
end;

A = eye(npts,npts) + (1/pi)*A;  
f = 2*soln;

%  Solve linear system.

mu = A\f;
    
%  Evaluate solution at a few points inside the the ellipse and 
%  determine the error.

testpt(1,:) = [0,0];  testpt(2,:) = [.5*x(1),.5*y(1)]; 
testpt(3,:) = input(' Enter a test point: ');
u = zeros(3,1);
for j=1:npts,
  for k=1:3,
    rx = x(j) - testpt(k,1);  ry = y(j) - testpt(k,2);
    rsq = rx^2 + ry^2;
    jp1 = mod(j,npts) + 1;
    u(k) = u(k) + mu(j)*((rx*nx(j) + ry*ny(j))/rsq)*.5*(h(j)+h(jp1));
  end;
end;
u = u/(2*pi);
for k=1:3,
  utrue(k,1) = 1 + 2*(testpt(k,1)-x0)/((testpt(k,1)-x0)^2+(testpt(k,2)-y0)^2);
end;
[u, utrue], diff = utrue - u, err=norm(utrue-u),