%  function [ results ] = deflection_scan( x0, y0, xp0, yp0, gamma, Bmap1D, alpha )
%  %deflection_scan deflection for a given field map
%  

c = 2.998e8; % speed of light (m/s)

%   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)
%   outray is list of [z, x, y, xp, yp] at output position z

x0 = 0.644e-3;
y0 = 0;
xp0 = 0;
yp0 = 0;
gamma = 4000. / 0.511;
Bmap1D = readfieldmap( 'xleap2-wiggler-data/xLEAP-II-onaxis/xLEAP-II-1.dat' );
alpha = 0.003;

%  x0 = 0.5e-3;
%  y0 = 0;
%  xp0 = 0;
%  yp0 = 0;
%  gamma = 4000. / 0.511;
%  Bmap1D = readfieldmap( 'xleap2-wiggler-data/xLEAP-II-onaxis/xLEAP-II-2.dat' );
%  alpha = 0.003;

dx = 1e-6; % delta offset
dxp = 1e-6; % delta angle

mat = [];

ray0 = output_ray( x0, y0, xp0, yp0, gamma, Bmap1D, alpha );
dray = output_ray( x0+dx, y0, xp0, yp0, gamma, Bmap1D, alpha ) - ray0;
mat = [mat; [dray(2), dray(4), dray(3), dray(5)]];
dray = output_ray( x0, y0, xp0+dxp, yp0, gamma, Bmap1D, alpha ) - ray0;
mat = [mat; [dray(2), dray(4), dray(3), dray(5)]];
dray = output_ray( x0, y0+dx, xp0, yp0, gamma, Bmap1D, alpha ) - ray0;
mat = [mat; [dray(2), dray(4), dray(3), dray(5)]];
dray = output_ray( x0, y0, xp0, yp0+dxp, gamma, Bmap1D, alpha ) - ray0;
mat = [mat; [dray(2), dray(4), dray(3), dray(5)]];

mat = transpose(mat)