%**************************************************************************
%This function models the shimmed phase error giving the expected results of
%the shimming.
%
%Input:
%plt_file, plot file
%param, parameter structure
%field, magnetic field structure
%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
%
%Zachary Wolf
%1/15/06

function [] = pid_shim_model_phaserr(plt_file, param, field, phaserr_shim)

%Extract the required field information
z = field.z;
bx = field.bx;
by = field.by;
nom_lambda_u = param.nom_lambda_u;
gamma = param.gamma;
num_poles = param.num_poles;
num_end_poles = param.num_end_poles;
nom_k_value = param.nom_k_value;
cell_length = param.cell_length;
z_peak = field.z_peak;

%Model the shims, add the shim fields to the measured fields
bx_m = bx;
by_m = model_phaserr_shim(nom_lambda_u, phaserr_shim, z, by);

%Calculate the field integrals
int_bx_dz_m = int_calc_int(z, bx_m);
int_by_dz_m = int_calc_int(z, by_m);
int_int_bx_dz_dz_m = int_calc_int(z, int_bx_dz_m);
int_int_by_dz_dz_m = int_calc_int(z, int_by_dz_m);

%Calculate trajectories
[x_m, x_p_m, x_pp_m] = traj_calc_x_traj(gamma, by_m, int_by_dz_m, int_int_by_dz_dz_m);
[y_m, y_p_m, y_pp_m] = traj_calc_y_traj(gamma, bx_m, int_bx_dz_m, int_int_bx_dz_dz_m);

%Calculate the slippage
slip = slip_calc_slip(gamma, z, x_p_m, y_p_m);

%Find the slippage at the field peaks
slip_peak = slip_calc_slip_at_peaks(z, slip, z_peak);

%Perform a linear fit to the slippage at the field peaks, neglect ends
z_peak_no_ends = z_peak(num_end_poles + 1 : length(z_peak) - num_end_poles);
slip_peak_no_ends = slip_peak(num_end_poles + 1 : length(slip_peak) - num_end_poles);
[slope, intercept, slip_resid_peak_no_ends] = slip_fit_slip_at_peaks(z_peak_no_ends, slip_peak_no_ends);

%Calculate Keff
k_eff = k_calc_keff(gamma, slope);

%Calculate Beff
b_eff = k_calc_beff(nom_lambda_u, k_eff);

%Calculate the radiation wavelength
lambda_r = k_calc_lambdar(nom_lambda_u, slope);

%Calculate the phase as a function of z, phase(z(1)) = 0 ref
phase = phase_calc_phase(lambda_r, slip);

%Calculate the phase at the field peaks
phase_peak = phase_calc_phase(lambda_r, slip_peak);

%Calculate the phase at the field peaks, neglect the ends
phase_peak_no_ends = phase_calc_phase(lambda_r, slip_peak_no_ends);

%Calculate the phase at the cell entrance
phase_cell_entrance = phase_calc_phase_cell_entrance(lambda_r, gamma, cell_length, z, phase);

%Calculate the phase at the cell exit
phase_cell_exit = phase_calc_phase_cell_exit(lambda_r, gamma, cell_length, z, phase);

%Make the cell entrance phase = 0 reference
phase = phase - phase_cell_entrance;
phase_peak = phase_peak - phase_cell_entrance;
phase_peak_no_ends = phase_peak_no_ends - phase_cell_entrance;
phase_cell_exit = phase_cell_exit  - phase_cell_entrance;

%Calculate the phase change through the cell
phase_cell = phase_cell_exit;

%Calculate the phase errors at the field peaks, neglect ends
phase_error_peak_no_ends = phase_calc_phase_error(lambda_r, slip_resid_peak_no_ends);

%Plot the modeled phase errors after shimming
model_phaserr_shim_plot_phaserr(plt_file, param, z_peak_no_ends, phase_error_peak_no_ends);

%Done
