%**************************************************************************
%virtual_shim_sim
%This program simulates the magnetic field from a virtual shim.  The field
%from the shimmed position has the field from the nominal position
%subtracted to get the net field of the virtual shim.
%Magnet block parameters are entered.  Equivalent charges are
%calculated to replace the magnet blocks.  The field from each charge is
%calculated and added to the field from the other charges.  Field integrals
%are calculated.
%
%The following structures are used.
%
%param, structure containing program parameters
%   output_file_dir     %Directory for output files
%   pm_br               %Remnant field in the blocks (T)
%   pm_dim(3)           %Magnet block dimensions for each block in coordinate system tied to block [w_x, w_y, w_z] (m) 
%   pm_x_nom_pos(1 : num_blocks)  %Nominal x positions of the blocks (m)
%   pm_y_nom_pos(1 : num_blocks)  %Nominal y positions of the blocks (m)
%   pm_z_nom_pos(1 : num_blocks)  %Nominal z positions of the blocks (m)
%   pm_m_dir(1 : num_blocks)      %Magnetization direction of the blocks, Convention: +x = 1, -x = 2, +y = 3, -y = 4, +z = 5, -z = 6
%   pm_div_size         %Permanent magnet segmentation size along each axis for equivalent charges (m)
%   shim_x_shift(1 : num_blocks)  %Virtual shim x shift (m)
%   shim_y_shift(1 : num_blocks)  %Virtual shim y shift (m)
%   shim_z_shift(1 : num_blocks)  %Virtual shim z shift (m)
%   field_pos(n, 3)     %Positions where field is calculated, [x(n), y(n), z(n)] coordinates (m)
%   field_int_pos(n, 2) %Positions where field integrals are calculated, [x(n), y(n)] coordinates (m)
%
%block_rect(i), structure array containing the rectangular permanent magnet block information
%   num                 %Number to identify the block
%   div_size            %Block segmentation size along each axis for equivalent charges (m)
%   br                  %Permanent magnet remnant field (T)
%   m(3)                %Magnetic moment per unit volume in coordinate system tied to block [m_x, m_y, m_z] (A/m)
%   dim(3)              %Magnet block dimension in each direction in coordinate system tied to block [w_x, w_y, w_z] (m)
%   rot(3x3)            %Rotation matrix, rotate the block about its center, direction cosines R_ij = e_rot_i . e_orig_j
%   trans(3)            %Translate the rotated block in x, y, and z in global coordinates parallel to the original block coordinates [t_x, t_y, t_z] (m)
%
%charge(i), structure containing the charge information
%   q                   %Magnetic charge, M.n dA (Am)
%   pos_block(3)        %Position of charge [x, y, z] in block coordinates (m)
%   pos(3)              %Position of charge [x, y, z] in global coordinates (m)
%   x                   %x-position of charge in global coordinates, pos(1), for vectorization (m)
%   y                   %y-position of charge in global coordinates, pos(2), for vectorization (m)
%   z                   %z-position of charge in global coordinates, pos(3), for vectorization (m)
%   block_num           %Block number of the parent block
%
%field, structure containing the calculated magnetic field
%   pos(nx3)            %Position [x(n), y(n), z(n)] of calculated field in global coordinates (m)
%   b(nx3)              %Magnetic field [Bx(n), By(n), Bz(n)] at position [x(n), y(n), z(n)] (T)
%
%Zachary Wolf, 2/4/11

%Clear all variables
clear;

%Set the working directory
cd('c:\Simulation\Virtual Shim Sim');

%Global
global mu_0;
mu_0 = 4 * pi * 10^-7;  %(Tm/A)

%Program parameters
[param] = virtual_shim_sim_param();

%Make file names
dat_file = strcat(lower(param.output_file_dir), 'virtual_shim_sim.dat');
plt_file = strcat(lower(param.output_file_dir), 'virtual_shim_sim.ps');
mat_file = strcat(lower(param.output_file_dir), 'virtual_shim_sim.mat');

%Initialize the data file
util_sim_dat_init(dat_file);

%Write a message to the data file
msg{1} = ' ';
msg{2} = 'Virtual Shim Simulation';
msg{3} = '=======================';
util_sim_dat_msg(dat_file, msg);

%Write the parameter values to the data file
virtual_shim_dat_param(dat_file, param);

%Initialize the plot file
util_sim_plt_init(plt_file, msg, param);

%Initialize the permanent magnet blocks
[block] = virtual_shim_block_init(param);

%Calculate charges on blocks
charge(length(block) * 2000) = struct('q', 0, 'pos_block', [0, 0, 0], 'pos', [0, 0, 0], 'x', 0, 'y', 0, 'z', 0, 'block_num', 0);
nn = 0;
for ii = 1 : length(block)
    [charge_block] = pm_block_rect_calc_charge(block(ii));
    for jj = 1 : length(charge_block)
        nn = nn + 1;
        charge(nn) = charge_block(jj);
    end
end
charge = charge(1 : nn);  %remove extra pre-allocated charges

%Calculate fields from the charges
[field] = pm_field_calc_b(param.field_pos, charge);

%Calculate the field integrals
[x_val, y_val, z, int_bx_dz, int_by_dz, int_int_bx_dz_dz, int_int_by_dz_dz] = pm_field_calc_int_bxy_dz(field, 0, 0);

%Plot the field
virtual_shim_plot_bv_vs_yz(plt_file, block, charge, field, 0);
virtual_shim_plot_bv_vs_xy(plt_file, block, charge, field, 0);
virtual_shim_plot_bx_vs_z(plt_file, block, charge, field, 0, 0);
virtual_shim_plot_by_vs_z(plt_file, block, charge, field, 0, 0);

%Plot the field integrals
virtual_shim_plot_intbxdz_vs_z(plt_file, block, charge, x_val, y_val, z, int_bx_dz);
virtual_shim_plot_intbydz_vs_z(plt_file, block, charge, x_val, y_val, z, int_by_dz);
virtual_shim_plot_intintbxdzdz_vs_z(plt_file, block, charge, x_val, y_val, z, int_int_bx_dz_dz);
virtual_shim_plot_intintbydzdz_vs_z(plt_file, block, charge, x_val, y_val, z, int_int_by_dz_dz);

%Write the results to a data file
%virtual_shim_dat_bxyz_vs_z(dat_file, field, 0, 0);
virtual_shim_dat_intbxydz(dat_file, x_val, y_val, z, int_bx_dz, int_by_dz, int_int_bx_dz_dz, int_int_by_dz_dz);

%Output the results to the output file directory
save(mat_file, 'param', 'block', 'charge', 'field');

%Done
disp('Done');
