Current and magnetic fields/ magnetic field center of a square loops

Dear All

Please help me to create the magnetic field at the center of a coil of square loop by using the Biot-Savart law:

Down the simulation is for circular loop, can you do the same thing but for square loop instead of circular loop?

% current_loop.m
clear all
close all
clc

% current loop ============================================
a = 1 ;                          % radius of current loop
N = 115;                       % number of elements in current loop
theta = zeros(1,N);       % angle of current loop element
xC = zeros(1,N);          % xyz coordinates for point current lop element
yC = zeros(1,N); 
zC = zeros(1,N);
dtheta = 360/N;      
theta(1) = dtheta/2;
theta(end) = 360-dtheta/2;
for c = 2 : N-1
    theta(c) = (c-1)*dtheta+theta(1);
end
xC = a.*cosd(theta); yC = a.*sind(theta);
L = 2*pi*a/N;                     % length of each current loop element
Lx = L.*cosd(90+theta); Ly = L.*sind(90+theta);

clear theta 

% Detector space (xP, yP, zP) where B is calculated ===================
NP = 217;                     % Detector points NP x N                      
xPmax = 8*a;                  % Dimensions of detector space
zPmin = 1*a/4; zPmax = 8*a;

xP = linspace(-xPmax,xPmax,NP);
zP = linspace(zPmin,zPmax,NP);
[xxP zzP] = meshgrid(xP,zP);
Bx = zeros(NP,NP);By = Bx; Bz = Bx;

% Calculation of magnetic field B: sum over each current element
for c = 1 : N
   rx= xxP - xC(c);
   rz = zzP - zC(c);
   ry = yC(c);
   r = sqrt(rx.^2 + ry.^2 + rz.^2);
   r3 = r.^3;
   Bx = Bx + Ly(c).*rz./r3;
   By = By - Lx(c).*rz./r3;
Bz = Bz + Ly(c).*rx./r3;
end
B = sqrt(Bx.^2 + By.^2 + Bz.^2);
B = B./max(max(B));            % normalize B to 1

% GRAPHICS =====================================================
figure(1)
pcolor(xxP,zzP,B.^0.2);
colormap(hot)
shading interp;
axis equal;
axis([-xPmax xPmax 0 zPmax]);
xlabel('xP');ylabel('zP');
set(gca,'XTick',[-6:2:6]);
set(gca,'YTick',[0:2:6]);
rectangle('Position',[-1 0 2 0.2],'FaceColor','k');
title('Magnetic field from current loop')
colorbar

figure(2);                               
surf(xxP,zzP,B,'FaceColor','interp',...
    'EdgeColor','none',...
    'FaceLighting','phong')
daspect([1 1 1])
axis tight
view(-122,36)
camlight left
colormap(jet)
grid off
axis off
colorbar
title('Magnetic field from current loop')

% B along z-axis: Biot-Savart & approx.
B_theory = abs(1./zP.^3);
B_theory = B_theory./max(B_theory);

figure(3)
index=find(B(1,:)==1);
plot(zP,B(:,index),'b');
hold on
plot(zP,B_theory,'r');
xlabel('zP'); ylabel('B');
legend('Biot-Savart','Approx');
title('Magnetic field from current loop: xP = 0')