%  function [ results ] = deflection_scan( x0, y0, xp0, yp0, gamma, Bmap, 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;
%  Bmap = 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;
%  Bmap = readfieldmap( 'xleap2-wiggler-data/xLEAP-II-onaxis/xLEAP-II-2.dat' );
%  alpha = 0.003;

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

mat = [];

% 2 by 2
ray0 = output_ray( x0, y0, xp0, yp0, gamma, Bmap, alpha );
dray = output_ray( x0+dx, y0, xp0, yp0, gamma, Bmap, alpha ) - ray0;
mat = [mat; [dray(2), dray(4)] / dx];
drayp = output_ray( x0, y0, xp0+dxp, yp0, gamma, Bmap, alpha ) - ray0;
mat = [mat; [drayp(2), drayp(4)] / dxp];

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

mat = transpose(mat)