%**************************************************************************
%This function performs an alalysis of the magnetic field.
%
%Input:
%param, parameter structure
%dat_file, output data file name
%plt_file, output plot file name
%zm, z positions at which the field measurements are made (m)
%bx, measured horizontal field (T)
%by, measured vertical field (T)
%temp, average undulator temperature (C)
%
%Output:
%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)
%
%Zachary Wolf
%12/23/05

function [field] = pid_anal_field(param, dat_file, plt_file, zm, bx, by, temp)

%Field -----

% %Find and fit the measured field peaks
% [zm_peak, by_peak] = field_find_field_peaks(zm, by);
% 
% %Find and fit the measured field zeros
% [zm_zero] = field_find_field_zeros(zm, by);

%Find and fit the measured Bx field peaks
[zm_peak_bx, bx_peak] = field_find_peaks(zm, bx);

%Find and fit the measured Bx field zeros
[zm_zero_bx] = field_find_zeros_btwn_peaks(zm, bx, zm_peak_bx);

%Find and fit the measured By field peaks
[zm_peak_by, by_peak] = field_find_peaks(zm, by);

%Find and fit the measured By field zeros
[zm_zero_by] = field_find_zeros_btwn_peaks(zm, by, zm_peak_by);
  
% %Checks
% %Required for phase matching
% if length(zm_peak) ~= param.num_poles
%     errordlg(sprintf(['pid_anal_field: found wrong number of peaks\nlength(zm_peak) = ', num2str(length(zm_peak)), '\nnum_poles = ', num2str(param.num_poles)]));
% end
% if length(zm_zero) ~= param.num_poles - 1
%     errordlg(sprintf(['pid_anal_field: found wrong number of zeros\nlength(zm_zero) = ', num2str(length(zm_zero)), '\nnum_poles - 1 = ', num2str(param.num_poles - 1)]));
% end
  
%Checks
if (length(zm_peak_bx) > 1) && (abs(length(zm_peak_bx) - param.num_poles) > 0)
    warndlg(sprintf(['pid_anal_field: found wrong number of Bx peaks\nlength(zm_peak_bx) = ', num2str(length(zm_peak_bx)), '\nnum_poles = ', num2str(param.num_poles)]));
end
if (length(zm_zero_bx) > 1) && (abs(length(zm_zero_bx) - (param.num_poles - 1)) > 0)
    warndlg(sprintf(['pid_anal_field: found wrong number of Bx zeros\nlength(zm_zero_bx) = ', num2str(length(zm_zero_bx)), '\nnum_poles - 1 = ', num2str(param.num_poles - 1)]));
end
if (length(zm_peak_by) > 1) && (abs(length(zm_peak_by) - param.num_poles) > 0)
    warndlg(sprintf(['pid_anal_field: found wrong number of By peaks\nlength(zm_peak_by) = ', num2str(length(zm_peak_by)), '\nnum_poles = ', num2str(param.num_poles)]));
end
if (length(zm_zero_by) > 1) && (abs(length(zm_zero_by) - (param.num_poles - 1)) > 0)
    warndlg(sprintf(['pid_anal_field: found wrong number of By zeros\nlength(zm_zero_by) = ', num2str(length(zm_zero_by)), '\nnum_poles - 1 = ', num2str(param.num_poles - 1)]));
end

% %Translate Z coordinates so the ID center is at z=0
% zm_peak_mean = mean(zm_peak);
% zm_zero_mean = mean(zm_zero);
% z = zm - zm_zero_mean;
% z_peak = zm_peak - zm_zero_mean;
% z_zero = zm_zero - zm_zero_mean;

%Translate Z coordinates so the center is at z=0
bx_max = max(abs(bx_peak));
by_max = max(abs(by_peak));
if (bx_max > by_max)
    zm_ctr = mean(zm_zero_bx);
else
    zm_ctr = mean(zm_zero_by);
end
z = zm - zm_ctr;
z_peak_bx = zm_peak_bx - zm_ctr;
z_zero_bx = zm_zero_bx - zm_ctr;
z_peak_by = zm_peak_by - zm_ctr;
z_zero_by = zm_zero_by - zm_ctr;
  
% %Check
% if abs(zm_peak_mean - zm_zero_mean) > 2.e-5
%     msgbox('pid_anal_field: problem finding undulator center', 'Warning');
% end

% %Neglect end data for some calculations
% z_peak_no_ends = z_peak(param.num_end_poles + 1 : length(z_peak) - param.num_end_poles);
% by_peak_no_ends = by_peak(param.num_end_poles + 1 : length(by_peak) - param.num_end_poles);
% z_zero_no_ends = z_zero(param.num_end_poles + 1 : length(z_peak) - param.num_end_poles);

%Neglect end data for some calculations
bx_peak_no_ends = 0;
z_peak_bx_no_ends = 0;
z_zero_bx_no_ends = 0;
if length(z_peak_bx) > 1
    bx_peak_no_ends = bx_peak(param.num_end_poles + 1 : length(bx_peak) - param.num_end_poles);
    z_peak_bx_no_ends = z_peak_bx(param.num_end_poles + 1 : length(z_peak_bx) - param.num_end_poles);
    z_zero_bx_no_ends = z_zero_bx(param.num_end_poles + 1 : length(z_zero_bx) - param.num_end_poles);
end
by_peak_no_ends = 0;
z_peak_by_no_ends = 0;
z_zero_by_no_ends = 0;
if length(z_peak_by) > 1
    by_peak_no_ends = by_peak(param.num_end_poles + 1 : length(by_peak) - param.num_end_poles);
    z_peak_by_no_ends = z_peak_by(param.num_end_poles + 1 : length(z_peak_by) - param.num_end_poles);
    z_zero_by_no_ends = z_zero_by(param.num_end_poles + 1 : length(z_zero_by) - param.num_end_poles);
end

%Define peak field positions of the larger of Bx or By
bx_max = max(abs(bx_peak_no_ends));
by_max = max(abs(by_peak_no_ends));
if (bx_max > by_max)
    z_peak = z_peak_bx;
    z_peak_no_ends = z_peak_bx_no_ends;
    z_zero = z_zero_bx;
    z_zero_no_ends = z_zero_bx_no_ends;
else
    z_peak = z_peak_by;
    z_peak_no_ends = z_peak_by_no_ends;
    z_zero = z_zero_by;
    z_zero_no_ends = z_zero_by_no_ends;
end

% %Calculate the undulator period
% [lambda_u] = field_calc_lambda_u(z_zero_no_ends);

%Calculate the undulator period
lambda_u_bx = 0;
if length(z_peak_bx) > 1
    [lambda_u_bx] = field_calc_lambda_u(z_zero_bx_no_ends);
end
lambda_u_by = 0;
if length(z_peak_by) > 1
    [lambda_u_by] = field_calc_lambda_u(z_zero_by_no_ends);
end
if lambda_u_bx == 0
    lambda_u = lambda_u_by;
elseif lambda_u_by == 0
    lambda_u = lambda_u_bx;
elseif bx_max > 10 * by_max
    lambda_u = lambda_u_bx;
elseif by_max > 10 * bx_max
    lambda_u = lambda_u_by;
else
    lambda_u = (lambda_u_bx + lambda_u_by) / 2;
end

%Plot the field data
field_plot_bxby_vs_zm(plt_file, zm, bx, by);
% field_plot_b(plt_file, param, 'Bx', z, bx);
% field_plot_b(plt_file, param, 'By', z, by);
% field_plot_bxy(plt_file, param, z, bx, by);
% field_plot_abs_bxy(plt_file, param, z, bx, by);

%Plot By with the peaks and zeros
% field_plot_by_peaks_zeros(plt_file, param, z, by, z_peak, by_peak, z_zero);
%field_plot_bx_vs_z_peaks_zeros(plt_file, z, bx, bx_peak, z_peak_bx, z_zero_bx);
%field_plot_by_vs_z_peaks_zeros(plt_file, z, by, by_peak, z_peak_by, z_zero_by);
field_plot_bxby_vs_z_peaks_zeros(plt_file, z, bx, bx_peak, z_peak_bx, z_zero_bx, by, by_peak, z_peak_by, z_zero_by);

%Plot the By peaks vs pole number
%field_plot_by_peaks(plt_file, param, by_peak);

%Plot the magnitudes of the peak fields excluding the ends
% field_plot_abs_fit_peaks_no_ends(plt_file, param, z_peak_no_ends, by_peak_no_ends);
field_plot_bxby_abs_peaks_no_ends(plt_file, z_peak_bx, z_peak_bx_no_ends, bx_peak_no_ends, z_peak_by, z_peak_by_no_ends, by_peak_no_ends);

%Plot the peak and zero field spacings, exclude the ends
% field_plot_peak_spacing(plt_file, param, z_peak_no_ends);
% field_plot_zero_spacing(plt_file, param, z_zero_no_ends);
field_plot_bxby_peak_zero_spacing_no_ends(plt_file, z_peak_bx, z_peak_bx_no_ends, z_zero_bx, z_zero_bx_no_ends, z_peak_by, z_peak_by_no_ends, z_zero_by, z_zero_by_no_ends);

%Write a summary of the peak field results to the data file
% field_dat_peaks_zeros_summary(param, dat_file, z_peak, by_peak, z_zero, lambda_u);
field_dat_bx_peaks_zeros_summary(dat_file, bx_peak, bx_peak_no_ends, z_peak_bx, z_peak_bx_no_ends, z_zero_bx, z_zero_bx_no_ends, lambda_u_bx);
field_dat_by_peaks_zeros_summary(dat_file, by_peak, by_peak_no_ends, z_peak_by, z_peak_by_no_ends, z_zero_by, z_zero_by_no_ends, lambda_u_by);

%Field integrals -----

%Calculate the first integral of Bx and By
int_bx_dz = int_calc_int(z, bx);
int_by_dz = int_calc_int(z, by);

% %Plot the first integrals of Bx and By
% field_plot_i1x(plt_file, z, int_bx_dz);
% field_plot_i1y(plt_file, z, int_by_dz);

%Calculate the second integral of Bx and By
int_int_bx_dz_dz = int_calc_int(z, int_bx_dz);
int_int_by_dz_dz = int_calc_int(z, int_by_dz);

% %Plot the second integrals of Bx and By
% field_plot_i2x(plt_file, z, int_int_bx_dz_dz);
% field_plot_i2y(plt_file, z, int_int_by_dz_dz);

%Filters -----

%Filter width
lambda = lambda_u;
%lambda = param.nom_lambda_u;

%Apply a lambda filter to Bx and By
bx_filt = filt_apply_lambda_filt(lambda, z, bx);
by_filt = filt_apply_lambda_filt(lambda, z, by);

%Apply a lambda filter to the first integral of Bx and By
int_bx_dz_filt = filt_apply_lambda_filt(lambda, z, int_bx_dz);
int_by_dz_filt = filt_apply_lambda_filt(lambda, z, int_by_dz);

%Apply a lambda filter to the second integral of Bx and By
int_int_bx_dz_dz_filt = filt_apply_lambda_filt(lambda, z, int_int_bx_dz_dz);
int_int_by_dz_dz_filt = filt_apply_lambda_filt(lambda, z, int_int_by_dz_dz);

%Plot the filtered first integrals of Bx and By
%field_plot_i1x_filt(plt_file, z, int_bx_dz, int_bx_dz_filt);
%field_plot_i1y_filt(plt_file, z, int_by_dz, int_by_dz_filt);

%Plot the filtered second integrals of Bx and By
%field_plot_i2x_filt(plt_file, z, int_int_bx_dz_dz, int_int_bx_dz_dz_filt);
%field_plot_i2y_filt(plt_file, z, int_int_by_dz_dz, int_int_by_dz_dz_filt);
field_plot_i2xi2y_filt(plt_file, z, int_int_bx_dz_dz, int_int_bx_dz_dz_filt, int_int_by_dz_dz, int_int_by_dz_dz_filt);

%Plot the field integrals
%int_plot_b_int(plt_file, 'Bx', z, bx, bx_filt, int_bx_dz, int_bx_dz_filt, int_int_bx_dz_dz, int_int_bx_dz_dz_filt);
%int_plot_b_int(plt_file, 'By', z, by, by_filt, int_by_dz, int_by_dz_filt, int_int_by_dz_dz, int_int_by_dz_dz_filt);

%Write the overall field integrals to the data file
int_dat_b_int(dat_file, int_bx_dz(end), int_int_bx_dz_dz(end), int_by_dz(end), int_int_by_dz_dz(end));

%Save Results -----

%Save quantities in the field structure
field.zm = zm;
field.z = z;
field.bx = bx;
field.by = by;
% field.z_peak = z_peak;
% field.by_peak = by_peak;
% field.z_zero = z_zero;
field.z_peak_bx = z_peak_bx;
field.z_peak_by = z_peak_by;
field.z_peak_bx_no_ends = z_peak_bx_no_ends;
field.z_peak_by_no_ends = z_peak_by_no_ends;
field.z_peak = z_peak;
field.z_peak_no_ends = z_peak_no_ends;
field.z_zero_bx = z_zero_bx;
field.z_zero_by = z_zero_by;
field.z_zero_bx_no_ends = z_zero_bx_no_ends;
field.z_zero_by_no_ends = z_zero_by_no_ends;
field.z_zero = z_zero;
field.z_zero_no_ends = z_zero_no_ends;
field.bx_peak = bx_peak;
field.by_peak = by_peak;
field.bx_peak_no_ends = bx_peak_no_ends;
field.by_peak_no_ends = by_peak_no_ends;
field.lambda_u = lambda_u;
field.int_bx_dz = int_bx_dz;
field.int_by_dz = int_by_dz;
field.int_int_bx_dz_dz = int_int_bx_dz_dz;
field.int_int_by_dz_dz = int_int_by_dz_dz;
field.bx_filt = bx_filt;
field.by_filt = by_filt;
field.int_bx_dz_filt = int_bx_dz_filt;
field.int_by_dz_filt = int_by_dz_filt;
field.int_int_bx_dz_dz_filt = int_int_bx_dz_dz_filt;
field.int_int_by_dz_dz_filt = int_int_by_dz_dz_filt;
field.temp = temp;

%Done
