function slice = load_comsol_slice(dataFile, z0, tol)
%LOAD_COMSOL_SLICE  Load COMSOL exported slice data.
% Supported formats (header line starting with % is ignored):
%   1) x  y  normE
%   2) x  y  z  normE
%   3) x  y  z  Ex_re  Ex_im  Ey_re  Ey_im  Ez_re  Ez_im  (normE computed)
if nargin < 3
    tol = 1e-3;
end
assert(isfile(dataFile), 'COMSOL data not found: %s', dataFile);

raw = readmatrix(dataFile, 'FileType', 'text', 'CommentStyle', '%');
if size(raw, 2) >= 9
    Ex = raw(:, 4) + 1i * raw(:, 5);
    Ey = raw(:, 6) + 1i * raw(:, 7);
    Ez = raw(:, 8) + 1i * raw(:, 9);
    val = sqrt(abs(Ex).^2 + abs(Ey).^2 + abs(Ez).^2);
    x = raw(:, 1);
    y = raw(:, 2);
    z = raw(:, 3);
elseif size(raw, 2) == 4
    x = raw(:, 1);
    y = raw(:, 2);
    z = raw(:, 3);
    val = raw(:, 4);
elseif size(raw, 2) == 3
    x = raw(:, 1);
    y = raw(:, 2);
    z = z0 * ones(size(x));
    val = raw(:, 3);
else
    error('Unsupported COMSOL file format (need 3 or 4+ columns).');
end

if nargin >= 2 && ~isempty(z0)
    mask = abs(z - z0) <= tol;
    if ~any(mask)
        warning('No points with z~=%g in COMSOL file; using all points.', z0);
        mask = true(size(x));
    end
    x = x(mask);
    y = y(mask);
    val = val(mask);
end

slice.x = x(:);
slice.y = y(:);
slice.val = val(:);
slice.elem = [];
slice.planeAxes = [1 2];
slice.z0 = z0;
end