function [ results ] = trajectory( x0, y0, gamma, map_zs, map_Bxs, map_Bys, 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);
%  t_span = linspace(0, t_max, nsteps);
t_span = linspace(0, 10000*dt, 100);

% initial coords
pos0 = [x0, y0, min(map_zs)];
vel0 = [0., 0., c*sqrt(1. - gamma^-2)];
theta0 = [pos0, vel0];

% solve ODE
%  [t, results] = ode45(@(t, theta) ODE_func(t, theta, map_zs, map_Bxs, map_Bys, alpha), t_span, theta0);

% integrate equations of motion
dvxdz = 


end