function [Y, Z, BY, BZ, FY, FZ] = Calc_Field_Wedge_Bar(Angle, dr, N_wedge, N_bar, Y, Z)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here

load('INT_SPACE.mat');

dy = dr*sin(Angle);dh = dr*cos(Angle);
R = [0:dy:(N_wedge*dy - dy), (N_wedge*dy - dy)*ones(1,N_bar), (N_wedge*dy - dy):-dy:0];
H = 0:-dh:(-2*(N_wedge*dh-dh)-N_bar*dh - dh);
%figure, plot(R, H, '-r*');

BY = zeros(size(Y));BZ = zeros(size(Y));

for k=1:(length(R)-1),
    if ~mod(k+1, (length(R))/10), disp([num2str((k+1)/length(R)*100), '% done']);end
    y_tmp = Y;
    z_tmp = Z - (H(k)+H(k+1))/2;
    r_tmp = (R(k)+R(k+1))/2;
    
    A_tmp = 2*r_tmp*y_tmp./(r_tmp^2 + y_tmp.^2 + z_tmp.^2);
    
    BY_tmp = r_tmp*z_tmp./(r_tmp^2 + y_tmp.^2 + z_tmp.^2).^(3/2).*INT1_space(floor(A_tmp*length(A_space))+1);
    
    BZ_tmp = r_tmp^2./(r_tmp^2 + y_tmp.^2 + z_tmp.^2).^(3/2).*INT2_space(floor(A_tmp*length(A_space))+1) - ...
        r_tmp*y_tmp./(r_tmp^2 + y_tmp.^2 + z_tmp.^2).^(3/2).*INT1_space(floor(A_tmp*length(A_space))+1);
    
    BY = BY + BY_tmp*dr;
    BZ = BZ + BZ_tmp*dr; 
end

Y = [-fliplr(Y(:, 2:end)), Y];
Z = [fliplr(Z(:, 2:end)), Z];
BY = [-fliplr(BY(:, 2:end)), BY];
BZ = [fliplr(BZ(:, 2:end)), BZ];
FY = zeros(size(BY));FZ = zeros(size(BZ));

BY_Z = BY(2:end, 1:(end-1)) - BY(1:(end-1), 1:(end-1));
BY_Y = BY(1:(end-1), 2:end) - BY(1:(end-1), 1:(end-1));
BZ_Z = BZ(2:end, 1:(end-1)) - BZ(1:(end-1), 1:(end-1));
BZ_Y = BZ(1:(end-1), 2:end) - BZ(1:(end-1), 1:(end-1));

FY(1:(end-1), 1:(end-1)) = BY(1:(end-1), 1:(end-1)).*BY_Y + BZ(1:(end-1), 1:(end-1)).*BY_Z;
FZ(1:(end-1), 1:(end-1)) = BY(1:(end-1), 1:(end-1)).*BZ_Y + BZ(1:(end-1), 1:(end-1)).*BZ_Z;


end

