clc
clear all

%set up field grid

Nsources = 35;
Nx = 200;
Ny = 200;
xmin = -10;
xmax = 10;
ymin = -10;
ymax = 10;
xmins = -5;
xmaxs = 5;
ymins = -5;
ymaxs = 5;

x = zeros(Nx,1);
y = zeros(Ny,1);
xx = zeros(Nx,Ny);
yy = zeros(Nx,Ny);

for I = 1:Nx
    x(I) = (xmax-xmin)/(Nx-1)*I + xmin;
    for J = 1:Ny
        y(J) = (ymax-ymin)/(Ny-1)*J + ymin;
        xx(I,J) = x(I);
        yy(I,J) = y(J);
    end
end

muf0 = [0 1 0];

%set up random source origins
for I = 1:Nsources
    rsrcx = rand()*(xmaxs-xmins)+xmins;
    rsrcy = rand()*(ymaxs-ymins)+ymins;
    rs(I,:) = [rsrcx,rsrcy,0];
    mufs(I,:) = muf0;
end

%populate field

rcutoff = 0.7;
BFx = zeros(Nx,Ny);
BFy = zeros(Nx,Ny);
BFz = zeros(Nx,Ny);
for I = 1:Nx
    for J = 1:Ny
        rr = [x(I),y(J),0];
        for K = 1:Nsources
            if(norm(rr-rs(K,:))>rcutoff)
                B = infdipole(rr,rs(K,:),mufs(K,:));
                BFx(I,J) = BFx(I,J) + B(1);
                BFy(I,J) = BFy(I,J) + B(2);
                BFz(I,J) = BFz(I,J) + B(3);
            end
        end
    end
end

%come up with streamline starting coords
Nstreams = 15;
for I = 1:Nstreams
    xst = (xmax-xmin)/(Nstreams-1)*I+xmin;
    rstx(I) = xst;
    rsty(I) = ymax;
end

strms = stream2(x,y,BFx,BFy,rstx,rsty);


figure(1)
quiver(xx,yy,BFx,BFy);

figure(2)
% for I = 1:size(strms)
%     hold on
%     plot(strms{I});
% end
contourf(xx,yy,cast(BFy>0,'double'));
