function [ results ] = trajectory( x0, y0, xp0, yp0, gamma, Bmap1D, alpha )
%getfield Get the trajectory
%   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 = 3e8; % speed of light (m/s)

% time steps
dz = map_zs(2) - map_zs(1);
Dz = max(map_zs) - min(map_zs);
t_max = Dz / c;
dt = dz / c;
nsteps = round(t_max / dt) * 0.01;
t_span = linspace(0, t_max, nsteps);
%  t_span = linspace(0, 10000*dt, 100);

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

% solve ODE
[t, results] = ode45(@(t, theta) ODE_func(t, theta, gamma, Bmap1D, alpha), t_span, theta0);

end