% Numerical inversion of Z(s) function.
% superceded by segment.invertZparam method.  Barnett 3/8/12

clear all classes; N = 100;
s = segment.smoothnonsym(N, 0.3, 0.3, 3);

% ========= some other shapes, without curvature terms ============
%a3=0.03; a4=0.85;    % crescent shape (disc w/ bump) shows explosion once
%s = segment(N, {@(t) exp(2i*pi*t)-a3./(a4+exp(2i*pi*t)), ...
%                @(t) 2i*pi*exp(2i*pi*t).*(1 + a3./(a4+exp(2i*pi*t)).^2)}, 'p');

% disc w/ shifted cot bump - no explosion
%a1=-0.01; a2=0.2; s = segment.radialfunc(N, {@(t) 1-a1*1i*cot((t + 1i*a2)/2),
%                    @(t) (1i*a1/2)./sin((t + 1i*a2)/2).^2});

% disc w/ gauss bump - shows explosion effect
%a1=0.1; a2=0.1; s = segment.radialfunc(N, {@(t) 1-a1*exp(-0.5*((t-pi)/a2).^2),
%                    @(t) (a1/a2^2)*(t-pi).*exp(-0.5*((t-pi)/a2).^2) });

%a=.3;b=.3;w=3; a1=0.1; a2=0.05; t0=3*pi/4; % smoothnonsym w/ gaussian bump @ t0
%s = segment.radialfunc(N, {@(q) 1 + a*cos(w*(q+b*cos(q))) - a1*exp(-0.5*((q-t0)/a2).^2), @(q) -a*sin(w*(q+b*cos(q))).*w.*(1-b*sin(q)) + (a1/a2^2)*(q-t0).*exp(-0.5*((q-t0)/a2).^2) });
% =====

%figure; s.plot;   % some plotting.........
if 0, taus = 0:0.01:2;    % animation
for i=1:numel(taus), t = [s.t; s.t(1)] + 1i*taus(i)/(2*pi);
  plot([s.x; s.x(1)], 'k+-'); hold on; plot(s.Z(t), 'b.-', 'linewidth', 2);
  hold off; axis equal; axis([-2 2 -2 2]); drawnow; pause(0.03); end
end

if 0, taus = 0:0.02:0.16; figure; plot([s.x; s.x(1)], 'k+-'); hold on;
for i=1:numel(taus), t = [s.t; s.t(1)] + 1i*taus(i)/(2*pi);
  plot(s.Z(t), 'b.-', 'linewidth', 2); end
axis equal;
end

% inverting Z(s) ...given segment and points z......
d = domain(s, 1);  % for inside test
if 0
 z = 0.5 + 0.0i + 0.1*exp(1i*2*pi*(0:0.01:1));  % circle of starting pts
else
  [xx yy] = meshgrid(-1.3:0.05:1.3); zg = xx(:)+1i*yy(:); ii=d.inside(zg);
  z = zg(ii);
end
tic; m = numel(z);
n=100;   % how many bdry pts to start from
imtmax = 0.5; % im part of t cutoff for iter
tol = 1e-14;
s.requadrature(n); si = nan(n,m); ni = nan(n,1);
if 1, % outer loop over starting pts, vectorize over z pts...
  for i=1:n
    x = s.t(i) + 0*z;        % starting pts
    kk = find(~isnan(x));    % start with all indices
    for j=1:100
      xold = x;
      %kkk = find(imag(x)>imtmax | imag(x)<0); % failed
      x(kk) = x(kk) + (z(kk) - s.Z(x(kk)))./s.Zp(x(kk));   % Newton iter
      %x(kkk) = xold(kkk); %freeze dead ones
      kk = find(abs(x-xold)>tol & imag(x)<imtmax & imag(x)>0);  % ignore Im t too large
      %jj = find(imag(x(kk))<0 | imag(x(kk))>imtmax);
      %x(kk(jj)) = real(x(kk(jj))); % bump up imag<0
      if numel(kk)==0, break; end          % group stopping criterion
    end
    ni(i) = j;   % save how many iters
    x(imag(x)>imtmax | imag(x)<=0) = nan;  % kill bad ones
    si(i,:) = x(:).';
  end
else  % outer loop over z pts, vectorize over starting pts... (no faster)
  ni = nan(m,1);
  for k=1:m
    x = s.t;
    for j=1:100
      xold = x;
      x = x + (z(k) - s.Z(x))./s.Zp(x);      % Newton iter
      if max(abs(x-xold))<1e-14, break; end  % group stopping criterion
    end
    ni(k) = j; si(:,k) = x(:);
  end
end
si = mod(real(si),1) + 1i*imag(si);  % wrap into [0,1] real part
[dummy,I] = sort(imag(si),1); sii = nan + 0*si;
for k=1:m, si(:,k) = si(I(:,k),k);   % as in help sort
  i = find(abs(diff([si(:,k); Inf]))>tol);
  sii(1:numel(i),k) = si(i,k); % keep unique vals
end
%sii = si(find(abs(diff([si; Inf*ones(1,m)],[],1))>tol));
toc
figure; subplot(1,2,1); s.requadrature(100);
s.plot; hold on; plot(real(z),imag(z), 'ko');
plot(s.Z(sii(find(~isnan(sii)))), 'r.'); axis(1.2*[-1 1 -1 1]);
subplot(1,2,2); plot(2*pi*si, 'r.'); hold on; plot(2*pi*s.t, 0*s.t, '.-');
axis([0 2*pi -2 2]); axis equal;

