current = 1;	 					% Fixes the current amplitude to 1 A
mu_0 = 4*pi*1e-7;					% Value of mu_0
step = 0.5;							% Resolution of the grid	
nbx = 15/step;						% number of points in x direction
nby = 20/step;						% number of points in y direction
x = linspace(-7.5,7.5,nbx);	% creates equally spaced points from -7.5 to 7.5 
y = linspace(-10,10,nby);		% creates equally spaced points from -10 to 10 

[X,Y] = meshgrid(x,y);			% Creates a grid of points

nb_dl = 10/step;					% number of elements of current
CurrentLine = linspace(-2.5,2.5,nb_dl); % Line of current

B = zeros(nby,nbx,3);			% Initiliases the B vector field
Btot = zeros(nby,nbx,3);		% Initiliases the total B vector field
for i = 1 : nb_dl					% For each small element of current
   dl = [0;step;0];				% Creates the dl vector
   for k = 1 : nby				% For each y coordinates
      for l = 1 : nbx			% For each x coordinates
         %Creates R vector over the grid
         R = [X(k,l) ; Y(k,l)- CurrentLine(i); 0];          
         % Caclulates the B value due to dl over the grid
         B(k,l,:) = mu_0/(4*pi)*current*cross(dl,R)/norm(R)^3; %Calc
      end	
   end
   z = zeros(nby,nbx); 	% Initialises z Value for the origin of B (always 0)
   figure(1);			% Create a figure
   quiver3(X,Y,z,B(:,:,1),B(:,:,2),B(:,:,3));	%Plots B
   xlabel('x in meters');
   ylabel('y in meters');
   zlabel('B in Tesla');
   title('B vector in Telsa');
   drawnow;					% Ask for a redraw;
   Btot = Btot + B;		% Generated the overall Fields to date
   figure(2);					% Create a figure
   quiver3(X,Y,z,Btot(:,:,1),Btot(:,:,2),Btot(:,:,3)); %plot Btot
   xlabel('x in meters');
   ylabel('y in meters');
   zlabel('B in Tesla');
   title('B vector in Telsa');
   drawnow;	
end
% Calculation is now finished
figure(3);	
surfl(X,Y,Btot(:,:,3)); %Plots the z value of B as Bx abd By are null
xlabel('x in meters');
ylabel('y in meters');
zlabel('B in Tesla');
title('Bz surface plot');
shading interp;		% Make it nicer
colormap(pink);		% idem
figure(4);			
contour(X,Y,Btot(:,:,3));	% Now get the contours of B
xlabel('x in meters');
ylabel('y in meters');
zlabel('B in Tesla');
title('Bz Contour plot');
BzdB = 20*log(abs(Btot(:,:,3))/1e-9); % Get Bz in dB
% Same as for B 
figure(5);	
surfl(X,Y,BzdB);
shading interp;
colormap(pink);
xlabel('x in meters');
ylabel('y in meters');
zlabel('Bz in dbNT');
title('Bz surface plot');
figure(6);
contour(X,Y,BzdB);
xlabel('x in meters');
ylabel('y in meters');
zlabel('Bz in dbNT');
title('Bz contour plot');
