function [ results ] = trajectory( x0, y0, xp0, yp0, gamma, Bmap, alpha )
%getfield Get the trajectory
%   x0, y0 is initial transverse position
%   gamma0 is beam's lorentz factor
%   Bmap 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)

if(iscell(Bmap))
    [nm, na] = size(Bmap);
    % get total length
    zs = [];
    for i = 1:nm
        zs = [zs, Bmap{i}{2}];
    end
    Dz = max(zs) - min(zs);
    z0 = min(zs);
    % get increment length
    Bm = Bmap{1}{1};
    dz = Bm(2,1) - Bm(1,1);
else
    % field map
    map_zs= Bmap(:,1);
    dz = map_zs(2) - map_zs(1);
    Dz = map_zs(end) - map_zs(1);
    z0 = map_zs(1);
end

% time steps
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, z0];
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, Bmap, alpha), t_span, theta0, opts);

end