Magnetic fields' 3D plotting in cylindrical coordinates

[Nikolaos Vergos]

Hello! One of my diploma thesis tasks is to calculate and plot some certain
magnetic fields.

My first task was to calculate and plot the magnetic field due to a current
loop, which includes elliptic integrals calculation. In the cylindrical
coordinates' system, we get two coefficients: Br and Bz. The code I wrote for
this is the following:
--------------------------------------------------
function [Bz,Br]=m_field(i,a,z0,r0,zmin,rmin,zmax,rmax,step)
m0=4*pi*10^-7;
B0=(i*m0)/(2*a);
in1=1;
for z=(zmin-z0):step:(zmax-z0)
in2=1;
for r=(rmin-r0):step:(rmax-r0)
if (r~=0)
al=abs(r/a);
be=abs(z/a);
ga=(z/r);
q=((1+al)^2+be^2);
k=sqrt(4*al/q);
[K,E] = ellipke(k^2);
Bz(in1,in2)=B0*(1/(pi*sqrt(q)))*(E*(1-al^2-be^2)/(q-4*al)+K);
Br(in1,in2)=B0*(ga/(pi*sqrt(q)))*(E*(1+al^2+be^2)/(q-4*al)-K);
else
Bz(in1,in2)=0;
Br(in1,in2)=0;
end
Z(in1,in2)=z;
R(in1,in2)=r;
in2=in2+1;
end
in1=in1+1;
end
s=size(Bz);
for i=2:s(1)-1
for j=2:s(2)-1
if (Bz(i,j)==0 & Br(i,j)==0)
Bz(i,j)=1/6*(Bz(i-1,j-1)+Bz(i-1,j+1)+Bz(i,j-1)+Bz(i,j+1)+Bz(i+1,j-
1)+Bz(i+1,j+1));
Br(i,j)=1/6*(Br(i-1,j-1)+Br(i-1,j+1)+Br(i,j-1)+Br(i,j+1)+Br(i+1,j-
1)+Br(i+1,j+1));
end
end
end
----------------------------------------

plus a script that, given initial values, it plots the off-axis field of the current
loop in a 2D graph.

Second task was to calculate and plot the magnetic field Bf due to a straight
current wire with Ampere's law, where f is the angle of the cylindrical
coordinates' system. In order to do this I had to change to cartesian
coordinates x and y, and in the end I also plotted the magnetic field which
follows the right hand grip rule, when the thumb points at the current (z)
direction.

My problem now is that I have somehow to create one 3D plot that will
superpose the two magnetic fields: Now the perpendicular current wire is put
along the z-axis, at the r=0 place. so the final image I have toget for the
magnetic field is a torus that comes from this superposition. The (Br,Bz) due
to the current loop are the "meridians" on the torus, and the Bq lines are the
parallels and the "equator".

Yet, I have no idea on how to put these fields together at one m file (and one
script) that will be generating the whole torus, ie (r,q,z) --> (Br,Bq,Bz)

Does anyone have the slightest idea? I would really appreciate it deeply.
Thanks in advance!