%**************************************************************************
%magic_fingers_sim
%This program simulates the magnetic field from a magic fingers array.
%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.
%
%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_num_x_pos        %Number of x positions (ie. upper blocks) in the magic finger array
%   pm_x_dis            %X distance between the centers of the magnet blocks (m)
%   pm_y_dis(1 : pm_num_x_pos)  %Y distance from each magnet block surface to the midplane (m)
%   pm_z_pos            %Z position of the magnet array (m)
%   pm_m_dir_up(1 : pm_num_x_pos)  %Magnetization direction of the upper blocks, Convention: +x = 1, -x = 2, +y = 3, -y = 4
%   pm_m_dir_lo(1 : pm_num_x_pos)  %Magnetization direction of the lower blocks, Convention: +x = 1, -x = 2, +y = 3, -y = 4
%   pm_div_size         %Permanent magnet segmentation size along each axis for equivalent charges (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)
%
%field_int, structure containing the calculated field integrals
%   pos(nx2)            %Position [x(n), y(n)] of calculated field integrals in global coordinates (m)
%   int_bx_dz(nx1)      %Field integral (int Bx dz) at position [x(n), y(n)] (Tm)
%   int_by_dz(nx1)      %Field integral (int By dz) at position [x(n), y(n)] (Tm)
%
%Zachary Wolf, 9/14/10

%Clear all variables
clear;

%Set the working directory
cd('c:\Simulation\Magic Fingers Sim');

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

%Program parameters
[param] = magic_fingers_sim_param();
r0 = 0.01;  %reference radius for field harmonics

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

%Initialize the data file
util_sim_dat_init(dat_file);

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

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

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

%Initialize the permanent magnet blocks
%disp('Initializing blocks...');
[block] = magic_fingers_block_init(param);

%Calculate charges on blocks
%disp('Calculating charges...');
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
%disp('Calculating fields...');
[field] = pm_field_calc_b(param.field_pos, charge);

%Calculate field integrals from the charges
%disp('Calculating field integrals...');
[field_int] = pm_field_calc_b_int(param.field_int_pos, charge);

%Debug
% clear field_int;
% field_int.pos(:, 1) = [1, 2, 3, 4, 5, 6, 7]';
% field_int.pos(:, 2) = [0, 0, 0, 0, 0, 0, 0]';
% field_int.int_bx_dz = [1, 2, 3, 4, 5, 6, 7]';
% %field_int.int_by_dz = [4, 4, 4, 4, 4, 4, 4]';
% field_int.int_by_dz = field_int.pos(:, 1).^2;

%Calculate the multipole fields
[al_n, bl_n] = magic_fingers_calc_multipoles(field_int, r0);

%Plot the field
magic_fingers_plot_bv_vs_xy(plt_file, block, charge, field, 0);
magic_fingers_plot_bv_vs_yz(plt_file, block, charge, field, 0);
magic_fingers_plot_bx_vs_z(plt_file, block, charge, field, 0, 0);
magic_fingers_plot_by_vs_z(plt_file, block, charge, field, 0, 0);
magic_fingers_plot_bx_vs_x(plt_file, block, charge, field, 0, 0);
magic_fingers_plot_by_vs_x(plt_file, block, charge, field, 0, 0);

magic_fingers_plot_by_vs_z(plt_file, block, charge, field, -.06, 0);

%Plot the field integrals
magic_fingers_plot_intbxdz_vs_x(plt_file, block, charge, field_int, 0);
magic_fingers_plot_intbydz_vs_x(plt_file, block, charge, field_int, 0);

%Plot the multipole fields
magic_fingers_plot_har_intbxydz(plt_file, al_n, bl_n, r0);

%Write the results to a data file
magic_fingers_dat_har_intbxydz(dat_file, al_n, bl_n, r0);
%magic_fingers_dat_bxyz_vs_z(dat_file, field, 0, 0);
magic_fingers_dat_intbxydz_vs_x(dat_file, field_int, 0);

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

%Debug
% x_val_req = .06;
% y_val_req = 0;
% [x_val, y_val, int_bx_dz, int_by_dz] = magic_fingers_calc_num_int_bxy_dz(field, x_val_req, y_val_req);
% disp(['Debug: int_by_dz = ', num2str(int_by_dz), ' (Tm) at x_val = ', num2str(x_val), ' (m)']);

%Done
%disp('Done');
