%**************************************************************************
%This function calculates the time derivative of the velocity components
%for a particle in the undulator.
%
%Input:
%z, z position along the undulator (m)
%vx, x component of the velocity (m/s)
%vy, y component of the velocity (m/s)
%vz, z component of the velocity (m/s)
%
%Output:
%ax, x component of the acceleration (m/s/s)
%ay, y component of the acceleration (m/s/s)
%az, z component of the acceleration (m/s/s)
%
%Zachary Wolf
%7/10/10

function [ax, ay, az] = rad_calc_acc(z, vx, vy, vz)

%Parameter
INTERP_FACTOR = 8;
  
%Check
if (length(z) ~= length(vx)) || (length(z) ~= length(vy)) || (length(z) ~= length(vz))
    msgbox('rad_calc_ang_spect: arrays have different lengths', 'Error');
    return;
end
if length(z) < 10
    msgbox('rad_calc_ang_spect: problem with array sizes', 'Error');
    return;
end

%Make a dense set of values for the calculation
zz = z(1) : (z(2) - z(1)) / INTERP_FACTOR : z(end);
vxx = spline(z, vx, zz);
vyy = spline(z, vy, zz);
vzz = spline(z, vz, zz);

%Initialize derivatives wrt z
dvxxdzz = zeros(1, length(zz));
dvyydzz = zeros(1, length(zz));
dvzzdzz = zeros(1, length(zz));

%Calculate first values
dvxxdzz(1) = (vxx(2) - vxx(1)) / (zz(2) - zz(1));
dvyydzz(1) = (vyy(2) - vyy(1)) / (zz(2) - zz(1));
dvzzdzz(1) = (vzz(2) - vzz(1)) / (zz(2) - zz(1));

%Calculate last values
dvxxdzz(end) = (vxx(end) - vxx(end - 1)) / (zz(end) - zz(end - 1));
dvyydzz(end) = (vyy(end) - vyy(end - 1)) / (zz(end) - zz(end - 1));
dvzzdzz(end) = (vzz(end) - vzz(end - 1)) / (zz(end) - zz(end - 1));

%Calculate derivatives wrt z
for ii = 2 : length(zz) - 1
    dvxxdzz = (vxx (ii + 1) - vxx(ii - 1)) / (zz(ii + 1) - zz(ii - 1));
    dvyydzz = (vyy (ii + 1) - vyy(ii - 1)) / (zz(ii + 1) - zz(ii - 1));
    dvzzdzz = (vzz (ii + 1) - vzz(ii - 1)) / (zz(ii + 1) - zz(ii - 1));
end
    
%Calculate derivatives wrt time
%df/dt = df/dz dz/dt
dvxxdt = dvxxdzz .* vzz;
dvyydt = dvyydzz .* vzz;
dvzzdt = dvzzdzz .* vzz;

%Return the acceleration components at the input z values
ax = spline(zz, dvxxdt, z);
ay = spline(zz, dvyydt, z);
az = spline(zz, dvzzdt, z);
    
%Done
