function [u, v] = rect2grid(z, zrect, theCorners, theSize)

% rect2grid -- Orthogonal grid from RECT result.
%  [u, v] = rect2grid(z, zrect, theCorners, theSize) produces
%   a curvilinear orthogonal grid by interpolating the complex
%   contour z, using zrect, the result of applying the conformal
%   mapper RECT to z for theCorners (indices).  If zrect is empty,
%   the RECT routine is called to compute it.  If zrect is a scalar,
%   that number of RECT iterations will be performed.  The size of
%   the u and v output grids (matrices), including the perimeter
%   is determined by theSize; it includes the perimeter.
%  rect2grid(nPoints) demonstrates itself with a random z contour
%   of nPoints (default = 20).
 
% Copyright (C) 1998 Dr. Charles R. Denham, ZYDECO.
%  All Rights Reserved.
%   Disclosure without explicit written consent from the
%    copyright owner does not constitute publication.
 
% Version of 21-Oct-1998 20:50:16.
% Revised    26-Oct-1998 14:29:11.

if nargin < 1, help(mfilename), z = 'demo'; end
if isequal(z, 'demo'), z = 20; end
if ischar(z), z = eval(z); end

if length(z) == 1
	n = z;
	z = rand(n, 1) + sqrt(-1)*rand(n, 1);
	z = z - mean(z);
	a = angle(z);
	[a, i] = sort(a);
	z = (1 + 0.1 * rand(size(a))) .* exp(sqrt(-1)*a);
	[ignore, nn] = sort(rand(1, length(z)-1));
	theCorners = sort([1 nn(1:3)+1]);
	theSize = [n n];
	zrect = [];
	tic
	[uu, vv] = rect2grid(z, zrect, theCorners, ceil(theSize/2));
	disp([' ## Elapsed time: ' num2str(toc) ' seconds.'])
	if ~isempty(uu)
		x = real(z); y = imag(z);
		x = [x; x(1)]; y = [y; y(1)];  % Make closed curve.
		uu1 = uu(:, 2:end-1);   % Trim the grid.
		vv1 = vv(:, 2:end-1);
		uu2 = uu(2:end-1, :).';
		vv2 = vv(2:end-1, :).';
		h = plot(uu1, vv1, 'g-', uu2, vv2, 'b-');
		hold on
		plot(x, y, 'r-', x(theCorners), y(theCorners), 'ro')
		hold off
		xlabel('x'), ylabel('y')
		theCommand = [mfilename ' ( ' int2str(n) ' )'];
		title(theCommand)
		set(gcf, 'ButtonDownFcn', theCommand)
		figure(gcf)
		axis equal
		zoomsafe
	end
	return
end

% Apply RECT until the straightness of the result
%  deviates from 1.0 by no more than 0.1 percent.

if length(zrect) < 2
	if length(zrect) == 1
		ntimes = zrect;
	else
		ntimes = ceil(sqrt(length(z)));
	end
	x = real(z); y = imag(z);
	zrect = x(:) + sqrt(-1).*y(:);
	for i = 1:ntimes
		[zrect, straight] = rect(zrect, 1, theCorners);
		if norm(1-straight) <= 0.001, break, end
	end
	if norm(1-straight) > 0.001
		disp([' ## rect2grid: Conformal mapping not successful'])
		disp(['               after ' int2str(ntimes) ' iterations.'])
		u = []; v = [];
		return
	end
end

if length(theSize) == 1
	theSize = theSize * [1 1];
end

theSize = theSize;

m = theSize(1); n = theSize(2);

% Get indices of matrix perimeter.

temp = zeros(theSize);
temp(:) = 1:prod(theSize);

ind = [];
ind = [ind; temp(1:m-1, 1)];
ind = [ind; temp(m, 1:n-1).'];
ind = [ind; temp(m:-1:2, n)];
ind = [ind; temp(1, n:-1:1).'];

% Interpolate around the boundary.

a = [0; cumsum(abs(diff([z(:); z(1)])))];
a = a - min(a);
a = a / max(a);

c = [theCorners length(a)];
d = cumsum([1 m-1 n-1 m-1 n-1]);

ai = zeros(size(a));
for i = 1:4
	cc = c(i):c(i+1); cc = cc - min(cc); cc = cc / max(cc);
	dd = d(i):d(i+1); dd = dd - min(dd); dd = dd / max(dd);
	ai(d(i):d(i+1)) = interp1(cc, a(c(i):c(i+1)), dd).';
end

x = real(z);
y = imag(z);
xi = interp1(a(:), [x(:); x(1)], ai(:));
yi = interp1(a(:), [y(:); y(1)], ai(:));

% Sprinkle interpolated values along the perimeter.

u = zeros(theSize); v = zeros(theSize);
u(ind) = xi; v(ind) = yi;

% Sample intervals: how to incorporate?

dx = abs(c(3)-c(2)) / n;
dy = abs(c(2)-c(1)) / m;
if (1), dx = 1; dy = 1; end

% Solve Laplace's equation inside the boundary.

u = fps(u);
v = fps(v);