%**************************************************************************
%epu_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)
%   a0_bx               %Fourier series coefficient a0, constant term, in Bx expansion (T)
%   an_bx               %Fourier series coefficient an, cosine terms, in Bx expansion (T)
%   bn_bx               %Fourier series coefficient bn, sine terms, in Bx expansion (T)
%   c0_bx               %Fourier series coefficient c0, constant term, in Bx magnitude and phase expansion (T)
%   cn_bx               %Fourier series coefficient cn, magnitude terms, in Bx expansion (T)
%   phin_bx             %Fourier series coefficient phin, phase terms, in Bx expansion (rad)
%   a0_by               %Fourier series coefficient a0, constant term, in By expansion (T)
%   an_by               %Fourier series coefficient an, cosine terms, in By expansion (T)
%   bn_by               %Fourier series coefficient bn, sine terms, in By expansion (T)
%   c0_by               %Fourier series coefficient c0, constant term, in By magnitude and phase expansion (T)
%   cn_by               %Fourier series coefficient cn, magnitude terms, in By expansion (T)
%   phin_by             %Fourier series coefficient phin, phase terms, in By expansion (rad)
%   deltan              %Phase of By relative to Bx (rad)
%   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] = epu_anal_field(param, dat_file, plt_file, zm, bx, by, bz, temp)

%Parameter
NUM_HAR = 20;

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

%Compensate for peaks at ends in EPU By in horizontal mode
if length(zm_peak_by) < 2 * param.num_end_poles + 6
    zm_peak_by =0;
    zm_zero_by = 0;
    by_peak = 0;
end
  
%Checks
%Number of peaks can change by 1 in and EPU
if (length(zm_peak_bx) > 1) && (abs(length(zm_peak_bx) - param.num_poles) > 1)
    warndlg(sprintf(['epu_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)) > 1)
    warndlg(sprintf(['epu_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) > 1)
    warndlg(sprintf(['epu_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)) > 1)
    warndlg(sprintf(['epu_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 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_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;
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;
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_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_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);

%Field harmonics -----

% %Debug
% lambda_u = .032;
% num_har = 10;
% z = 0 : 10.e-6 : lambda_u;
% z_start = 0;
% z_end = z(end);
% bx = 1. * cos(2 * pi / lambda_u * z) + 2. * sin(4 * pi / lambda_u * z);
% by = 3. * sin(2 * pi / lambda_u * z) + 1. * cos(4 * pi / lambda_u * z);

%Calculate the field harmonics
num_har = NUM_HAR;
if mean(abs(z_peak_bx_no_ends)) > mean(abs(z_peak_by_no_ends))
    z_start = z_peak_bx_no_ends(1);
    z_end = z_peak_bx_no_ends(end);
else
    z_start = z_peak_by_no_ends(1);
    z_end = z_peak_by_no_ends(end);
end
[a0_bx, an_bx, bn_bx] = field_calc_har(z, bx, z_start, z_end, lambda_u, num_har);
[a0_by, an_by, bn_by] = field_calc_har(z, by, z_start, z_end, lambda_u, num_har);

%Calculate the magnitudes and phases of each harmonic
c0_bx = a0_bx;
cn_bx = sqrt(an_bx.^2 + bn_bx.^2);
phin_bx = atan2(-bn_bx, an_bx);
c0_by = a0_by;
cn_by = sqrt(an_by.^2 + bn_by.^2);
phin_by = atan2(-bn_by, an_by);

%Calculate the phase by by relative to bx
deltan = phin_by - phin_bx;

% %Debug
% index_array = find((z >= z_start) & (z < z_end));
% z_fit = z(index_array);
% bx_meas = bx(index_array);
% by_meas = by(index_array);
% bx_fit = a0_bx;
% by_fit = a0_by;
% for ii = 1 : num_har
%     bx_fit = bx_fit + an_bx(ii) * cos(ii * (2 * pi / lambda_u) * z_fit) + bn_bx(ii) * sin(ii * (2 * pi / lambda_u) * z_fit);
%     by_fit = by_fit + an_by(ii) * cos(ii * (2 * pi / lambda_u) * z_fit) + bn_by(ii) * sin(ii * (2 * pi / lambda_u) * z_fit);
% end
% field_plot_bx_fit_vs_z(plt_file, z_fit, bx_meas - bx_fit);
% field_plot_by_fit_vs_z(plt_file, z_fit, by_meas - by_fit);
% bx_fit = c0_bx;
% by_fit = c0_by;
% for ii = 1 : num_har
%     bx_fit = bx_fit + cn_bx(ii) * cos(ii * (2 * pi / lambda_u) * z_fit + phin_bx(ii));
%     by_fit = by_fit + cn_by(ii) * cos(ii * (2 * pi / lambda_u) * z_fit + phin_by(ii));
% end
% field_plot_bx_fit_vs_z(plt_file, z_fit, bx_meas - bx_fit);
% field_plot_by_fit_vs_z(plt_file, z_fit, by_meas - by_fit);

%Plot the field harmonics
field_plot_har_bx(plt_file, a0_bx, an_bx, bn_bx);
field_plot_har_by(plt_file, a0_by, an_by, bn_by);

%Write the field harmonics to the data file
field_dat_har_bx(dat_file, a0_bx, an_bx, bn_bx);
field_dat_har_by(dat_file, a0_by, an_by, bn_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);

%Filters -----

%Filter width
lambda = 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_bx_int(plt_file, z, bx, bx_filt, int_bx_dz, int_bx_dz_filt, int_int_bx_dz_dz, int_int_bx_dz_dz_filt);
%int_plot_by_int(plt_file, 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_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.a0_bx = a0_bx;
field.an_bx = an_bx;
field.bn_bx = bn_bx;
field.c0_bx = c0_bx;
field.cn_bx = cn_bx;
field.phin_bx = phin_bx;
field.a0_by = a0_by;
field.an_by = an_by;
field.bn_by = bn_by;
field.c0_by = c0_by;
field.cn_by = cn_by;
field.phin_by = phin_by;
field.deltan = deltan;
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
