%**************************************************************************
%This function performs an analysis related to undulator radiation.
%The slippage, the distance between the light and the electrons, is calculated.
%Undulator parameters, such as the K value, are calculated.
%The phase is analyzed, including rms phase errors and total phase through
%the device.
%
%Input:
%param, parameter structure
%dat_file, output data file name
%plt_file, output plot file name
%field, field structure
%traj, trajectory structure
%
%Output:
%phas, structure containing phase information
%	z                   %Translated Z, from traj structure (m)
%	z_peak              %z at the field peaks (m)
%	z_peak_no_ends      %z at the field peaks excluding the ends (m)
%   lambda_u            %Measured average undulator period (m)
%	slip                %Slippage, slip = 0 at z = 0 (m)
%	slip_peak           %Slippage at the field peaks (m)
%	slip_peak_no_ends	%Slippage at the field peaks excluding the ends (m)
%	slope               %Slope of the slip_peak vs z_peak (no ends) linear fit
%   slip_resid_peak_no_ends  %Residuals from linear fit of slip_peak vs z_peak (m)
%	k_eff               %K value for undulator
%	kx_eff              %K value for undulator using only the x trajectory
%	ky_eff              %K value for undulator using only the y trajectory
%	b_eff               %Effective field strength in undulator (T)
%	lambda_r            %Radiation wavelength (m)
%	k_eff_sc            %K value for undulator using the slope corrected trajectories
%	b_eff_sc            %Effective field strength in undulator using slope corrected trajectories (T)
%	lambda_r_sc         %Radiation wavelength using slope corrected trajectories (m)
%	phase               %Phase, phase = 0 at cell boundary (rad)
%	phase_peak          %Phase at the field peaks (rad)
%	phase_peak_no_ends	%Phase at the field peaks excluding the ends (rad)
%	phase_error_peak_no_ends	 %Phase errors at the field peaks excluding the ends (rad)
%	phase_error_peak_no_ends_ave %Average value of phase errors at the field peaks excluding the ends (rad)
%	phase_error_peak_no_ends_rms %RMS variation of phase errors at the field peaks excluding the ends (rad)
%	phase_sc            %Phase, phase = 0 at cell boundary, use slope corrected trajectories (rad)
%	phase_peak_sc       %Phase at the field peaks, use slope corrected trajectories (rad)
%	phase_peak_no_ends_sc  %Phase at the field peaks excluding the ends, use slope corrected trajectories (rad)
%	phase_error_peak_no_ends_sc	 %Phase errors at the field peaks excluding the ends, use slope corrected trajectories (rad)
%	phase_error_peak_no_ends_ave_sc  %Average value of phase errors at the field peaks excluding the ends, use slope corrected trajectories (rad)
%	phase_error_peak_no_ends_rms_sc  %RMS variation of phase errors at the field peaks excluding the ends, use slope corrected trajectories (rad)
%	phase_cell          %Total phase change between the light and the electrons in the cell (rad)
%	phase_nom           %Phase relative to nominal radiation wave, phase = 0 at cell boundary (rad)
%	phase_nom_peak      %Phase relative to nominal radiation wave at the field peaks (rad)
%	phase_nom_peak_no_ends  %Phase relative to nominal radiation wave at the field peaks excluding the ends (rad)
%	phase_nom_error_peak_no_ends     %Phase errors relative to nominal radiation wave (rad)
%	phase_nom_error_peak_no_ends_ave %Average value of phase errors relative to nominal radiation wave (rad)
%	phase_nom_error_peak_no_ends_rms %RMS variation of phase errors relative to nominal radiation wave (rad)
%	phase_nom_error_peak_no_ends_enter_ave %Average value of phase errors relative to nominal radiation wave at the entrance (rad)
%	phase_nom_error_peak_no_ends_enter_rms %RMS variation of phase errors relative to nominal radiation wave at the entrance (rad)
%	phase_nom_error_peak_no_ends_exit_ave  %Average value of phase to exit errors relative to nominal radiation wave at the exit (rad)
%	phase_nom_error_peak_no_ends_exit_rms  %RMS variation of phase to exit errors relative to nominal radiation wave at the exit (rad)
%	phase_nom_cell      %Total phase change between the nominal radiation wave and the electrons in the cell (rad)
%   phase_integral      %Int_-inf^z [(Int_-inf^z'' Bx(z') dz')^2 + (Int_-inf^z'' By(z') dz')^2] dz'' (T^2m^3)
%   phase_integral_end  %Int_-inf^+inf [(Int_-inf^z'' Bx(z') dz')^2 + (Int_-inf^z'' By(z') dz')^2] dz'' (T^2m^3)
%   phase_integral_cell %Int_z_enter^z_exit [(Int_-inf^z'' Bx(z') dz')^2 + (Int_-inf^z'' By(z') dz')^2] dz'' (T^2m^3)
%   l_eff               %Effective length of undulator contributing radiation (m)
%
%Zachary Wolf
%7/8/10

function [phas] = epu_anal_phase(param, dat_file, plt_file, field, traj)

%Extract required information
gamma = param.gamma;
num_end_poles = param.num_end_poles;
nom_lambda_r = param.nom_lambda_r;
cell_length = param.cell_length;
nom_phase_cell = param.nom_phase_cell;
num_pole_stat = param.num_pole_stat;
lambda_u = field.lambda_u;
z_peak = field.z_peak;
z_peak_no_ends = field.z_peak_no_ends;
int_bx_dz = field.int_bx_dz;
int_by_dz = field.int_by_dz;
z = traj.z;
x_p = traj.x_p;
y_p = traj.y_p;
x_p_slope_cor = traj.x_p_slope_cor;
y_p_slope_cor = traj.y_p_slope_cor;
vz = traj.vz;

%Slippage -----

%Calculate the slippage as a function of z
slip_vz = slip_calc_slip_vz(z, vz);
slip_slope_o1 = slip_calc_slip_slope_o1(gamma, z, x_p, y_p);  %trajectory slope, order 1
slip_slope = slip_calc_slip_slope(gamma, z, x_p, y_p);  %calculated from trajectory slope
slip = slip_vz;

%Calculate the slippage using x_p and y_p alone
slip_xp = slip_calc_slip_slope(gamma, z, x_p, zeros(1, length(y_p)));
slip_yp = slip_calc_slip_slope(gamma, z, zeros(1, length(x_p)), y_p);

%Calculate the slippage using slope corrected trajectories
slip_sc = slip_calc_slip_slope(gamma, z, x_p_slope_cor, y_p_slope_cor);

%Debug
% slip_diff = slip_vz - slip_slope_o1;
% slip_plot_slip_diff(plt_file, z, slip_diff);
% slip_diff = slip_vz - slip_slope;
% slip_plot_slip_diff(plt_file, z, slip_diff);
% slip_diff = slip_slope - slip_slope_o1;
% slip_plot_slip_diff(plt_file, z, slip_diff);

%Find the slippage at the field peaks
slip_peak = slip_calc_slip_at_peaks(z, slip, z_peak);
slip_xp_peak = slip_calc_slip_at_peaks(z, slip_xp, z_peak);
slip_yp_peak = slip_calc_slip_at_peaks(z, slip_yp, z_peak);
slip_peak_sc = slip_calc_slip_at_peaks(z, slip_sc, z_peak);

%Perform a linear fit to the slippage at the field peaks, neglect ends
z_peak_no_ends = z_peak(num_end_poles + 1 : length(z_peak) - num_end_poles);
slip_peak_no_ends = slip_peak(num_end_poles + 1 : length(slip_peak) - num_end_poles);
[slope, intercept, slip_resid_peak_no_ends] = slip_linfit_slip(z_peak_no_ends, slip_peak_no_ends);
slip_xp_peak_no_ends = slip_xp_peak(num_end_poles + 1 : length(slip_xp_peak) - num_end_poles);
[slope_xp, intercept, slip_xp_resid_peak_no_ends] = slip_linfit_slip(z_peak_no_ends, slip_xp_peak_no_ends);
slip_yp_peak_no_ends = slip_yp_peak(num_end_poles + 1 : length(slip_yp_peak) - num_end_poles);
[slope_yp, intercept, slip_yp_resid_peak_no_ends] = slip_linfit_slip(z_peak_no_ends, slip_yp_peak_no_ends);
slip_peak_no_ends_sc = slip_peak_sc(num_end_poles + 1 : length(slip_peak_sc) - num_end_poles);
[slope_sc, intercept, slip_resid_peak_no_ends_sc] = slip_linfit_slip(z_peak_no_ends, slip_peak_no_ends_sc);

%Plot the slippage
%slip_plot_slip(plt_file, z, slip);
%slip_plot_slip_pole(plt_file, z, slip, z_peak);

%Plot the slippage residuals
%slip_plot_slip_resid(plt_file, z_peak_no_ends, slip_resid_peak_no_ends);
%slip_plot_slip_resid_pole(plt_file, z_peak_no_ends, slip_resid_peak_no_ends, z_peak);

%Keff, Beff, Lambda_rad -----

%Calculate Keff
k_eff = k_calc_keff(gamma, slope);
kx_eff = k_calc_keff(gamma, slope_xp);
ky_eff = k_calc_keff(gamma, slope_yp);
k_eff_sc = k_calc_keff(gamma, slope_sc);

%Calculate Beff
b_eff = k_calc_beff(lambda_u, k_eff);
b_eff_sc = k_calc_beff(lambda_u, k_eff_sc);

%Calculate the radiation wavelength
lambda_r = k_calc_lambdar(lambda_u, slope);
lambda_r_sc = k_calc_lambdar(lambda_u, slope_sc);  %Slope corrected trajectories

%Write the undulator parameters to the data file
k_dat_keff_beff_lambdar(dat_file, k_eff, b_eff, lambda_r);
k_dat_keff_beff_lambdar_slope_corr_traj(dat_file, k_eff_sc, b_eff_sc, lambda_r_sc);

%Phase -----

%Calculate the phase as a function of z, phase(z(1)) = 0 ref
phase = phase_calc_phase(lambda_r, slip);
phase_sc = phase_calc_phase(lambda_r_sc, slip_sc);

%Calculate the phase at the field peaks
phase_peak = phase_calc_phase(lambda_r, slip_peak);
phase_peak_sc = phase_calc_phase(lambda_r_sc, slip_peak_sc);

%Calculate the phase at the field peaks, neglect the ends
phase_peak_no_ends = phase_calc_phase(lambda_r, slip_peak_no_ends);
phase_peak_no_ends_sc = phase_calc_phase(lambda_r_sc, slip_peak_no_ends_sc);

%Calculate the phase at the cell entrance
phase_cell_entrance = phase_calc_phase_cell_entrance(lambda_r, gamma, cell_length, z, phase);

%Calculate the phase at the cell exit
phase_cell_exit = phase_calc_phase_cell_exit(lambda_r, gamma, cell_length, z, phase);

%Make the cell entrance phase = 0 reference
phase = phase - phase_cell_entrance;
phase_peak = phase_peak - phase_cell_entrance;
phase_peak_no_ends = phase_peak_no_ends - phase_cell_entrance;
phase_cell_exit = phase_cell_exit  - phase_cell_entrance;
phase_cell_entrance = 0;  %Keep at end

%Calculate the phase change through the cell
phase_cell = phase_cell_exit;

%Calculate the phase errors at the field peaks, neglect ends
phase_error_peak_no_ends = phase_calc_phase_error(lambda_r, slip_resid_peak_no_ends);
phase_error_peak_no_ends_sc = phase_calc_phase_error(lambda_r_sc, slip_resid_peak_no_ends_sc);

%Calculate statistics of the phase errors
[phase_error_peak_no_ends_ave, phase_error_peak_no_ends_rms] = phase_calc_phase_error_stats(phase_error_peak_no_ends);
[phase_error_peak_no_ends_ave_sc, phase_error_peak_no_ends_rms_sc] = phase_calc_phase_error_stats(phase_error_peak_no_ends_sc);

%Plot the phase
phase_plot_phase(plt_file, z, phase);

%Plot the phase at the field peaks
%phase_plot_phase_at_peaks(plt_file, z_peak, phase_peak);

%Plot the phase errors
%phase_plot_phase_errors(plt_file, z_peak_no_ends, phase_error_peak_no_ends);
phase_plot_phase_errors_pole(plt_file, z_peak_no_ends, phase_error_peak_no_ends, z_peak);
phase_plot_phase_errors_pole_traj_slope_corr(plt_file, z_peak_no_ends, phase_error_peak_no_ends_sc, z_peak);


%Phase To Nominal Wave -----

%Use the calculated lambda_r if nom_lambda_r = 0
if nom_lambda_r == 0
    nom_lambda_r = lambda_r;
end

%Calculate the phase relative to the nominal radiation wave as a function of z, phase(z(1)) = 0 ref
phase_nom = phase_calc_phase(nom_lambda_r, slip);

%Calculate the phase relative to the nominal radiation wave at the field peaks
phase_nom_peak = phase_calc_phase(nom_lambda_r, slip_peak);

%Calculate the phase relative to the nominal radiation wave at the field peaks, neglect the ends
phase_nom_peak_no_ends = phase_calc_phase(nom_lambda_r, slip_peak_no_ends);

%Calculate the phase relative to the nominal radiation wave at the cell entrance
phase_nom_cell_entrance = phase_calc_phase_cell_entrance(nom_lambda_r, gamma, cell_length, z, phase_nom);

%Calculate the phase relative to the nominal radiation wave at the cell exit
phase_nom_cell_exit = phase_calc_phase_cell_exit(nom_lambda_r, gamma, cell_length, z, phase_nom);

%Make the cell entrance phase_nom = 0 reference
phase_nom = phase_nom - phase_nom_cell_entrance;
phase_nom_peak = phase_nom_peak - phase_nom_cell_entrance;
phase_nom_peak_no_ends = phase_nom_peak_no_ends - phase_nom_cell_entrance;
phase_nom_cell_exit = phase_nom_cell_exit  - phase_nom_cell_entrance;
phase_nom_cell_entrance = 0;  %Keep at end

%Calculate the phase change through the cell
phase_nom_cell = phase_nom_cell_exit;

%Calculate the phase error relative to the nominal radiation wave in the undulator
calc_phase_nom_peak_no_ends = phase_calc_phase_nom_peak_no_ends(param, length(phase_nom_peak_no_ends));
phase_nom_error_peak_no_ends = phase_nom_peak_no_ends - calc_phase_nom_peak_no_ends;

%Calculate the phase to the cell exit
phase_nom_to_exit_peak_no_ends = phase_nom_cell - phase_nom_peak_no_ends;

%Calculate the error in the phase to the cell exit
calc_phase_nom_to_exit_peak_no_ends = nom_phase_cell - calc_phase_nom_peak_no_ends;
phase_nom_error_to_exit_peak_no_ends = phase_nom_to_exit_peak_no_ends - calc_phase_nom_to_exit_peak_no_ends;

%Calculate statistics of the phase errors relative to the nominal radiation wave
[phase_nom_error_peak_no_ends_ave, phase_nom_error_peak_no_ends_rms] = phase_calc_phase_nom_error_stats(phase_nom_error_peak_no_ends);
[phase_nom_error_peak_no_ends_enter_ave, phase_nom_error_peak_no_ends_enter_rms] = phase_calc_phase_nom_error_enter_stats(num_pole_stat, phase_nom_error_peak_no_ends);
[phase_nom_error_peak_no_ends_exit_ave, phase_nom_error_peak_no_ends_exit_rms] = phase_calc_phase_nom_error_exit_stats(num_pole_stat, phase_nom_error_to_exit_peak_no_ends);

%Plot the phase error relative to the nominal radiation wave in the undulator
%phase_plot_phase_nom_errors(plt_file, z_peak_no_ends, phase_nom_error_peak_no_ends);
phase_plot_phase_nom_errors_pole(plt_file, z_peak_no_ends, phase_nom_error_peak_no_ends, z_peak);

%Phase Integral -----

%Calculate the phase integral
phase_integral = int_calc_int_y_dx(z, int_bx_dz.^2 + int_by_dz.^2);
phase_integral_end = phase_integral(end);
phase_integral_cell_entrance = spline(z, phase_integral, - cell_length / 2);
phase_integral_cell_exit = spline(z, phase_integral, cell_length / 2);
phase_integral_cell = phase_integral_cell_exit - phase_integral_cell_entrance;

%Plot the phase integral
% phase_plot_phase_integral(plt_file, z, phase_integral);

%Write the phase integral to the data file
phase_dat_phase_integral(dat_file, phase_integral_end, phase_integral_cell);

%Effective Length -----

%Fit the phase integral in the center of the undulator
index_array = find((z > z_peak_no_ends(1)) & (z < z_peak_no_ends(end)));
z_fit = z(index_array);
pi_fit = phase_integral(index_array);
p = polyfit(z_fit, pi_fit, 3);

%Find where the phase integral diverges from the fit at the entrance
index_array = find(z < z_peak_no_ends(1));
z_enter = z(index_array);
pi_enter = phase_integral(index_array);
pi_enter_fit = polyval(p, z_enter);
index_array = find(abs(pi_enter - pi_enter_fit) > .1 * (pi_enter(end) - pi_enter(1)));
z_start = z_enter(index_array(end));

%Find where the phase integral diverges from the fit at the exit
index_array = find(z > z_peak_no_ends(end));
z_exit = z(index_array);
pi_exit = phase_integral(index_array);
pi_exit_fit = polyval(p, z_exit);
index_array = find(abs(pi_exit - pi_exit_fit) > .1 * (pi_exit(end) - pi_exit(1)));
z_end = z_exit(index_array(1));

%Phase length
l_eff = z_end - z_start;

%Write the effective length to the data file
phase_dat_effective_length(dat_file, l_eff);

%Plot the phase integral and effective length
index_array = find((z >= z_start) & (z <= z_end));
z_fit = z(index_array);
pi_fit = polyval(p, z_fit);
phase_plot_phase_integral_effective_length(plt_file, z, phase_integral, z_fit, pi_fit, l_eff);

%Results -----

%Save values in the phas structure
phas.z = z;
phas.z_peak = z_peak;
phas.z_peak_no_ends = z_peak_no_ends;
phas.lambda_u = lambda_u;
phas.slip = slip;
phas.slip_peak = slip_peak;
phas.slip_peak_no_ends = slip_peak_no_ends;
phas.slope = slope;
phas.slip_resid_peak_no_ends = slip_resid_peak_no_ends;
phas.k_eff = k_eff;
phas.kx_eff = kx_eff;
phas.ky_eff = ky_eff;
phas.b_eff = b_eff;
phas.lambda_r = lambda_r;
phas.k_eff_sc = k_eff_sc;
phas.b_eff_sc = b_eff_sc;
phas.lambda_r_sc = lambda_r_sc;
phas.phase = phase;
phas.phase_peak = phase_peak;
phas.phase_peak_no_ends = phase_peak_no_ends;
phas.phase_error_peak_no_ends = phase_error_peak_no_ends;
phas.phase_error_peak_no_ends_ave = phase_error_peak_no_ends_ave;
phas.phase_error_peak_no_ends_rms = phase_error_peak_no_ends_rms;
phas.phase_sc = phase_sc;
phas.phase_peak_sc = phase_peak_sc;
phas.phase_peak_no_ends_sc = phase_peak_no_ends_sc;
phas.phase_error_peak_no_ends_sc = phase_error_peak_no_ends_sc;
phas.phase_error_peak_no_ends_ave_sc = phase_error_peak_no_ends_ave_sc;
phas.phase_error_peak_no_ends_rms_sc = phase_error_peak_no_ends_rms_sc;
phas.phase_cell = phase_cell;
phas.phase_nom = phase_nom;
phas.phase_nom_peak = phase_nom_peak;
phas.phase_nom_peak_no_ends = phase_nom_peak_no_ends;
phas.phase_nom_error_peak_no_ends = phase_nom_error_peak_no_ends;
phas.phase_nom_error_peak_no_ends_ave = phase_nom_error_peak_no_ends_ave;
phas.phase_nom_error_peak_no_ends_rms = phase_nom_error_peak_no_ends_rms;
phas.phase_nom_error_peak_no_ends_enter_ave = phase_nom_error_peak_no_ends_enter_ave;
phas.phase_nom_error_peak_no_ends_enter_rms = phase_nom_error_peak_no_ends_enter_rms;
phas.phase_nom_error_peak_no_ends_exit_ave = phase_nom_error_peak_no_ends_exit_ave;
phas.phase_nom_error_peak_no_ends_exit_rms = phase_nom_error_peak_no_ends_exit_rms;
phas.phase_nom_cell = phase_nom_cell;
phas.phase_integral = phase_integral;
phas.phase_integral_end = phase_integral_end;
phas.phase_integral_cell = phase_integral_cell;
phas.l_eff = l_eff;

%Write a summary of the results to the data file
phase_dat_phase_summary(dat_file, param, phas);

%Done
