function [ results ] = field_integral( x0, y0, xp0, yp0, gamma, Bmap1D, 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)

c = 2.998e8; % speed of light (m/s)
m = 9.11e-31; % electtron mass (kg)
q = -1.6e-19; % electron charge (C)

% field map
map_zs= Bmap1D(:,1);
map_Bxs= Bmap1D(:,2);
map_Bys= Bmap1D(:,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);

% initial coords
pos0 = [x0, y0, min(map_zs)];
%  vel0 = [0, 0, c*sqrt(1. - gamma^-2)];
v = c*sqrt(1. - gamma^-2);
vx = c * xp0;
vy = c * yp0;
vz = sqrt(v^2 - vx^2 - vy^2);
vel0 = [vx, vy, vz];

% get field at params
[ Bx, By ] = getfield( z_span, 0, 0, Bmap1D, alpha );

% integrate equations of motion
qbymgamma = q / m / gamma;
dvxdz = qbymgamma * - By;
dvydz = qbymgamma * Bx;
vx = vel0(1) + cumsum(dvxdz * dz);
vy = vel0(2) + cumsum(dvydz * dz);
x = x0 + cumsum(vx * dz / c);
y = y0 + cumsum(vy * dz / c);

results = [z_span; x; y; vx/c; vy/c]';

end