%**************************************************************************
%This program takes calibrated field measurements from a 2-probe array.
%It calculates the position of probe 1 (on the probe center line) relative to
%the magnetic center.  It also calculates parameters describing the
%undulaltor field.
%
%Input:
%param, structure containing program parameters
%dat_file, file to write data to
%plt_file, file to send plots to
%b_meas, structure containing measured fields
%
%Output:
%field, structure containing field analysis information
%
%Zachary Wolf
%7/3/13

function [field] = hall2_anal_field(param, dat_file, plt_file, b_meas)

%Extract the fields from the b_meas structure
z = b_meas.z;
bx1 = b_meas.bx1;
by1 = b_meas.by1;
bz1 = b_meas.bz1;
bx2 = b_meas.bx2;
by2 = b_meas.by2;
bz2 = b_meas.bz2;

%Plot the fields
halln_field_plot_bxyz_vs_z(plt_file, 'Probe 1', z, bx1, by1, bz1);
halln_field_plot_bxyz_vs_z(plt_file, 'Probe 2', z, bx2, by2, bz2);
% halln_field_plot_bxy_vs_z(plt_file, 'Probe 1', z, bx1, by1);
% halln_field_plot_bxy_vs_z(plt_file, 'Probe 2', z, bx2, by2);
% halln_field_plot_byz_vs_z(plt_file, 'Probe 1', z, by1, bz1);
% halln_field_plot_byz_vs_z(plt_file, 'Probe 2', z, by2, bz2);
% halln_field_plot_bxz_vs_z(plt_file, 'Probe 1', z, bx1, bz1);
% halln_field_plot_bxz_vs_z(plt_file, 'Probe 2', z, bx2, bz2);
hall2_field_plot_b1b2_vs_z(plt_file, 'Bx', z, bx1, bx2);
hall2_field_plot_b1b2_vs_z(plt_file, 'By', z, by1, by2);
hall2_field_plot_b1b2_vs_z(plt_file, 'Bz', z, bz1, bz2);

%Find the peaks and zeros of the fields
[field_peaks_zeros] = hall2_field_find_peaks_zeros(dat_file, plt_file, param, z, bx1, by1, bz1, bx2, by2, bz2);

%Extract field peaks and zeros
bx1_peak_ne = field_peaks_zeros.bx1_peak_ne;
by1_peak_ne = field_peaks_zeros.by1_peak_ne;
bz1_peak_ne = field_peaks_zeros.bz1_peak_ne;
bx2_peak_ne = field_peaks_zeros.bx2_peak_ne;
by2_peak_ne = field_peaks_zeros.by2_peak_ne;
bz2_peak_ne = field_peaks_zeros.bz2_peak_ne;
z_peak_bx1_ne = field_peaks_zeros.z_peak_bx1_ne;
z_peak_by1_ne = field_peaks_zeros.z_peak_by1_ne;
z_peak_bz1_ne = field_peaks_zeros.z_peak_bz1_ne;
z_peak_bx2_ne = field_peaks_zeros.z_peak_bx2_ne;
z_peak_by2_ne = field_peaks_zeros.z_peak_by2_ne;
z_peak_bz2_ne = field_peaks_zeros.z_peak_bz2_ne;
z_zero_bx1_ne = field_peaks_zeros.z_zero_bx1_ne;
z_zero_by1_ne = field_peaks_zeros.z_zero_by1_ne;
z_zero_bx2_ne = field_peaks_zeros.z_zero_bx2_ne;
z_zero_by2_ne = field_peaks_zeros.z_zero_by2_ne;

%Plot the field peaks
halln_field_plot_b12_peaks_no_ends(plt_file, 'Bx', z_peak_bx1_ne, bx1_peak_ne, z_peak_bx2_ne, bx2_peak_ne);
halln_field_plot_b12_peaks_no_ends(plt_file, 'By', z_peak_by1_ne, by1_peak_ne, z_peak_by2_ne, by2_peak_ne);
halln_field_plot_b12_peaks_no_ends(plt_file, 'Bz', z_peak_bz1_ne, bz1_peak_ne, z_peak_bz2_ne, bz2_peak_ne);

%Initialize
%Linear polarization, vertical field
dbx_at_by1_peak_ne = 0;
dby_at_by1_peak_ne = 0;
bz1_at_by1_zero_ne = 0;
bz2_at_by1_zero_ne = 0;
%Linear polarization, horizontal field
dby_at_bx1_peak_ne = 0;
dbx_at_bx1_peak_ne = 0;
bz1_at_bx1_zero_ne = 0;
bz2_at_bx1_zero_ne = 0;
%Circular polarization
bz1_at_bx1_peak_ne = 0;
bz1_at_by1_peak_ne = 0;
bz2_at_bx2_peak_ne = 0;
bz2_at_by2_peak_ne = 0;

%Linear vertical field mode-----

%Linear vertical field mode, required field values
if strcmp(param.und_pol_mode, 'LIN_POL_VERT_FIELD')
    
    %Check
    if length(by1_peak_ne) ~= length(by2_peak_ne)
        error('hall2_anal_field: probe 1 and 2 have different numbers of peaks');
    end
    
    %Find dBx at the By1 peaks
    %[bx1_at_by1_peak_ne] = halln_field_fit_peak_at_pos(z, bx1, z_peak_by1_ne);
    %[bx2_at_by2_peak_ne] = halln_field_fit_peak_at_pos(z, bx2, z_peak_by2_ne);
    %dbx_at_by1_peak_ne = bx2_at_by2_peak_ne - bx1_at_by1_peak_ne;
    %
    [bx1_at_by1_peak_ne] = halln_field_fit_field_at_pos(z, bx1, z_peak_by1_ne);
    [bx2_at_by1_peak_ne] = halln_field_fit_field_at_pos(z, bx2, z_peak_by1_ne);
    dbx_at_by1_peak_ne = bx2_at_by1_peak_ne - bx1_at_by1_peak_ne;
    %
    %dbx = bx2 - bx1;
    %[dbx_at_by1_peak_ne] = halln_field_fit_field_at_pos(z, dbx, z_peak_by1_ne);
    
    %Find dBy at the By1 peaks
    dby_at_by1_peak_ne = by2_peak_ne - by1_peak_ne;
    dby_at_by1_peak_ne = dby_at_by1_peak_ne - mean(dby_at_by1_peak_ne);  %Symmetrize the difference, corrects small offset errors
    %
    %dby = by2 - by1;
    %[dby_at_by1_peak_ne] = halln_field_fit_field_at_pos(z, dby, z_peak_by1_ne);
    %dby_at_by1_peak_ne = dby_at_by1_peak_ne - mean(dby_at_by1_peak_ne);  %Symmetrize the difference, corrects small offset errors
    
    %Find Bz at the By zeros
    %[bz1_at_by1_zero_ne] = halln_field_fit_peak_at_pos(z, bz1, z_zero_by1_ne);
    %[bz2_at_by1_zero_ne] = halln_field_fit_peak_at_pos(z, bz2, z_zero_by1_ne);
    %
    [bz1_at_by1_zero_ne] = halln_field_fit_field_at_pos(z, bz1, z_zero_by1_ne);
    [bz2_at_by1_zero_ne] = halln_field_fit_field_at_pos(z, bz2, z_zero_by1_ne);
    
    %Define a common set of z positions
    z_peak_ctr = z_peak_by1_ne;
    z_zero_ctr = z_zero_by1_ne;

    %Calculate the undulator period
    lambda_u = mean(z_zero_by1_ne(2 : end) - z_zero_by1_ne(1 : end - 1)) * 2;
    
    %Filter
    [dbx_at_by1_peak_ne_filt] = halln_field_calc_filter_alt_values(dbx_at_by1_peak_ne);
    [dby_at_by1_peak_ne_filt] = halln_field_calc_filter_alt_values(dby_at_by1_peak_ne);
    [bz1_at_by1_zero_ne_filt] = halln_field_calc_filter_alt_values(bz1_at_by1_zero_ne);
    [bz2_at_by1_zero_ne_filt] = halln_field_calc_filter_alt_values(bz2_at_by1_zero_ne);
    
    %Plot required fields
    hall2_field_plot_b1b2_db_vs_z(plt_file, 'Bx vs z, Bx at By Peak', z, bx1, bx2, by1, z_peak_ctr, bx1_at_by1_peak_ne, bx2_at_by1_peak_ne, dbx_at_by1_peak_ne, dbx_at_by1_peak_ne_filt);
    hall2_field_plot_b1b2_db_vs_z(plt_file, 'By vs z, By at Peak', z, by1, by2, by1, z_peak_ctr, by1_peak_ne, by2_peak_ne, dby_at_by1_peak_ne, dby_at_by1_peak_ne_filt);
	halln_field_plot_b_peaks_filt(plt_file, 'Probe 1 Bz, Bz at By Zeros', z, bz1, z_zero_ctr, bz1_at_by1_zero_ne, bz1_at_by1_zero_ne_filt);
	halln_field_plot_b_peaks_filt(plt_file, 'Probe 2 Bz, Bz at By Zeros', z, bz2, z_zero_ctr, bz2_at_by1_zero_ne, bz2_at_by1_zero_ne_filt);
    
    %Use filtered values
    dbx_at_by1_peak_ne = dbx_at_by1_peak_ne_filt;
    dby_at_by1_peak_ne = dby_at_by1_peak_ne_filt;
    bz1_at_by1_zero_ne = bz1_at_by1_zero_ne_filt;
    bz2_at_by1_zero_ne = bz2_at_by1_zero_ne_filt;
end

%Linear horizontal field mode-----

%Linear horizontal field mode, required field values
if strcmp(param.und_pol_mode, 'LIN_POL_HORIZ_FIELD')
    
    %Check
    if length(bx1_peak_ne) ~= length(bx2_peak_ne)
        error('hall2_anal_field: probe 1 and 2 have different numbers of peaks');
    end
    
    %Find dBy at the Bx1 peaks
    %[by1_at_bx1_peak_ne] = halln_field_fit_peak_at_pos(z, by1, z_peak_bx1_ne);
    %[by2_at_bx1_peak_ne] = halln_field_fit_peak_at_pos(z, by2, z_peak_bx1_ne);
    %dby_at_bx1_peak_ne = by2_at_bx1_peak_ne - by1_at_bx1_peak_ne;
    %
    [by1_at_bx1_peak_ne] = halln_field_fit_field_at_pos(z, by1, z_peak_bx1_ne);
    [by2_at_bx1_peak_ne] = halln_field_fit_field_at_pos(z, by2, z_peak_bx1_ne);
    dby_at_bx1_peak_ne = by2_at_bx1_peak_ne - by1_at_bx1_peak_ne;
    %
    %dby = by2 - by1;
    %[dby_at_bx1_peak_ne] = halln_field_fit_peak_at_pos(z, dby, z_peak_bx1_ne);
    
    %Find dBx at the Bx1 peaks
    dbx_at_bx1_peak_ne = bx2_peak_ne - bx1_peak_ne;
    dbx_at_bx1_peak_ne = dbx_at_bx1_peak_ne - mean(dbx_at_bx1_peak_ne);  %Symmetrize the difference, corrects small offset errors
    %
    %dbx = bx2 - bx1;
    %[dbx_at_bx1_peak_ne] = halln_field_fit_field_at_pos(z, dbx, z_peak_bx1_ne);
    %dbx_at_bx1_peak_ne = dbx_at_bx1_peak_ne - mean(dbx_at_bx1_peak_ne);  %Symmetrize the difference, corrects small offset errors
    
    %Find Bz at the Bx zeros
    %[bz1_at_bx1_zero_ne] = halln_field_fit_peak_at_pos(z, bz1, z_zero_bx1_ne);
    %[bz2_at_bx1_zero_ne] = halln_field_fit_peak_at_pos(z, bz2, z_zero_bx1_ne);
    %
    [bz1_at_bx1_zero_ne] = halln_field_fit_field_at_pos(z, bz1, z_zero_bx1_ne);
    [bz2_at_bx1_zero_ne] = halln_field_fit_field_at_pos(z, bz2, z_zero_bx1_ne);
    
    %Define a common set of z positions
    z_peak_ctr = z_peak_bx1_ne;
    z_zero_ctr = z_zero_bx1_ne;

    %Calculate the undulator period
    lambda_u = mean(z_zero_bx1_ne(2 : end) - z_zero_bx1_ne(1 : end - 1)) * 2;
    
    %Filter
    [dby_at_bx1_peak_ne_filt] = halln_field_calc_filter_alt_values(dby_at_bx1_peak_ne);
    [dbx_at_bx1_peak_ne_filt] = halln_field_calc_filter_alt_values(dbx_at_bx1_peak_ne);
    [bz1_at_bx1_zero_ne_filt] = halln_field_calc_filter_alt_values(bz1_at_bx1_zero_ne);
    [bz2_at_bx1_zero_ne_filt] = halln_field_calc_filter_alt_values(bz2_at_bx1_zero_ne);
    
    %Plot required fields
    hall2_field_plot_b1b2_db_vs_z(plt_file, 'By vs z, By at Bx Peak', z, by1, by2, bx1, z_peak_ctr, by1_at_bx1_peak_ne, by2_at_bx1_peak_ne, dby_at_bx1_peak_ne, dby_at_bx1_peak_ne_filt);
    hall2_field_plot_b1b2_db_vs_z(plt_file, 'Bx vs z, Bx at Peak', z, bx1, bx2, bx1, z_peak_ctr, bx1_peak_ne, bx2_peak_ne, dbx_at_bx1_peak_ne, dbx_at_bx1_peak_ne_filt);
	halln_field_plot_b_peaks(plt_file, 'Probe 1 Bz, Bz at Bx Zeros', z, bz1, z_zero_ctr, bz1_at_bx1_zero_ne);
	halln_field_plot_b_peaks(plt_file, 'Probe 2 Bz, Bz at Bx Zeros', z, bz2, z_zero_ctr, bz2_at_bx1_zero_ne);
    
    %Use filtered values
    dby_at_bx1_peak_ne = dby_at_bx1_peak_ne_filt;
    dbx_at_bx1_peak_ne = dbx_at_bx1_peak_ne_filt;
    bz1_at_bx1_zero_ne = bz1_at_bx1_zero_ne_filt;
    bz2_at_bx1_zero_ne = bz2_at_bx1_zero_ne_filt;
end

%Circular polarization-----

%Circular polarization, required field values
if strcmp(param.und_pol_mode, 'CIRC_POL_RIGHT_HAND') || strcmp(param.und_pol_mode, 'CIRC_POL_LEFT_HAND')
    
    %Check
    if length(bx1_peak_ne) ~= length(bx2_peak_ne)
        error('hall2_anal_field: probe 1 and 2 have different numbers of peaks');
    end
    if length(by1_peak_ne) ~= length(by2_peak_ne)
        error('hall2_anal_field: probe 1 and 2 have different numbers of peaks');
    end

    %For circular polarization use Bz at the Bx or By peak positions
    [bz1_at_bx1_peak_ne] = halln_field_fit_field_at_pos(z, bz1, z_peak_bx1_ne);
    [bz1_at_by1_peak_ne] = halln_field_fit_field_at_pos(z, bz1, z_peak_by1_ne);
    [bz2_at_bx2_peak_ne] = halln_field_fit_field_at_pos(z, bz2, z_peak_bx2_ne);
    [bz2_at_by2_peak_ne] = halln_field_fit_field_at_pos(z, bz2, z_peak_by2_ne);
    
    %Define a common set of z positions
    z_peak_ctr = z_peak_bx1_ne;

    %Calculate the undulator period
    lambda_u = mean(z_zero_bx1_ne(2 : end) - z_zero_bx1_ne(1 : end - 1)) * 2;
    
    %Filter
    [bz1_at_bx1_peak_ne_filt] = halln_field_calc_filter_alt_values(bz1_at_bx1_peak_ne);
    [bz1_at_by1_peak_ne_filt] = halln_field_calc_filter_alt_values(bz1_at_by1_peak_ne);
    [bz2_at_bx2_peak_ne_filt] = halln_field_calc_filter_alt_values(bz2_at_bx2_peak_ne);
    [bz2_at_by2_peak_ne_filt] = halln_field_calc_filter_alt_values(bz2_at_by2_peak_ne);

    %Plot the data used for the calculations
    halln_field_plot_bx12_abs_peaks_no_ends(plt_file, z_peak_bx1_ne, bx1_peak_ne, z_peak_bx2_ne, bx2_peak_ne);
    halln_field_plot_by12_abs_peaks_no_ends(plt_file, z_peak_by1_ne, by1_peak_ne, z_peak_by2_ne, by2_peak_ne);
    halln_field_plot_bz12_at_bx12_peaks_no_ends_filt(plt_file, z_peak_bx1_ne, bz1_at_bx1_peak_ne, bz1_at_bx1_peak_ne_filt, z_peak_bx2_ne, bz2_at_bx2_peak_ne, bz2_at_bx2_peak_ne_filt);
    halln_field_plot_bz12_at_by12_peaks_no_ends_filt(plt_file, z_peak_by1_ne, bz1_at_by1_peak_ne, bz1_at_by1_peak_ne_filt, z_peak_by2_ne, bz2_at_by2_peak_ne, bz2_at_by2_peak_ne_filt);
    
    %Use filtered values
    bz1_at_bx1_peak_ne = bz1_at_bx1_peak_ne_filt;
    bz1_at_by1_peak_ne = bz1_at_by1_peak_ne_filt;
    bz2_at_bx2_peak_ne = bz2_at_bx2_peak_ne_filt;
    bz2_at_by2_peak_ne = bz2_at_by2_peak_ne_filt;
    
    %Account for bx_peak and by_peak having different lengths
    len = min(length(bx1_peak_ne), length(by1_peak_ne));
    bx1_peak_ne = bx1_peak_ne(1 : len);
    by1_peak_ne = by1_peak_ne(1 : len);
    bx2_peak_ne = bx2_peak_ne(1 : len);
    by2_peak_ne = by2_peak_ne(1 : len);
    bz1_at_bx1_peak_ne = bz1_at_bx1_peak_ne(1 : len);
    bz1_at_by1_peak_ne = bz1_at_by1_peak_ne(1 : len);
    bz2_at_bx2_peak_ne = bz2_at_bx2_peak_ne(1 : len);
    bz2_at_by2_peak_ne = bz2_at_by2_peak_ne(1 : len);
    z_peak_ctr = z_peak_ctr(1 : len);
end

%Linear polarization with the field at +-45 degrees-----

%Circular polarization, required field values
if strcmp(param.und_pol_mode, 'LIN_45_DEG')
    
    %Check
    if length(bx1_peak_ne) ~= length(bx2_peak_ne)
        error('hall2_anal_field: probe 1 and 2 have different numbers of peaks');
    end
    if length(by1_peak_ne) ~= length(by2_peak_ne)
        error('hall2_anal_field: probe 1 and 2 have different numbers of peaks');
    end
    
    %Define a common set of z positions
    z_peak_ctr = z_peak_bx1_ne;

    %Calculate the undulator period
    lambda_u = mean(z_zero_bx1_ne(2 : end) - z_zero_bx1_ne(1 : end - 1)) * 2;
    
    %Account for bx_peak and by_peak having different lengths
    len = min(length(bx1_peak_ne), length(by1_peak_ne));
    bx1_peak_ne = bx1_peak_ne(1 : len);
    by1_peak_ne = by1_peak_ne(1 : len);
    bx2_peak_ne = bx2_peak_ne(1 : len);
    by2_peak_ne = by2_peak_ne(1 : len);
    z_peak_ctr = z_peak_ctr(1 : len);
end

%End of analyzing different polarization modes-----

%Return the results in the field structure
field.z = z;
field.bx1 = bx1;
field.by1 = by1;
field.bz1 = bz1;
field.bx2 = bx2;
field.by2 = by2;
field.bz2 = bz2;
field.field_peaks_zeros = field_peaks_zeros;
field.lambda_u = lambda_u;
field.z_peak_ctr = z_peak_ctr;
%Linear polarization, vertical field
field.by1_peak_ne = by1_peak_ne;
field.dbx_at_by1_peak_ne = dbx_at_by1_peak_ne;
field.dby_at_by1_peak_ne = dby_at_by1_peak_ne;
field.bz1_at_by1_zero_ne = bz1_at_by1_zero_ne;
field.bz2_at_by1_zero_ne = bz2_at_by1_zero_ne;
%Linear polarization, horizontal field
field.bx1_peak_ne = bx1_peak_ne;
field.dby_at_bx1_peak_ne = dby_at_bx1_peak_ne;
field.dbx_at_bx1_peak_ne = dbx_at_bx1_peak_ne;
field.bz1_at_bx1_zero_ne = bz1_at_bx1_zero_ne;
field.bz2_at_bx1_zero_ne = bz2_at_bx1_zero_ne;
%Circular polarization
field.bx1_peak_ne = bx1_peak_ne;
field.by1_peak_ne = by1_peak_ne;
field.bx2_peak_ne = bx2_peak_ne;
field.by2_peak_ne = by2_peak_ne;
field.bz1_at_bx1_peak_ne = bz1_at_bx1_peak_ne;
field.bz1_at_by1_peak_ne = bz1_at_by1_peak_ne;
field.bz2_at_bx2_peak_ne = bz2_at_bx2_peak_ne;
field.bz2_at_by2_peak_ne = bz2_at_by2_peak_ne;

%Done