%**************************************************************************
%This function performs an analysis of the phase between the light and
%the electron trajectory.  It consists of three parts.  In the first part,
%the slippage, the distance between the light and the electrons, is analyzed.
%In the second part, the undulator parameters are calculated.
%In the third part, the phase is analyzed.
%
%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_norm           %Normalized slippage, slip_norm = 0 at z = 0 (m)
%	slip_norm_peak      %Normalized slippage at the field peaks (m)
%	slip_norm_peak_no_ends	%Normalized slippage at the field peaks excluding the ends (m)
%	slope_norm          %Slope of the slip_norm_peak vs z_peak (no ends) linear fit
%	k_eff               %K value for undulator
%	b_eff               %Effective field strength in undulator (T)
%	lambda_r            %Radiation wavelength (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, error at first pole after end defines zero error (rad)
%   ideal_phase_match   %Ideal phase of pole after end poles (rad)
%   phase_match_error_enter_ave  %Average phase matching error at entrance (rad)
%   phase_match_error_enter_rms  %RMS phase matching error at entrance (rad)
%   phase_match_error_center_ave  %Average phase matching error in center (rad)
%   phase_match_error_center_rms  %RMS phase matching error in center (rad)
%   phase_match_error_exit_ave  %Average phase matching error at in exit (rad)
%   phase_match_error_exit_rms  %RMS phase matching error in exit (rad)
%	phase_cell          %Total phase change between the light and the electrons in the cell (rad)
%   nom_phase_cell      %Ideal phase change in cell (rad)
%   phase_error_cell    %Difference between measured cell phase and ideal cell phase (rad)
%   pme_enter           %Phase matching error to correct at the entrance (rad)
%   phase_shifter_enter %Phase shifter setting to correct the error at the entrance (rad)
%   phase_shim_enter    %Phase shim setting to correct the error at the entrance (rad)
%   pme_exit            %Phase matching error to correct at the exit (rad)
%   phase_shifter_exit  %Phase shifter setting to correct the error at the exit (rad)
%   phase_shim_exit     %Phase shim setting to correct the error at the exit (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, error at first pole after end defines zero error (rad)
%	phase_nom_cell      %Total phase change between the nominal radiation wave and the electrons in the cell (rad)
%
%Zachary Wolf
%12/26/05

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

%Parameter
NUM_POLES_AVE = 10;  %number of poles to average over for phase matching

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

%Slippage -----

%Calculate the slippage as a function of z
slip = slip_calc_slip(gamma, z, x_p, y_p);
%Use the corrected trajectory slopes
%slip = slip_calc_slip(gamma, z, x_p_slope_cor, y_p_slope_cor);

%Find the slippage at the field peaks
slip_peak = slip_calc_slip_at_peaks(z, slip, z_peak);

%Perform a linear fit to the slippage at the field peaks, neglect ends
slip_peak_no_ends = slip_peak(num_end_poles + 1 : length(slip_peak) - num_end_poles);
[slope, intercept, slip_resid_peak_no_ends] = slip_fit_slip_at_peaks(z_peak_no_ends, slip_peak_no_ends);

%Plot the slippage
%slip_plot_slip_polenum(plt_file, param, z, slip);

%Plot the slippage residuals
%slip_plot_slip_resid_polenum(plt_file, param, z_peak_no_ends, slip_resid_peak_no_ends);

%Normalized slippage -----

%Calculate the normalized slippage as a function of z
slip_norm = slip_calc_slip_norm(z, int_bx_dz, int_by_dz);

%Find the normalized slippage at the field peaks
slip_norm_peak = slip_calc_slip_norm_at_peaks(z, slip_norm, z_peak);

%Perform a linear fit to the normalized slippage at the field peaks, neglect ends
slip_norm_peak_no_ends = slip_norm_peak(num_end_poles + 1 : length(slip_norm_peak) - num_end_poles);
[slope_norm, intercept_norm, slip_norm_resid_peak_no_ends] = slip_fit_slip_at_peaks(z_peak_no_ends, slip_norm_peak_no_ends);

%Plot the normalized slippage
%slip_plot_slip_norm(plt_file, param, z, slip_norm);

%Keff, Beff, Lambda_rad -----

%Calculate Keff
k_eff = k_calc_keff(gamma, slope);

%Calculate Keff from the normalized slippage
k_eff_norm = k_calc_keff_norm(slope_norm);

%Check
if abs((k_eff - k_eff_norm) / k_eff) > 1.e-5
    warning('Problem with the Keff calculation.');
end

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

%Calculate the radiation wavelength
lambda_r = k_calc_lambdar(lambda_u, slope);

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

%Calculate Keff at 20 C, it the temperature information is available
if field.temp > 10
    %Calculated Keff at 20C
    k_eff_20c = k_calc_keff_20c(field.temp, k_eff);
    %Write the result to the data file
    k_dat_keff_20c(dat_file, field.temp, k_eff, k_eff_20c);
end

%Write check results of the k_eff calculation to the data file
%k_dat_keff_check(dat_file, gamma, z, slip, z_peak_no_ends, slip_peak_no_ends);

%Phase -----

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

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

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

%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;

%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_ave = mean(phase_error_peak_no_ends);
phase_error_peak_no_ends_rms = std(phase_error_peak_no_ends);

%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_error(plt_file, param, z_peak_no_ends, phase_error_peak_no_ends);

%Write the phase error statistics to the data file
phase_dat_phase_error(dat_file, phase_error_peak_no_ends_ave, phase_error_peak_no_ends_rms);

%Phase matching -----

%The number of poles must be correct in order for the phase matching to work
if num_poles ~= length(z_peak)
    errordlg(sprintf(['pid_anal_phase: problem, wrong number of peaks\nlength(z_peak) = ', num2str(length(z_peak)), '\nnum_poles = ', num2str(num_poles)]));
    field_dat_num_poles_err(dat_file, length(z_peak), num_poles);
    num_poles = length(z_peak);  %Proceed after documenting error
end

%Calculate the ideal phase through the undulator
%ideal_phase_match = ideal phase at the pole after the end poles
%Ninterior = NUM_POLES - (2 * NUM_END_POLES)
%Nbetween_interior = Ninterior - 1
%ideal_phase_match + (NUM_POLES - (2 * NUM_END_POLES) - 1) * pi + ideal_phase_match = NOM_PHASE_CELL
%2 * ideal_phase_match = NOM_PHASE_CELL - (NUM_POLES - (2 * NUM_END_POLES) - 1) * pi
%ideal_phase_match = (NOM_PHASE_CELL - (NUM_POLES - (2 * NUM_END_POLES) - 1) * pi) / 2
ideal_phase_match = (nom_phase_cell - (num_poles - (2 * num_end_poles) - 1) * pi) / 2;
ideal_phase_peak_no_ends = zeros(1, num_poles - 2 * num_end_poles);
ideal_phase_peak_no_ends(1) = ideal_phase_match;
for ii = 2 : length(ideal_phase_peak_no_ends)
    ideal_phase_peak_no_ends(ii) = ideal_phase_peak_no_ends(ii - 1) + pi;
end

%Phase matching at the entrance
phase_advance_enter = phase_peak_no_ends - 0;  %Phase advance from cell boundary to a pole, phase of cell boundary is zero
ideal_phase_advance_enter = ideal_phase_peak_no_ends - 0;  %Ideal phase advance from cell boundary to a pole, phase of cell boundary is zero
phase_match_error_enter = phase_advance_enter - ideal_phase_advance_enter;
phase_match_error_enter_ave = mean(phase_match_error_enter(1 : NUM_POLES_AVE));
phase_match_error_enter_rms = std(phase_match_error_enter(1 : NUM_POLES_AVE));

%Phase matching in the center
% phase_advance_center = phase_peak_no_ends - phase_peak_no_ends(1);  %Phase advance from first pole after ends to a center pole
% ideal_phase_advance_center = ideal_phase_peak_no_ends - ideal_phase_peak_no_ends(1);  %Ideal phase advance from first pole after ends to a center pole
phase_advance_center = phase_peak_no_ends;  %Phase advance from cell boundary
ideal_phase_advance_center = ideal_phase_peak_no_ends;  %Ideal phase advance from cell boundary
phase_match_error_center = phase_advance_center - ideal_phase_advance_center;
phase_match_error_center_ave = mean(phase_match_error_center);
phase_match_error_center_rms = std(phase_match_error_center);

%Phase matching at the exit
phase_advance_exit = phase_cell_exit - phase_peak_no_ends;  %Phase advance from a pole to the cell exit
ideal_phase_advance_exit = nom_phase_cell - ideal_phase_peak_no_ends;  %Ideal phase advance from a pole to the cell exit
phase_match_error_exit = phase_advance_exit - ideal_phase_advance_exit;
phase_match_error_exit_ave = mean(phase_match_error_exit(end - NUM_POLES_AVE + 1 : end));
phase_match_error_exit_rms = std(phase_match_error_exit(end - NUM_POLES_AVE + 1 : end));

%Phase error in cell
phase_error_cell = phase_cell - nom_phase_cell;

%Calculate entrance corrections, use error through the center to determine the entrance correction
pme_enter = phase_match_error_center_ave;
for ii = -100 : 100
    if (ii * 2 * pi > pme_enter) && ((ii - 1) * 2 * pi <= pme_enter)
        phase_shifter_enter = ii * 2 * pi - pme_enter;
        phase_shim_enter = pme_enter - (ii - 1) * 2 * pi;
        break;
    end
end

%Calculate exit corrections, use cell error minus the entrance error
pme_exit = phase_error_cell - phase_match_error_center_ave;
for ii = -100 : 100
    if (ii * 2 * pi > pme_exit) && ((ii - 1) * 2 * pi <= pme_exit)
        phase_shifter_exit = ii * 2 * pi - pme_exit;
        phase_shim_exit = pme_exit - (ii - 1) * 2 * pi;
        break;
    end
end

%Plot the phase matching errors
phase_match_error_peak_no_ends = phase_peak_no_ends - ideal_phase_peak_no_ends;
phase_plot_phase_match_error(plt_file, z_peak_no_ends, phase_match_error_peak_no_ends);

%Write the phase matching method to the data file
phase_dat_phase_match_method(dat_file, nom_phase_cell, num_poles, num_end_poles, ideal_phase_match);

%Write the phase matching results to the data file
phase_dat_phase_match_results(dat_file, NUM_POLES_AVE, phase_match_error_enter_ave, phase_match_error_enter_rms, phase_match_error_center_ave, phase_match_error_center_rms, phase_match_error_exit_ave, phase_match_error_exit_rms);

%Write the cell phase results
phase_dat_phase_cell_results(dat_file, phase_cell, nom_phase_cell, phase_error_cell);

%Write the phase matching correction results
phase_dat_phase_match_corr(dat_file, pme_enter, phase_shifter_enter, phase_shim_enter, pme_exit, phase_shifter_exit, phase_shim_exit);

%Phase To Nominal Wave -----

%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;

%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
%pole_num_no_ends = 1 : length(phase_nom_peak_no_ends);
%phase_nom_error_peak_no_ends = (phase_nom_peak_no_ends - phase_nom_peak_no_ends(1)) - (pi * (pole_num_no_ends - 1));

%Calculate the phase error relative to the nominal radiation wave in the undulator
phase_nom_error_peak_no_ends = phase_nom_peak_no_ends - ideal_phase_peak_no_ends;
phase_nom_error_peak_no_ends_ave = mean(phase_nom_error_peak_no_ends);
phase_nom_error_peak_no_ends_rms = std(phase_nom_error_peak_no_ends);

%Plot the phase error relative to the nominal radiation wave in the undulator
%phase_plot_phase_nom_error(plt_file, param, z_peak_no_ends, phase_nom_error_peak_no_ends);

%Write the nominal wave phase error statistics to the data file
%phase_dat_phase_nom_error(dat_file, phase_nom_error_peak_no_ends_ave, phase_nom_error_peak_no_ends_rms);

%Write the cell nominal wave phase results
%phase_dat_phase_nom_cell_results(dat_file, phase_nom_cell, nom_phase_cell);

%Phase matching to nominal wave -----

%Phase matching at the entrance
phase_nom_advance_enter = phase_nom_peak_no_ends - 0;  %Phase advance from cell boundary to a pole, phase of cell boundary is zero
ideal_phase_advance_enter = ideal_phase_peak_no_ends - 0;  %Ideal phase advance from cell boundary to a pole, phase of cell boundary is zero
phase_nom_match_error_enter = phase_nom_advance_enter - ideal_phase_advance_enter;
phase_nom_match_error_enter_ave = mean(phase_nom_match_error_enter(1 : NUM_POLES_AVE));
phase_nom_match_error_enter_rms = std(phase_nom_match_error_enter(1 : NUM_POLES_AVE));

%Phase matching in the center
% phase_nom_advance_center = phase_nom_peak_no_ends - phase_nom_peak_no_ends(1);  %Phase advance from first pole after ends to a center pole
% ideal_phase_advance_center = ideal_phase_peak_no_ends - ideal_phase_peak_no_ends(1);  %Ideal phase advance from first pole after ends to a center pole
phase_nom_advance_center = phase_nom_peak_no_ends;  %Phase advance from cell boundary
ideal_phase_advance_center = ideal_phase_peak_no_ends;  %Ideal phase advance from cell boundary
phase_nom_match_error_center = phase_nom_advance_center - ideal_phase_advance_center;
phase_nom_match_error_center_ave = mean(phase_nom_match_error_center);
phase_nom_match_error_center_rms = std(phase_nom_match_error_center);

%Phase matching at the exit
phase_nom_advance_exit = phase_nom_cell_exit - phase_nom_peak_no_ends;  %Phase advance from a pole to the cell exit
ideal_phase_advance_exit = nom_phase_cell - ideal_phase_peak_no_ends;  %Ideal phase advance from a pole to the cell exit
phase_nom_match_error_exit = phase_nom_advance_exit - ideal_phase_advance_exit;
phase_nom_match_error_exit_ave = mean(phase_nom_match_error_exit(end - NUM_POLES_AVE + 1 : end));
phase_nom_match_error_exit_rms = std(phase_nom_match_error_exit(end - NUM_POLES_AVE + 1 : end));

%Write the phase matching method to the data file
%phase_dat_phase_match_method(dat_file, nom_phase_cell, num_poles, num_end_poles, ideal_phase_match);

%Write the phase matching results to the data file
%phase_dat_phase_nom_match_results(dat_file, NUM_POLES_AVE, phase_nom_match_error_enter_ave, phase_nom_match_error_enter_rms, phase_nom_match_error_center_ave, phase_nom_match_error_center_rms, phase_nom_match_error_exit_ave, phase_nom_match_error_exit_rms);

%Temperature effects -----

%Calculate the phase in the cell at 20 C, it the temperature information is available
if temp > 10
    %Calculated phase_nom_cell at 20C
    phase_nom_cell_20c = phase_calc_phase_nom_cell_20c(temp, k_eff, phase_nom_cell);
    %Write the result to the data file
    phase_dat_phase_nom_cell_20c(dat_file, temp, phase_nom_cell, phase_nom_cell_20c);
end

%Save results -----

%Save the quantities 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_norm = slip_norm;
phas.slip_norm_peak = slip_norm_peak;
phas.slip_norm_peak_no_ends = slip_norm_peak_no_ends;
phas.slope_norm = slope_norm;
phas.k_eff = k_eff;
phas.b_eff = b_eff;
phas.lambda_r = lambda_r;
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.ideal_phase_match = ideal_phase_match;
phas.phase_match_error_enter_ave = phase_match_error_enter_ave;
phas.phase_match_error_enter_rms = phase_match_error_enter_rms;
phas.phase_match_error_center_ave = phase_match_error_center_ave;
phas.phase_match_error_center_rms = phase_match_error_center_rms;
phas.phase_match_error_exit_ave = phase_match_error_exit_ave;
phas.phase_match_error_exit_rms = phase_match_error_exit_rms;
phas.phase_cell = phase_cell;
phas.nom_phase_cell = nom_phase_cell;
phas.phase_error_cell = phase_error_cell;
phas.pme_enter = pme_enter;
phas.phase_shifter_enter = phase_shifter_enter;
phas.phase_shim_enter = phase_shim_enter;
phas.pme_exit = pme_exit;
phas.phase_shifter_exit = phase_shifter_exit;
phas.phase_shim_exit = phase_shim_exit;
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_cell = phase_nom_cell;

%Results -----

%Write a summary of the results to the data file
% param.nom_phase_match = ideal_phase_match;
% phase_dat_phase_results(dat_file, param, phas);

%Write the results to the data file
%phase_dat_phase(dat_file, param, phas);

%Done
