%**************************************************************************
%This program analyzes arbitrary magnetic field data.  The trajectory of
%a beam of specified energy is calculated.  Using a specified radiation
%wavelength, the phase is calculated.
%
%The following structures are used:
%
%param, structure containing required program parameters
%	input_file_dir      %Input file directory
%   input_file_name     %Input file name
%   input_file_type     %Input file type, gives file format
%   gamma               %Lorentz factor
%   nom_zm_probe_at_dev_ctr  %Nominal z position in bench coordinates when the probe is at the device center (m)
%   nom_lambda_u        %Nominal undulator period (m)
%   num_poles           %Number of poles
%   nom_lambda_r        %Nominal radiation wavelength (m)
%
%field, structure containing magnetic field information
%   zm                  %Z positions at which the field measurements are made (m)
%   z                   %Z positions relative to the start position (m)
%   bx                  %Measured horizontal field (T)
%   by                  %Measured vertical field (T)
%   int_bx_dz           %First integral of Bx (Tm)
%   int_by_dz           %First integral of By (Tm)
%   int_int_bx_dz_dz    %Second integral of Bx (Tm^2)
%   int_int_by_dz_dz    %Second integral of By (Tm^2)
%   bx_filt             %Filtered measured horizontal field (T)
%   by_filt             %Filtered measured vertical field (T)
%   int_bx_dz_filt      %Filtered first integral of Bx (Tm)
%   int_by_dz_filt      %Filtered first integral of By (Tm)
%   int_int_bx_dz_dz_filt  %Filtered second integral of Bx (Tm^2)
%   int_int_by_dz_dz_filt  %Filtered second integral of By (Tm^2)
%
%traj, structure containing trajectory information
%   z                   %Z positions relative to the start position (m)
%   x                   %Trajectory x position (m)
%   x_p                 %First derivative dx/dz
%   x_pp                %Second derivative d(dx/dz)/dz (m^-1)
%   y                   %Trajectory y position (m)
%   y_p                 %First derivative dy/dz
%   y_pp                %Second derivative d(dy/dz)/dz (m^-1)
%   x_filt              %Filtered trajectory x position (m)
%   x_p_filt            %Filtered first derivative dx/dz
%   x_pp_filt           %Filtered second derivative d(dx/dz)/dz (m^-1)
%   y_filt              %Filtered trajectory y position (m)
%   y_p_filt            %Filtered first derivative dy/dz
%   y_pp_filt           %Filtered second derivative d(dy/dz)/dz (m^-1)
%
%rad, structure containing radiation information
%   z                   %Z positions relative to the start position (m)
%	slip                %Slippage, slip = 0 at z = 0 (m)
%   slip_drift_diff     %Slippage difference, measured minus slippage with x' = y' = 0, field free
%	phase               %Phase, (2 pi / lambda_r) * slip (rad)
%	phase_drift_diff    %Phase difference, (2 pi / lambda_r) * slip_drift_diff (rad)
%
%Calculations:
%
%X Trajectory:
%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
%x-trajectory
%gamma m d2x/dt2 = - q vz by
%d2x/dt2 = x'' vz^2
%gamma m x'' vz^2 = - q vz by
%x'' = - (q / gamma m vz) by
%x' = - (q / gamma m vz) int by dz
%x = - (q / gamma m vz) int int by dz dz
%
%Y Trajectory:
%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
%y-trajectory
%gamma m d2y/dt2 = q vz bx
%d2y/dt2 = y'' vz^2
%gamma m y'' vz^2 = q vz bx
%y'' = (q / gamma m vz) bx
%y' = (q / gamma m vz) int bx dz
%y = (q / gamma m vz) int int bx dz dz
%
%Electron speed along z axis:
%gamma = 1 / sqrt(1 - v^2/c^2)
%v^2 = c^2 (1 - 1/gamma^2)
%v^2 = vz^2 + vx^2 + vy^2 ~= vz^2 (1 + x'^2 + y'^2)
%vz^2 = v^2 / (1 + x'^2 + y'^2) ~= c^2 (1 - 1/gamma^2) (1 - x'^2 - y'^2)
%vz^2 ~= c^2 (1 - 1/gamma^2 - x'^2 - y'^2)
%vz ~= c (1 - (1/2) (1/gamma^2) - (1/2) x'^2 - (1/2) y'^2)
%
%Slippage, distance between electrons and light:
%d_slip = (c - vz) d_t = (1 - vz/c) d_z
%d_slip = ((1/2) (1/gamma^2) + (1/2) x'^2 + (1/2) y'^2) d_z
%slip = (1/2) int (1/gamma^2 + x'^2 + y'^2) d_z
%
%Phase:
%phase = (2 pi / lambda_r) slip
%
%Zachary Wolf, 2/10/06

function [] = pid_traj_anal()

%Set up program parameters using the user interface
%Create the updated pid_traj_anal_param.m file which is run to set the parameter values
%pid_traj_anal_ui;

%Clear the workspace
clear

%Initialize global physical constants used throughout the program
global Q_ELEC M_ELEC C_LIGHT;
Q_ELEC = -1.60217733e-19;    %Electron charge (C)
M_ELEC = 9.1093897e-31;      %Electron mass (kg)
C_LIGHT = 2.99792458e8;      %Speed of light (m/s)

%Get the program parameters
[param] = pid_traj_anal_get_param();

%Get the measured data
[header, zm, bx, by] = pid_traj_anal_get_meas(param);

%Open output files, write headers
[dat_file, plt_file] = pid_traj_anal_init(param, header);

%Apply a field correction, if required
if param.corr_required == true
	[bx, by] = pid_traj_anal_field_corr(param, dat_file, plt_file, zm, bx, by);
end

%Perform a magnetic field analysis
[field] = pid_traj_anal_field(param, dat_file, plt_file, zm, bx, by);

%Perform the trajectory analysis
[traj] = pid_traj_anal_traj(param, dat_file, plt_file, field);

%Perform the phase analysis
[rad] = pid_traj_anal_phase(param, dat_file, plt_file, field, traj);

%Output the results to the input file directory
mat_file = strcat(param.input_file_dir, 'pid_traj_anal.mat');
save(mat_file, 'param', 'header', 'field', 'traj', 'rad');

%Done
