%**************************************************************************
%This function calculates the velocity components of a particle in an EPU.
%The calculation is done to first order.  Second order effects cannot be
%included since the gradient of the field must be known for second order
%calculations.  See LCLS-TN-15-3.
%
%Input:
%gamma, E / m c^2
%int_bx_dz, first integral of Bx (Tm)
%int_by_dz, first integral of By (Tm)
%
%Output:
%vx, x component of the velocity (m/s)
%vy, y component of the velocity (m/s)
%vz, z component of the velocity (m/s)
%
%Zachary Wolf
%3/23/15

%Method
%d/dt(gamma m V) = q V x B
%gamma m d/dt(V) = q (-vz by X + vz bx Y)
%gamma m d/dt(vx) = - q vz by
%gamma m d/dt(vy) = q vz bx
%d/dt(vx) = - (q / gamma m) vz by
%d/dt(vy) = (q / gamma m) vz bx
%d/dt = (dz/dt) d/dz = vz d/dz
%d/dz(vx) = - (q / gamma m) by
%d/dz(vy) = (q / gamma m) bx
%vx = - (q / gamma m) (int by dz)
%vy = (q / gamma m) (int bx dz)
%gamma^2 = 1 / [1 - (vx^2 + vy^2 +vz^2) / c^2]
%(1 / gamma^2) = 1 - (vx^2 + vy^2 +vz^2) / c^2
%vx^2 + vy^2 + vz^2 = c^2 (1 - (1 / gamma^2))
%vz^2 = c^2 (1 - (1 / gamma^2)) - vx^2 - vy^2
%vz = sqrt[c^2 (1 - (1 / gamma^2)) - vx^2 - vy^2]

function [vx, vy, vz] = traj_calc_vel(gamma, int_bx_dz, int_by_dz)

%Physical constants
global Q_ELEC;               %Electron charge (C)
global M_ELEC;               %Electron mass (kg)
global C_LIGHT;              %Speed of light (m/s)

%Calculate vx
vx = - (Q_ELEC / (gamma * M_ELEC)) * int_by_dz;

%Calculate vy
vy = (Q_ELEC / (gamma * M_ELEC)) * int_bx_dz;

%Calculate vz
vz = sqrt(C_LIGHT^2 * (1 - (1 / gamma^2)) - vx.^2 - vy.^2);
    
%Done
