%**************************************************************************
%pm_block_rect_sim
%This program simulates the magnetic field from a permanent magnet
%block.  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_div_size         %Permanent magnet segmentation size along each axis for equivalent charges (m)
%   pm_br               %Permanent magnet remnant field (T)
%   pm_m(3)             %Magnetic moment per unit volume in coordinate system tied to block [m_x, m_y, m_z] (A/m)
%   pm_dim(3)           %Magnet block dimension in each direction, in coordinate system tied to block [w_x, w_y, w_z] (m)
%   pm_rot(3x3)         %Rotation matrix, rotate the block about its center, direction cosines R_ij = e_rot_i . e_orig_j
%   pm_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)
%   field_pos(nx3)      %Positions where field is calculated, [x(n), y(n), z(n)] coordinates (m)
%
%block(i), structure array containing the rectangular permanent magnet block information
%   num                 %Number to identify the block
%   div_size            %Magnet 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)
%   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, 4/11/08

%Clear all variables
clear;

%Set the working directory
cd('c:\Simulation\PM Block Rect Sim');

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

%Program parameters
[param] = pm_block_rect_sim_param();

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

%Initialize the data file
util_sim_dat_init(dat_file);

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

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

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

%Initialize the permanent magnet block structure
block.num = 1;
block.div_size = param.pm_div_size;
block.br = param.pm_br;  
block.m = param.pm_m;
block.dim = param.pm_dim;
block.rot = param.pm_rot;
block.trans = param.pm_trans;

%Calculate magnetic charges from the block information
[charge] = pm_block_rect_calc_charge(block);

%Add radiation damage
for ii = 1 : length(charge)
    if charge(ii).pos(1)^2 + charge(ii).pos(2)^2 < 0.005^2
        charge(ii).q = 0;
    end
end

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

%Plot the field
% pm_block_rect_plot_bv_vs_xy(plt_file, block, charge, field, 0);
pm_block_rect_plot_bv_vs_yz(plt_file, block, charge, field, 0);
% pm_block_rect_plot_bv_vs_zx(plt_file, block, charge, field, 0);
% pm_block_rect_plot_bz_vs_xy(plt_file, block, charge, field, .007);
% pm_block_rect_plot_bx_vs_yz(plt_file, block, charge, field, 0);
% pm_block_rect_plot_by_vs_zx(plt_file, block, charge, field, 0);
% pm_block_rect_plot_bx_vs_x(plt_file, block, charge, field, 0, 0);
% pm_block_rect_plot_bx_vs_y(plt_file, block, charge, field, 0, 0);
% pm_block_rect_plot_bx_vs_z(plt_file, block, charge, field, 0, 0);
% pm_block_rect_plot_by_vs_x(plt_file, block, charge, field, 0, 0);
pm_block_rect_plot_by_vs_y(plt_file, block, charge, field, 0, 0);
pm_block_rect_plot_by_vs_z(plt_file, block, charge, field, 0, block.dim(2) + .001);
% pm_block_rect_plot_bz_vs_x(plt_file, block, charge, field, 0, 0);
% pm_block_rect_plot_bz_vs_y(plt_file, block, charge, field, 0, 0);
% pm_block_rect_plot_bz_vs_z(plt_file, block, charge, field, 0, 0);

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

%Done
