%**************************************************************************
%Planar Insertion Device Analysis
%
%This program analyzes the data from a planar insertion device.
%It uses the convention that the primary field is in the Y direction,
%the beam is moving in the Z direction, and the main trajectory is in
%the X direction.  Error fields in By give errors in the the X trajectory.
%Error fields in Bx give errors in the Y trajectory.
%
%First the peak and zero fields and their positions and statistics are determined.
%The trajectory of the beam is then calculated.  The phase shift between
%the light and the electron beam (for a given energy) is determined.
%Results are written to a file.
%
%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
%   corr_required       %Logical variable, apply a correction to the measured fields, true or false
%   corr_bx_mult        %Multiplicative correction factor for Bx
%   corr_by_mult        %Multiplicative correction factor for By
%   corr_bx_add         %Additive constant correction for Bx (T)
%   corr_by_add         %Additive constant correction for By (T)
%   corr_bxy_add_file   %Correction file name including directory, additive for Bx and By, 'none' for no correction
%   corr_bxy_add_file_type  %Correction file type, gives format
%   set_i1xy_required   %Logical variable whether to apply a correction to set the first field integrals to a specified value
%   set_i1x             %Value of I1x after additive correction (Tm)
%   set_i1y             %Value of I1y after additive correction (Tm)
%   gamma               %Lorentz factor
%   nom_lambda_u        %Nominal undulator period (m)
%   nom_k_value         %Nominal undulator K parameter value
%   nom_lambda_r        %Nominal radiation wavelength (m)
%   num_poles           %Number of poles in the undulator
%   num_end_poles       %Number of poles used to determine end effects
%   cell_length         %Cell length, undulator length plus drift length at each end
%   nom_phase_cell      %Nominal phase change through the cell (rad)
%   save_as_final_result%Save output as a final result, true or false
%   result_file_dir     %Directory where final results are saved
%   result_file_name    %Name of file for final result
%   fig_close_on_exit   %Close figures when program exits, true or false
%   fig_vis             %Figures visible on screen, true or false
%
%field, structure containing magnetic field information
%   zm                  %Z positions at which the field measurements are made (m)
%   z                   %Translated Z, with z=0 at undulator center (m)
%   bx                  %Measured horizontal field (T)
%   by                  %Measured vertical field (T)
%   z_peak_bx           %Z positions of the Bx peak fields (m)
%   z_peak_by           %Z positions of the By peak fields (m)
%   z_peak_bx_no_ends   %Z positions of the Bx peak fields, end data removed (m)
%   z_peak_by_no_ends   %Z positions of the By peak fields, end data removed (m)
%   z_peak              %Z positions of the peak fields of the larger of Bx or By (m)
%   z_peak_no_ends      %Z positions of the peak fields of the larger of Bx or By, end data removed (m)
%   z_zero_bx           %Z positions of the Bx field zeros (m)
%   z_zero_by           %Z positions of the By field zeros (m)
%   z_zero_bx_no_ends   %Z positions of the Bx zero fields, end data removed (m)
%   z_zero_by_no_ends   %Z positions of the By zero fields, end data removed (m)
%   z_zero              %Z positions of the zero points of the larger of Bx or By (m)
%   z_zero_no_ends      %Z positions of the zero points of the larger of Bx or By, end data removed (m)
%   bx_peak             %Fitted Bx peak values (T)
%   by_peak             %Fitted By peak values (T)
%   bx_peak_no_ends     %Fitted Bx peak values, end data removed (T)
%   by_peak_no_ends     %Fitted By peak values, end data removed (T)
%   lambda_u            %Measured average undulator period (m)
%   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)
%   temp                %Average undulator temperature (C)
%
%traj, structure containing trajectory information
%   z                   %Translated Z, with z=0 at undulator center (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)
%   x_slope_cor         %X trajectory minus linear fit
%   x_p_slope_cor       %X trajectory slope minus linear fit slope
%   y_slope_cor         %Y trajectory minus linear fit
%   y_p_slope_cor       %Y trajectory slope minus linear fit slope
%
%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)
%
%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
%
%Pure sinusoidal field:
%by = B0 cos(2 pi z / lambda_u)
%x' = - (q / gamma m vz) int by dz
%x' = - (q / gamma m vz) (B0 lambda_u/ 2 pi) sin(2 pi z / lambda_u)
%x' = (K / gamma) sin(2 pi z / lambda_u)    (define K)
%K := - q B0 lambda_u / m vz 2 pi        (can approximate vz by v, or c)
%
%Slippage in a planar sinusoidal field, y' = 0:
%slip = (1/2) int (1/gamma^2 + x'^2) d_z
%slip = (1/2) int (1/gamma^2 + (K^2/gamma^2) sin^2(k z)) d_z
%sin^2(k z) = (1/2) (1 - cos(2 k z))
%slip = (1/2) (1/gamma^2) z + (1/2) (K^2/gamma^2) (1/2) z - (1/2) (K^2/gamma^2) (1/2) int cos(2 k z) dz
%slip = (1/2) (1/gamma^2) (1 + (1/2) K^2) z - (1/4) (K^2/gamma^2) (1/2k) sin(2 k z)
%
%For a sinusoidal field, at the field peaks:
%by = B0 cos(2 pi z / lambda_u) = B0 cos(k z)
%at the peaks, k z = 0, pi, 2 pi, ...
%slip_peak = (1/2) (1/gamma^2) (1 + (1/2) K^2) z
%
%Fit slippage at peaks:
%define slope of slippage at peak fields vs z := slope
%
%Assume the measured slope came from a sinusoidal field:
%calculate Keff, the K value giving the measured slope
%(1/2) (1/gamma^2) (1 + (1/2) Keff^2) = slope
%(1 + (1/2) Keff^2) = 2 gamma^2 slope
%Keff^2 = 2 (2 gamma^2 slope - 1)
%Keff = sqrt(2 (2 gamma^2 slope - 1))
%
%Determine Beff from Keff:
%K := - q B0 lambda_u / m vz 2 pi
%Keff := - q Beff lambda_u / m vz 2 pi
%Beff = Keff (m vz 2 pi / (- q) lambda_u)
%
%Radiation wavelength:
%radiation has slippage = lambda_r when delta_z = lambda_u
%delta_slippage / delta_z = lambda_r / lambda_u
%slope = lambda_r / lambda_u
%lambda_r = lambda_u slope
%lambda_r = (1/2) (1/gamma^2) (1 + (1/2) Keff^2) lambda_u
%
%Phase:
%phase = (2 pi / lambda_r) slip
%in the undulator,
%slip = slip_0 + slope (z - z_0) + slip_resid
%slope = lambda_r / lambda_u
%phase = (2 pi / slope lambda_u) slip
%phase = (2 pi / lambda_r) slip_0 + (2 pi / lambda_u) (z - z_0) + (2 pi / lambda_r) slip_resid
%phase = phase_0 + (2 pi / lambda_u) (z - z_0) + (2 pi / lambda_r) slip_resid
%
%Phase error:
%phase_error = phase - [phase_0 + (2 pi / lambda_u) (z - z_0)]
%phase_error = (2 pi / lambda_r) slip_resid
%
%Phase relative to nominal radiation wave:
%phase_nom = (2 pi / nom_lambda_r) slip
%
%Phase error relative to nominal radiation wave (in undulator):
%phase_nom_error = (phase_nom - phase_nom(start)) - pi (pole_num - pole_num(start))
%
%Normalized slippage, alternate technique, remove gamma dependence from calculations
%slip = (1/2) int (1/gamma^2 + x'^2 + y'^2) d_z
%slip = (1/2) int {1/gamma^2 + [- (q / gamma m vz) int by dz]^2 + [(q / gamma m vz) int bx dz]^2} d_z
%slip = (1/2) (1/gamma^2) int {1 + [- (q / m vz) int by dz]^2 + [(q / m vz) int bx dz]^2} d_z
%slip = (1/2) (1/gamma^2) int {1 + (q / m vz)^2 [(int by dz)^2 + (int bx dz)^2]} d_z
%slip_norm := gamma^2 slip
%slip_norm = (1/2) int {1 + (q / m vz)^2 [(int by dz)^2 + (int bx dz)^2]} d_z
%normalized slippage in a planar sinusoidal field, y' = 0:
%slip = (1/2) (1/gamma^2) (1 + (1/2) K^2) z - (1/4) (K^2/gamma^2) (1/2k) sin(2 k z)
%slip_norm = (1/2) (1 + (1/2) K^2) z - (1/4) K^2 (1/2k) sin(2 k z)
%for a sinusoidal field, at the field peaks:
%slip_peak = (1/2) (1/gamma^2) (1 + (1/2) K^2) z
%slip_norm_peak = (1/2) (1 + (1/2) K^2) z
%define slope of normalized slippage at peak fields vs z := slope_norm
%assume the measured slope came from a sinusoidal field:
%calculate Keff, the K value giving the measured slope
%(1/2) (1 + (1/2) Keff^2) = slope_norm
%(1 + (1/2) Keff^2) = 2 slope_norm
%Keff^2 = 2 (2 slope_norm - 1)
%Keff = sqrt(2 (2 slope_norm - 1))
%
%Zachary Wolf, 7/25/03, 12/23/05, 08/04/07

function [] = pid_anal()

%Set up program parameters using the user interface
%Create the updated pid_anal_param.m file which is run to set the parameter values
%pid_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_anal_param();

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

%Make file names
dat_file = strcat(lower(param.input_file_dir), 'pid_anal.dat');
plt_file = strcat(lower(param.input_file_dir), 'pid_anal.ps');
mat_file = strcat(param.input_file_dir, 'pid_anal.mat');

%Initialize the data file
util_dat_init(dat_file);

%Write a message to the data file
msg{1} = ' ';
msg{2} = 'Planar Insertion Device Analysis';
msg{3} = '--------------------------------';
util_dat_msg(dat_file, msg);

%Write a header to the data file
util_dat_header(dat_file, header);

%Write the parameter values to the data file
pid_anal_dat_param(dat_file, param);

%Initialize the plot file
util_plt_init(plt_file, msg, param, header);

%Apply a field correction, if required
if param.corr_required == true
    %Apply the field correction
	[bxc, byc] = field_corr_apply_field_corr(param, zm, bx, by);

    %Write a message to the data file that the field was corrected
    field_corr_dat_field_corr(dat_file, param);

    %Plot the correction
    field_corr_plot_field_corr(plt_file, zm, bxc - bx, byc - by);

    %Use the corrected fields
    bx = bxc;
    by = byc;
end

%!!!!!!!!!!!!
%Special cases-----

%Print a warning if values are changed
% warndlg('Special case values are changed');

%Select part of a scan
%  zm_min = 2.00;
%  zm_max = 2.64;
%  index_array = find((zm >= zm_min) & (zm <= zm_max));
%  zm = zm(index_array);
%  bx = bx(index_array);
%  by = by(index_array);

%Reverse the measurement direction
zm = - fliplr(zm);
bx = - fliplr(bx);
by = fliplr(by);
field_corr_dat_meas_dir_rev(dat_file);

% %Add a shim in software
% zm_shim = zm(1) + (zm(end) - zm(1)) * 0.4;
% i1x_shim = -12.e-6; %(Tm)
% i1y_shim = 0;  %(Tm)
% dz_shim = 0.01;  %(m)
% bx_shim = i1x_shim / dz_shim;
% by_shim = i1y_shim / dz_shim;
% index_array = find((zm > zm_shim - dz_shim / 2) & (zm < zm_shim + dz_shim / 2));
% zm_apply_shim = zeros(1, length(zm));
% zm_apply_shim(index_array) = 1;
% bx = bx + bx_shim * zm_apply_shim;
% by = by + by_shim * zm_apply_shim;
% field_corr_dat_apply_soft_shim(dat_file, zm_shim, i1x_shim, i1y_shim);
% 
% %Add a shim in software
% zm_shim = zm(1) + (zm(end) - zm(1)) * 0.9;
% i1x_shim = 12.e-6; %(Tm)
% i1y_shim = 0;  %(Tm)
% dz_shim = 0.01;  %(m)
% bx_shim = i1x_shim / dz_shim;
% by_shim = i1y_shim / dz_shim;
% index_array = find((zm > zm_shim - dz_shim / 2) & (zm < zm_shim + dz_shim / 2));
% zm_apply_shim = zeros(1, length(zm));
% zm_apply_shim(index_array) = 1;
% bx = bx + bx_shim * zm_apply_shim;
% by = by + by_shim * zm_apply_shim;
% field_corr_dat_apply_soft_shim(dat_file, zm_shim, i1x_shim, i1y_shim);

%!!!!!!!!!!!!

%Set the field integrals, if required
if param.set_i1xy_required == true
    %Apply the field correction
	[bx_add, by_add, bxc, byc] = field_corr_apply_field_int_corr(param, zm, bx, by);

    %Write a message to the data file that the field was corrected
    field_corr_dat_field_int_corr(dat_file, param, bx_add, by_add);

    %Use the corrected fields
    bx = bxc;
    by = byc;
end

%Get the undulator gap and temperature
[gap] = pid_anal_get_und_gap(header);
[temp] = pid_anal_get_und_temp(header);
pid_anal_dat_und_gap(dat_file, gap);
pid_anal_dat_und_temp(dat_file, temp);

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

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

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

%Output the results to the input file directory
save(mat_file, 'param', 'header', 'gap', 'temp', 'field', 'traj', 'phas');

%Save the results to the results folder, if desired
if param.save_as_final_result == true
%     %Dat file
%     result_dat_file = strcat(param.result_file_dir, param.result_file_name);
% 	pid_anal_dat_results(result_dat_file, param, header, field, traj, phas);
%     %Mat file
%     index = strfind(param.result_file_name, '.');
%     result_mat_file = strcat(param.result_file_dir, param.result_file_name(1 : index), 'mat');
%     save(result_mat_file, 'param', 'header', 'field', 'traj', 'phas');
end

%Write the results to the workspace to share with other programs
wrk_file = 'c:\analysis\workspace\pid_anal.wrk';
pid_anal_dat_vals(wrk_file, param, field, traj, phas);

%Remove all figures, if desired
if param.fig_vis == false
    close('all','hidden');
else
    if param.fig_close_on_exit == true
        close('all','hidden');
    end
end

%Done
