function [ Bx, By ] = getfield( z, x, y, 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';
    display(Bmap)
    display([nz,nd])

if(iscell(Bmap)) % chaining field maps

    % figure out which wiggler
    [nq, nm] = size(Bmap);
    is = [1,1]; zs = [Bmap{1}{2}];
    for i = 2:nm
        is = [is, i, i]; zs = [zs, Bmap{i}{2}+[eps,0]];
    end
    iwigg = interp1(zs, is, z, 'nearest');
    
    % grab this wiggler
    Bm = Bmap{iwigg}{1};
    zs = Bmap{iwigg}{2};
    xs = Bmap{iwigg}{3};
    xoffset = interp1(zs, xs, z, 'linear', 0);

    % field map
    zs= Bm(:,1);
    Bxs= Bm(:,2);
    Bys= Bm(:,3);
    
%      lambda_wiggler = 0.00055;

    fieldrolloff = 1. - alpha * ((x-xoffset) / 0.001)^2;
%      fieldydependence = cosh(2*pi*y/lambda_wiggler);

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


else
    [nz, nd] = size(Bmap);

    if(nd == 3) % 1D map

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

        fieldrolloff = 1. - alpha * (x / 0.001)^2;
    %      fieldydependence = cosh(2*pi*y/lambda_wiggler);

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

    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 = 1e-3*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,2) * (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

end