function [ results ] = field_integral_raw( x0, y0, Bmap, alpha )
%field_integral Get the trajectory from the field integrals
%   x0, y0 is initial transverse position
%   gamma0 is beam's lorentz factor
%   map_zs, map_Bxs, map_Bys is the fieldmap data
%   alpha is the field rolloff fraction at 1 mm: alpha = 1 - B(x=1mm) / B(x=0)

% field map
map_zs= Bmap(:,1);
map_Bxs= Bmap(:,2);
map_Bys= Bmap(:,3);

% time steps
dz = map_zs(2) - map_zs(1);
Dz = max(map_zs) - min(map_zs);
nsteps = round(Dz / dz);
z_span = linspace(min(map_zs), max(map_zs), nsteps);

% get field at params
[ Bx, By ] = getfield( z_span, x0, y0, Bmap, alpha );

% integrate equations of motion
I1x = cumsum(Bx * dz);
I1y = cumsum(By * dz);
I2x = cumsum(I1x * dz);
I2y = cumsum(I1y * dz);

results = [z_span; I2x; I2y; I1x; I1y]';

end