function [ Bx, By ] = getfield( z, x, Bmap, alpha )
%getfield Get the field at a point
%   z is longitudinal position queried
%   x is transverse position queried
%   Bmap is the fieldmap data (1D or 2D)
%   alpha is the field rolloff fraction at 1 mm: alpha = 1 - B(x=1mm) / B(x=0)
    
%  interpmethod = 'spline';
%  interpmethod = 'linear';
interpmethod = 'nearest';

[nz, nd] = size(Bmap);

if(nd == 3) % 1D map

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

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

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

else % 2D field map

    % 2d interp isn't working...
    
%      zs= reshape(Bmap(:,1),[nz/11,11]);
%      Bxs= reshape(Bmap(:,2),[nz/11,11]);
%      Bys= reshape(Bmap(:,3),[nz/11,11]);
%      xs= reshape(Bmap(:,4),[nz/11,11]);
%  
%      Bx = interp2(zs, xs, Bxs, z, x, 'spline', 0);
%      By = interp2(zs, xs, Bys, z, x, 'spline', 0);

    % ... so just average two nearby ( I checked that the zs are identical for each scan)

    % field map
    
    xs = linspace(-1,1,11);
    xi = interp1(xs, 1:11, x, 'linear');
    xifloor = floor(xi);
    
    ilo = 1+(xifloor-1)*(nz/11); ilop1 = 1+(xifloor+1-1)*(nz/11);
    ihi = xifloor*(nz/11); ihip1 = (xifloor+1)*(nz/11);
    zs = Bmap(ilo:ihi,1);
    Bxs = Bmap(ilo:ihi,2) * (xifloor + 1 - xi) + Bmap(ilop1:ihip1,3) * (xi - xifloor);
    Bys = Bmap(ilo:ihi,3) * (xifloor + 1 - xi) + Bmap(ilop1:ihip1,3) * (xi - xifloor);
    
    Bx = interp1(zs, Bxs, z, interpmethod, 0);
    By = interp1(zs, Bys, z, interpmethod, 0);

end

end