%**************************************************************************
%This function models the effects of the phase error gap and phase shims.
%It calculates a new array of expected field measurements.
%
%Input:
%lambda_u, undulator period (m)
%phaserr_shim, structure containing the shim information
%   num_gap_shim        %Number of gap shims
%   gap_shim_z_start    %Z position the gap shim starts at
%   gap_shim_z_end      %Z position the gap shim ends at
%   gap_shim_pole_start %Pole number the gap shim starts at
%   gap_shim_pole_end   %Pole number the gap shim ends at
%   gap_shim_slope      %Phase error slope change from the gap shim
%   gap_shim_dbrel      %Fractional field change from the gap shim
%   gap_shim_dim        %Gap size change from the shim
%   num_phase_shim      %Number of phase shims
%   phase_shim_z        %Z position of the phase shim
%   phase_shim_pole_start  %Pole number the phase shim starts at
%   phase_shim_pole_end %Pole number the phase shim ends at
%   phase_shim_delta    %Phase error change from the phase shim
%   phase_shim_dbrel    %Fractional field change from the phase shim
%   phase_shim_dim      %Dimension of the phase shim
%z, position the field was measured at (m)
%b, measured field strength (T)
%
%Output:
%b_shim, expected field strength after shims (T)
%
%Zachary Wolf, 4/12/04, 1/15/06

function b_shim = model_phaserr_shim(lambda_u, phaserr_shim, z, b)

%Compute gap shim fields
b_gap_shim = zeros(1, length(b));
for ii = 1 : phaserr_shim.num_gap_shim
    index_array = find((z > phaserr_shim.gap_shim_z_start(ii)) &(z < phaserr_shim.gap_shim_z_end(ii)));
    add_b_gap_shim = zeros(1, length(b_gap_shim));
    add_b_gap_shim(index_array) = b(index_array) * phaserr_shim.gap_shim_dbrel(ii);
    b_gap_shim = b_gap_shim + add_b_gap_shim;
end

%Compute phase shim fields
b_phase_shim = zeros(1, length(b));
for ii = 1 : phaserr_shim.num_phase_shim
    index_array = find(abs(z - phaserr_shim.phase_shim_z(ii)) < lambda_u / 2);  %add shim over lambda_u
    add_b_phase_shim = zeros(1, length(b_phase_shim));
    add_b_phase_shim(index_array) = b(index_array) * phaserr_shim.phase_shim_dbrel(ii);
    b_phase_shim = b_phase_shim + add_b_phase_shim;
end

%Compute the shimmed field
b_shim = b + b_gap_shim + b_phase_shim;

%Done
