%  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 
%  arclength.  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.
%

%  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
soln = zeros(npts,1);                         % Dirichlet data

%  Will use MATLAB routine fzero (which solves a scalar nonlinear 
%  equation via bisection) to determine placement of boundary points.
%  Function 'spacing' determines the difference between the arclength
%  from theta_min to theta_max and the desired spacing in arclength.
%  It uses function 'inteval' to evaluate the integrand.
%  Spacing of points is equal to only about 1.e-8.  Hence we should
%  see superalgebraic convergence to about this level of accuracy but
%  not further.

global aglob                                  % Used in 'inteval'
global bglob
global hglob                                  % Used in 'spacing'
global theta_min

aglob = a;  bglob = b;
L = quad8('inteval', 0, 2*pi, 1.e-8);         % Perimeter of ellipse
h = L/npts;                                   % Spacing between points
hglob = h;

theta = 0;                                    % Set initial point.
cost = cos(theta);  sint = sin(theta);
x(1) = a*cost;  y(1) = b*sint;
sig = sqrt(a^2*sint^2 + b^2*cost^2);
nx(1) = b*cost/sig;  ny(1) = a*sint/sig;
curv(1) = a*b/(sig^3);
soln(1) = 1 + 2*(x(1)-x0)/((x(1)-x0)^2+(y(1)-y0)^2);
theta_min = theta;

for i=2:npts,
  theta_ends = [theta; 2*pi];
  theta = fzero('spacing', theta_ends, optimset('TolX', 2.e-8));
  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);
  soln(i) = 1 + 2*(x(i)-x0)/((x(i)-x0)^2+(y(i)-y0)^2);
  theta_min = theta;
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)*h;
  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)*h;
  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),