[Bx,Br]=m_field(i,a,x0,r0,xmin,rmin,xmax,rmax,step)

function [Bx,Br]=m_field(i,a,x0,r0,xmin,rmin,xmax,rmax,step)
%%%%%%%%%%%%%%%%%%%%%%%
% [Bx,Br]=m_field(i,a,x0,r0,xmin,rmin,xmax,rmax)
% Bx - axial comp. of mag. field
% Br - radial comp. of mag. field
% i=coil current;
% a=coil radius;
% x0, y0 - coil coordinate
% xmin,ymin,xmax,ymax - interest domain for magnetic field
% succes!
% 
% created by
% Iulian Teliban
%%%%%%%%%%%%%%%%%%%%%%%
miu0=4*pi*10^-7;
B0=(i*miu0)/(2*a);
in1=1;
for x=(xmin-x0):step:(xmax-x0)
    in2=1;
    for r=(rmin-r0):step:(rmax-r0)
        if (r~=0)
            al=abs(r/a);
            be=abs(x/a);
            ga=(x/r);
            q=((1+al)^2+be^2);
            k=sqrt(4*al/q);
            [K,E] = ellipke(k^2);
            Bx(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
            Bx(in1,in2)=0;
            Br(in1,in2)=0;
        end
        X(in1,in2)=x;
        R(in1,in2)=r;
        in2=in2+1;
    end
    in1=in1+1;
end
s=size(Bx);
for i=2:s(1)-1
    for j=2:s(2)-1
        if (Bx(i,j)==0 & Br(i,j)==0)
            Bx(i,j)=1/6*(Bx(i-1,j-1)+Bx(i-1,j+1)+Bx(i,j-1)+Bx(i,j+1)+Bx(i+1,j-1)+Bx(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