%**************************************************************************
%traj_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_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)
%   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)
%   temp                %Average undulator temperature (C)
%
%Zachary Wolf
%7/23/11

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

%Parameters
int_b_dz_fringe_lim = param.int_b_dz_fringe_lim;

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

%Determine the device center
if (param.nom_zm_dev_enter == 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
else
    zm_ctr = param.nom_zm_dev_ctr;
end

%Translate Z coordinates so the center is at z=0
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;
if (param.nom_zm_dev_enter == 0) 
    z_dev_enter = z_peak_by(1);
else
    z_dev_enter = param.nom_zm_dev_enter - zm_ctr;
end
if (param.nom_zm_dev_exit == 0) 
    z_dev_exit = z_peak_by(end);
else
    z_dev_exit = param.nom_zm_dev_exit - zm_ctr;
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);

%Write bench coordinates to the data file
field_dat_bench_coord(dat_file, zm_peak_by, zm_zero_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);

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

%Fringe Fields -----

%Bx
% [z_field_start, z_field_end, len_fringe_enter, len_fringe_exit] = field_calc_fringe_fields(z, bx, int_bx_dz, z_dev_enter, z_dev_exit, int_b_dz_fringe_lim);
% field_plot_fringe_enter(plt_file, 'Bx', z, bx, int_bx_dz, z_dev_enter, z_field_start);
% field_plot_fringe_exit(plt_file, 'Bx', z, bx, int_bx_dz, z_dev_exit, z_field_end);
% field_dat_fringe_fields(dat_file, 'Bx', int_b_dz_fringe_lim, z_field_start, z_field_end, z_dev_enter, z_dev_exit, len_fringe_enter, len_fringe_exit);

%By
[z_field_start, z_field_end, len_fringe_enter, len_fringe_exit] = field_calc_fringe_fields(z, by, int_by_dz, z_dev_enter, z_dev_exit, int_b_dz_fringe_lim);
field_plot_fringe_enter(plt_file, 'By', z, by, int_by_dz, z_dev_enter, z_field_start);
field_plot_fringe_exit(plt_file, 'By', z, by, int_by_dz, z_dev_exit, z_field_end);
field_dat_fringe_fields(dat_file, 'By', int_b_dz_fringe_lim, z_field_start, z_field_end, z_dev_enter, z_dev_exit, len_fringe_enter, len_fringe_exit);

%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_zero_bx = z_zero_bx;
field.z_zero_by = z_zero_by;
field.bx_peak = bx_peak;
field.by_peak = by_peak;
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.temp = temp;

%Done
