function [ Bx, By ] = getfield( z, x, Bmap1D, alpha )
%getfield Get the field at a point
%   z is longitudinal position queried
%   x is transverse position queried
%   [zs, Bxs, Bys] is the fieldmap data
%   alpha is the field rolloff fraction at 1 mm: alpha = 1 - B(x=1mm) / B(x=0)

% field map
zs= Bmap1D(:,1);
Bxs= Bmap1D(:,2);
Bys= Bmap1D(:,3);

fieldrolloff = 1. - alpha * (x / 0.001)^2;

Bx = interp1(zs, Bxs, z, 'spline', 0) * fieldrolloff;
By = interp1(zs, Bys, z, 'spline', 0) * fieldrolloff;

end