function ueg = refeval(p, zz, ii, side)
% Reference evaluation of u for a single domain with smooth segment boundary.
% Note: scattering BVP, ie u = layerpot + u_inc.
% Uses MPSpack, & FMMLIB (one bad strip call done with all pts in it together).
% Inputs: p = scattering object (one domain, one segment with layerpots on it)
%         zz = list of grid points
%         side = 'i' or 'e' for interior/exterior
% Barnett 11/3/12

N = numel(p.co);
de = p.doms(1); s = p.segs(1);

refine = 10;  % native eval outside bad strip, just to check
NN = refine*N;
co = p.co; p.updateN(NN); coNN = quadr.interptrig(co,NN); p.co = coNN;
clear coNN;
fprintf('uN via FMM...'); ueg = nan*zz; % u_native
tic, o.FMM=1; ueg(ii) = p.pointsolution(pointset(zz(ii)),o); toc
ueg(ii) = ueg(ii) + p.ui(zz(ii));
p.updateN(N); p.co = co; % return the segment and coeffs to normal (frees RAM)

% regenerate zzc=subset of zz in bad strip, ic = their indices
tic, fprintf('refeval: making zzc and ic...\n');
als = 5.0/N; % imag strip width wrt [0,1]-scaled seg, for acc to 1e-13
if side=='e', als = -als; end
v = s.Z(s.t + 1i*als); % imag-shifted N bdry pts. Now clip off loopy ones:
vsd = abs(repmat(s.x,[1 N]) - repmat(v.',[N 1])); % find dists v to pO, O(N^2)
[vsdm, ivsd] = min(vsd,[],1); clear vsd; % get indices of nearest pts
oi.toz = s.t(ivsd);  % nearest on-pO initial guesses to help invertZ
vs = s.invertZparam(v,oi); % slow if wobbly curve
v = v(find(abs(imag(vs))>0.99*abs(als))); % snip off too-close self-int loops
fc = ~utils.inpolywrapper(zz(ii), v); % for ii pts, true if in strip
if side=='e', fc = ~fc; end % remember inside the Gamma-curve for ext
ic = ii(find(fc)); clear fc; zzc = zz(ic);      % pts to correct

% setup for refined native eval in bad strip
refine = 1e3;
NN = refine*N;
co = p.co;
fprintf('updating N...\n'); p.updateN(NN);
fprintf('trig interp...\n'); coNN = quadr.interptrig(co,NN); p.co = coNN;
clear coNN;

% setup dyadic extrapolation weight vector ev for data at points h...
m=10; % number of data points % m=10
V = zeros(m,m); h = 2.^(0:(m-1));
for j=1:m, V(:,j) = h.^(j-1); end, iV = inv(V); ev = iV(:,1); % get const term
%sin(h*0.05*.25)*ev % test model func w/ known val at zero, of correct freq

nap = 100*N; vj = s.Z((1:nap)/nap - 5i/NN); % Gamma curve where NN eval bad
icj = find(utils.inpolywrapper(zzc, vj)); n = numel(icj); % pts to interp to
zi = nan(m,n);
fprintf('setup interp pts for %d bad pts... ',n)
for i=1:n         % loop over (few hundred) bad pts & set up their interp pts
  z = zzc(icj(i)); % bad pt @ z
  a = find(abs(z-s.x)==min(abs(z-s.x))); na = s.nx(a); % nearest surf normal
  mind = s.speed(a)*5/NN; % local closest dist while still accurate
  zi(:,i) = z + mind*h.'*na; % line of interp pts pointing out (exterior case)
end
if numel(find(isnan(zi)))>0, warning('warning: zi has some NaNs!\n'); end
nout = numel(find(~de.inside(zi)));
if nout>0, warning(sprintf('ouch: %d interp pts outside domain!\n',nout)); end
toc, fprintf('FMM %d sources to %d+%d targets... ',NN,numel(zzc),numel(zi))
tic, o.FMM=1; data = p.pointsolution(pointset([zzc; zi(:)]),o); toc
p.updateN(N); p.co = co; % return the segment and coeffs to normal (frees RAM)
uec = data(1:numel(zzc)); % get refined native eval data
% now do all interps via weighted sums with left-matvec, overwrite bad pts...
uec(icj) = ev.' * reshape(data(numel(zzc)+1:end),m,n);
clear data;
ueg(ic) = uec + p.ui(zzc); % overwrite bad strip: add u_inc to crazy-refined uec