function [Dir,Neu] = Mie_traces(k0,kd,curve0)
%
% determines the solution for the scattered field from the Fourier-Hankel coefficients
% for the expansion of the TM wave scattered by a circular dielectric cylinder of radius 1
%
%    k      : wavenumber in the exterior
%    kd     : wavenumber in the dielectric
%    curve0 : curve points where the solutions are evaluated
%    c      : vector with coefficients of the Fourier-Hankel coefficients
%    Dir    : Exterior Dirichlet trace
%    Neu    : Exterior Neumann trace
%
%   Please send an email to cjerez@sam.math.ethz.ch if you find
%   errors in the code or have problems using it. 

x = curve0(1,:);
y = curve0(2,:);

[theta,r0] = cart2pol(x,y);

Z   = kd/k0;         % contrast 

N   = round(5*kd); % Number of modes being retained

BR  = besselj(0:(N+1), k0);
YR  = bessely(0:(N+1), k0);
HR  = BR + 1i*YR;              % Hankel polynomial

BN  = besselj(0:(N+1), kd);

% Derivatives of Bessel functions (use derbes function)

DBR = derbes(BR); 
DBN = derbes(BN); 
DHR = derbes(HR); 


c = zeros(N+1,1);

for m=1:(N+1),

 c(m) = (BN(m)*DBR(m)-Z*DBN(m)*BR(m))/(Z*DBN(m)*HR(m)-BN(m)*DHR(m));

end

% Dirichlet trace over the boundary

Dir = zeros(size(theta));

BR0  = besselj(0:(N+1), k0*r0);
YR0  = bessely(0:(N+1), k0*r0);
HR0  = BR0 + 1i*YR0;              % Hankel polynomial
DHR0 = derbes(HR0); 



Dir = c(1)*HR0(1);

for m=2:N+1,
    Dir = Dir + (1i)^(-(m-1))*c(m)*HR0(m)*exp(1i*(m-1)*theta);
    Dir = Dir + (-1)^((m-1))*(1i)^((m-1))*c(m)*HR0(m)*exp(-1i*(m-1)*theta);
end

% Neumann trace over the boundary 

Neu = zeros(size(theta));

Neu = c(1)*DHR0(1);

for m=2:N,
    Neu = Neu + (1i)^(-(m-1))*c(m)*DHR0(m)*exp(1i*(m-1)*theta);
    Neu = Neu + (-1)^((m-1))*(1i)^((m-1))*c(m)*DHR0(m)*exp(-1i*(m-1)*theta);
end