%**************************************************************************
%phase_anal_field
%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)
%bz, measured longitudinal field (T)
%temp, average undulator temperature (C)
%
%Output:
%field, structure containing magnetic field information
%   zm                  %Z positions at which the field measurements are made, bench coordinates (m)
%   zm_ctr              %Z position of the undulator center, bench coordinates (m)
%   z                   %Translated Z, with z=0 at undulator center (m)
%   bx                  %Measured horizontal field (T)
%   by                  %Measured vertical field (T)
%   bz                  %Measured longitudinal 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              %Z positions of the peak fields of the larger of Bx or By (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_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)
%   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_bx         %Undulator period calculated from Bx (m)
%   lambda_u_by         %Undulator period calculated from By (m)
%   lambda_u            %Undulator period, average from Bx and By (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)
%   b_thresh            %Field strength threshold for calculating the magnetic length (T)
%   mag_length_bx       %Magnetic length where |Bx|_end > b_thresh (m)
%   mag_length_by       %Magnetic length where |By|_end > b_thresh (m)
%   intbdz_thresh       %Field integral threshold for calculating the magnetic length (Tm)
%   mag_length_i1x      %Magnetic length where |I1x|_end > intbdz_thresh (m)
%   mag_length_i1y      %Magnetic length where |I1y|_end > intbdz_thresh (m)
%   temp                %Average undulator temperature (C)
%
%Zachary Wolf
%5/20/11

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

%Field -----

%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(zm, 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(zm, by);
  
%Checks
if (length(zm_peak_bx) > 1) && (abs(length(zm_peak_bx) - param.num_poles) > 0)
    warndlg(sprintf(['epu_anal_field: found wrong number of 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(['epu_anal_field: found wrong number of 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(['epu_anal_field: found wrong number of 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(['epu_anal_field: found wrong number of zeros\nlength(zm_zero_by) = ', num2str(length(zm_zero_by)), '\nnum_poles - 1 = ', num2str(param.num_poles - 1)]));
end

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

%Remove end data
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));
by_max = max(abs(by_peak));
if (bx_max > by_max)
    z_peak = z_peak_bx;
    z_peak_no_ends = z_peak_bx_no_ends;
else
    z_peak = z_peak_by;
    z_peak_no_ends = z_peak_by_no_ends;
end

%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;
else
    lambda_u = (lambda_u_bx + lambda_u_by) / 2;
end

%Calculate the magnetic length
b_thresh = 0.00004;  %(T)
%Bx
index_array = find(abs(bx) > b_thresh);
z_array = z(index_array);
if (length(z_array) > 1) && (z_array(1) > z(1)) && (z_array(end) < z(end))
    mag_length_bx = z_array(end) - z_array(1);
else
    mag_length_bx = 0;
end
%By
index_array = find(abs(by) > b_thresh);
z_array = z(index_array);
if (length(z_array) > 1) && (z_array(1) > z(1)) && (z_array(end) < z(end))
    mag_length_by = z_array(end) - z_array(1);
else
    mag_length_by = 0;
end

%Plot the field data
field_plot_bxbybz_vs_zm(plt_file, zm, bx, by, bz);
%field_plot_bxbybz_vs_z(plt_file, z, bx, by, bz);
field_plot_bx_vs_z(plt_file, z, bx);
%field_plot_bx_vs_z_pole(plt_file, z, bx, z_peak_bx);
field_plot_by_vs_z(plt_file, z, by);
%field_plot_by_vs_z_pole(plt_file, z, by, z_peak_by);
%field_plot_bz_vs_z(plt_file, z, bz);
field_plot_bxy_vs_z(plt_file, z, bx, by);
%field_plot_bxy_vs_z_pole(plt_file, z, bx, by, z_peak_bx, z_peak_by);

%Plot Bx, By with the peaks and zeros
%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 magnitudes of the peak fields excluding the ends
%field_plot_bx_abs_peaks_no_ends(plt_file, z_peak_bx, z_peak_bx_no_ends, bx_peak_no_ends);
%field_plot_by_abs_peaks_no_ends(plt_file, z_peak_by, z_peak_by_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_bx_peak_spacing_no_ends(plt_file, z_peak_bx, z_peak_bx_no_ends);
%field_plot_bx_zero_spacing_no_ends(plt_file, z_zero_bx, z_zero_bx_no_ends);
%field_plot_bx_peak_zero_spacing_no_ends(plt_file, z_peak_bx, z_peak_bx_no_ends, z_zero_bx, z_zero_bx_no_ends);
%field_plot_by_peak_spacing_no_ends(plt_file, z_peak_by, z_peak_by_no_ends);
%field_plot_by_zero_spacing_no_ends(plt_file, z_zero_by, z_zero_by_no_ends);
%field_plot_by_peak_zero_spacing_no_ends(plt_file, z_peak_by, z_peak_by_no_ends, z_zero_by, z_zero_by_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_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);

%Write the magnetic length to the data file
field_dat_magnetic_length(dat_file, b_thresh, mag_length_bx, mag_length_by);

%Field integrals -----

%Calculate the first integral of Bx and By
int_bx_dz = int_calc_int_y_dx(z, bx);
int_by_dz = int_calc_int_y_dx(z, by);
    
%Calculate the second integral of Bx and By
int_int_bx_dz_dz = int_calc_int_y_dx(z, int_bx_dz);
int_int_by_dz_dz = int_calc_int_y_dx(z, int_by_dz);

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

%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);
field_plot_i2xi2y(plt_file, z, int_int_bx_dz_dz, int_int_by_dz_dz);

%Calculate the magnetic length using the first field integrals
intbdz_thresh = 3.e-6;  %(Tm)
%Bx
index_array = find(abs(int_bx_dz) > intbdz_thresh);
z_array_enter = z(index_array);
index_array = find(abs(int_bx_dz - int_bx_dz(end)) > intbdz_thresh);
z_array_exit = z(index_array);
if (length(z_array_enter) > 1) && (length(z_array_exit) > 1) && (z_array_enter(1) > z(1)) && (z_array_exit(end) < z(end))
    mag_length_i1x = z_array_exit(end) - z_array_enter(1);
else
    mag_length_i1x = 0;
end
%By
index_array = find(abs(int_by_dz) > intbdz_thresh);
z_array_enter = z(index_array);
index_array = find(abs(int_by_dz - int_by_dz(end)) > intbdz_thresh);
z_array_exit = z(index_array);
if (length(z_array_enter) > 1) && (length(z_array_exit) > 1) && (z_array_enter(1) > z(1)) && (z_array_exit(end) < z(end))
    mag_length_i1y = z_array_exit(end) - z_array_enter(1);
else
    mag_length_i1y = 0;
end

%Write the magnetic length to the data file
field_dat_field_integral_length(dat_file, intbdz_thresh, mag_length_i1x, mag_length_i1y);

%Write the overall field integrals to the data file
int_dat_bxby_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.zm_ctr = zm_ctr;
field.z = z;
field.bx = bx;
field.by = by;
field.bz = bz;
field.z_peak_bx = z_peak_bx;
field.z_peak_by = z_peak_by;
field.z_peak = z_peak;
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_no_ends = z_peak_no_ends;
field.z_zero_bx = z_zero_bx;
field.z_zero_by = z_zero_by;
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_bx = lambda_u_bx;
field.lambda_u_by = lambda_u_by;
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.b_thresh = b_thresh;
field.mag_length_bx = mag_length_bx;
field.mag_length_by = mag_length_by;
field.intbdz_thresh = intbdz_thresh;
field.mag_length_i1x = mag_length_i1x;
field.mag_length_i1y = mag_length_i1y;
field.temp = temp;

%Done
