%**************************************************************************
%This function plots the z component of the magnetic field in the
%x-y plane whose z position is specified.
%The closest z position with data points is used.
%
%Input:
%plt_file, plot file for results
%block, structure containing the block information
%charge, structure containing the magnetic charge information
%field, structure containing the magnetic field information
%z_val_req, requested z position of slice (m)
%
%Zachary Wolf, 4/19/10

function pm_block_rect_plot_bz_vs_xy(plt_file, block, charge, field, z_val_req)

%Initialize
figure;
hold on;

%Draw the projection of the blocks
for ii = 1 : length(block)
    p1 = [- block(ii).dim(1) / 2,   block(ii).dim(2) / 2,   block(ii).dim(3) / 2];
    p2 = [  block(ii).dim(1) / 2,   block(ii).dim(2) / 2,   block(ii).dim(3) / 2];
    p3 = [- block(ii).dim(1) / 2, - block(ii).dim(2) / 2,   block(ii).dim(3) / 2];
    p4 = [  block(ii).dim(1) / 2, - block(ii).dim(2) / 2,   block(ii).dim(3) / 2];
    p5 = [- block(ii).dim(1) / 2,   block(ii).dim(2) / 2, - block(ii).dim(3) / 2];
    p6 = [  block(ii).dim(1) / 2,   block(ii).dim(2) / 2, - block(ii).dim(3) / 2];
    p7 = [- block(ii).dim(1) / 2, - block(ii).dim(2) / 2, - block(ii).dim(3) / 2];
    p8 = [  block(ii).dim(1) / 2, - block(ii).dim(2) / 2, - block(ii).dim(3) / 2];
    p1g = p1 * block(ii).rot + block(ii).trans;
    p2g = p2 * block(ii).rot + block(ii).trans;
    p3g = p3 * block(ii).rot + block(ii).trans;
    p4g = p4 * block(ii).rot + block(ii).trans;
    p5g = p5 * block(ii).rot + block(ii).trans;
    p6g = p6 * block(ii).rot + block(ii).trans;
    p7g = p7 * block(ii).rot + block(ii).trans;
    p8g = p8 * block(ii).rot + block(ii).trans;
    line([p1g(1), p2g(1)], [p1g(2), p2g(2)]);
    line([p1g(1), p3g(1)], [p1g(2), p3g(2)]);
    line([p1g(1), p5g(1)], [p1g(2), p5g(2)]);
    line([p2g(1), p4g(1)], [p2g(2), p4g(2)]);
    line([p2g(1), p6g(1)], [p2g(2), p6g(2)]);
    line([p3g(1), p4g(1)], [p3g(2), p4g(2)]);
    line([p3g(1), p7g(1)], [p3g(2), p7g(2)]);
    line([p4g(1), p8g(1)], [p4g(2), p8g(2)]);
    line([p5g(1), p6g(1)], [p5g(2), p6g(2)]);
    line([p5g(1), p7g(1)], [p5g(2), p7g(2)]);
    line([p6g(1), p8g(1)], [p6g(2), p8g(2)]);
    line([p7g(1), p8g(1)], [p7g(2), p8g(2)]);
end

%Draw the projection of the charges
for ii = 1 : length(charge)
    if (charge(ii).q > 0)
        plot3(charge(ii).pos(1), charge(ii).pos(2), 0, 'r+');
    end
    if (charge(ii).q < 0)
        plot3(charge(ii).pos(1), charge(ii).pos(2), 0, 'kx');
    end
end

%Plot the field
[z_close, iz_close] = min(abs(field.pos(:, 3) - z_val_req));
z_val = field.pos(iz_close(1), 3);
index_array = find(field.pos(:, 3) == z_val);
x_array = field.pos(index_array, 1);
y_array = field.pos(index_array, 2);
bz_array = field.b(index_array, 3);
x_vals = unique(x_array);
y_vals = unique(y_array);
bz_vals = zeros(length(y_vals), length(x_vals));
for ii = 1 : length(bz_array)
    ix = find(x_vals == x_array(ii));
    iy = find(y_vals == y_array(ii));
    bz_vals(iy, ix) = bz_array(ii);
end
mesh(x_vals, y_vals, bz_vals);

%Labels
title(['Bz vs XY At Z = ', num2str(z_val), ' m']);
xlabel('X Position (m)');
ylabel('Y Position (m)');
zlabel('Bz (T)');

%Set axes scale
x_limits = xlim;
y_limits = ylim;
lim = max(abs([x_limits, y_limits]));
xlim([-lim, lim]);
ylim([-lim, lim]);

%Finalize
grid on;
zoom on;
hold off;

%Print
orient landscape;
print('-append', '-dpsc', plt_file);

%Done
