function p = plotct2D( phifile, N )

A = load(phifile);
A = unique( A, 'rows' );          % remove duplicates and sort

if nargin < 2
  Nx = sqrt(size(A,1));           % assume Nx = Ny
  Ny = Nx;
else
  Nx = N(1) + 1;
  Ny = N(2) + 1;
end

X = reshape( A(:,1), Ny, Nx );
Y = reshape( A(:,2), Ny, Nx );
Z = reshape( A(:,3), Ny, Nx );

p = surf(X, Y, Z);

%fill3(X, Y, Z, Z)
%shading interp
%axis off