try
    uzs = unique(Bmap(:,1));
    nzs = numel(uzs);
    uxs = unique(Bmap(:,4));
    nxs = numel(uxs);
    zs = reshape(Bmap(:,1),[nzs,nxs])';
    xs = reshape(Bmap(:,4),[nzs,nxs])';
    Bxs = reshape(Bmap(:,2),[nzs,nxs])';
    Bys = reshape(Bmap(:,3),[nzs,nxs])';
catch
    fprintf('plot_2dfields_contours.m - WARNING: non-square grid?\n');
    [nzs,~] = size(Bmap);
    uxs = unique(Bmap(:,4));
    nxs = numel(uxs);
    nzs = floor(nzs/nxs);
    zs = reshape(Bmap(:,1),[nzs,nxs])';
    xs = reshape(Bmap(:,4),[nzs,nxs])';
    Bxs = reshape(Bmap(:,2),[nzs,nxs])';
    Bys = reshape(Bmap(:,3),[nzs,nxs])';
end

figure
contourf(zs, 1e3*xs, Bxs);
xlabel('z position (m)'); ylabel('x position (mm)');
colorbar(); title('Bx field (T)')
enhance_plot('arial',20,2,8)

figure
contourf(zs, 1e3*xs, Bys);
xlabel('z position (m)'); ylabel('x position (mm)');
colorbar(); title('By field (T)')
enhance_plot('arial',20,2,8)