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
%   Bmap1D 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)

% 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);
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
%  opts = odeset('RelTol',1e-6,'AbsTol',1e-9);
opts = odeset('RelTol',1e-5,'AbsTol',1e-8);
[t, results] = ode45(@(t, theta) ODE_func(t, theta, gamma, Bmap1D, alpha), t_span, theta0, opts);

end