%Magnetfeld Berechnung mit Biot-Savart Gesetz

tic
clear all; clc; clf

skalierungsfaktor = 10^-6;                  % vom m-Bereich zum mm-Bereich wg. T = Vs / m^2

I = 1 ;                                     % in A

radius = 3.5;                               % Spulenradius in mm

konstante = 10^-10;                         % mu0 / (4 pi) in Vs / (A mm) % Faktor 10^3 kleiner wg. mm statt m

windungen = 1;                              %Windungen


%spulenlaenge = 0;                           % Spulenlaenge in mm
%steigung = spulenlaenge / windungen

spulenlaenge = 2*pi*radius; % Spulenlaenge in mm
steigung = 0;

matrix = 30;
B = zeros(200,200,200,5);

zaehler = 0;                                % für Zeitmessung
xmax = matrix;
ymax = matrix;
zmax = 1;

genauigkeit = 5; % 1 = 6.28 Elemente pro Kreis; 2 = 12.56 Elemente pro Kreis; 4 bis 5 scheint gut zu sein
leiterelemente = ceil(2 * pi * genauigkeit * windungen) % 2 * Pi * genauigkeit für eine Umrundung
leiter = repmat(leiterelemente+1,3);

windungen = leiterelemente / genauigkeit / 2 / pi;


%for i = 1:leiterelemente+1;
 %   leiter(i,1) = i * steigung / (genauigkeit * 2 * pi) + matrix/2 - spulenlaenge/2;  % x-Ausrichtung
 %  leiter(i,2) = radius * cos(i/genauigkeit) + matrix/2;  % y-Ausrichtung
 %   leiter(i,3) = radius * sin(i/genauigkeit);   % z-Ausrichtung
%end

ii = 1:leiterelemente+1;
leiter(ii,1) = ii * steigung / (genauigkeit * 2 * 3.1415297) + matrix/2 - spulenlaenge/2;  % x-Ausrichtung
leiter(ii,2) = radius * cos(ii/genauigkeit) + matrix/2;                                    % y-Ausrichtung
leiter(ii,3) = radius * sin(ii/genauigkeit);                                               % z-Ausrichtung

for x = 1:xmax
zaehler = zaehler + 1;                                                                     % für Zeitmessung
hhh = waitbar(0,num2str(zaehler*100/matrix));                                              %  Wartebalken
waitbar(zaehler/matrix)                                                                    % Wartebalken
for y = 1:ymax                                                                             % wenn streamslice nicht genutzt wird, nur einen y-Wert berechnen
    for z = 1:zmax
        for i = 1:leiterelemente
            dl(1) = leiter(i+1,1)-leiter(i,1);
            dl(2) = leiter(i+1,2)-leiter(i,2);
            dl(3) = leiter(i+1,3)-leiter(i,3);
            vecs = [(leiter(i,1)+leiter(i+1,1))/2, ...
                (leiter(i,2)+leiter(i+1,2))/2, ...
                (leiter(i,3)+leiter(i+1,3))/2];
            vecr = [x y z];
            vecrminusvecs = vecr - vecs;
            einheitsvecr = vecrminusvecs./norm(vecrminusvecs);                            % ok
           r = sqrt(vecrminusvecs(1).^2 + vecrminusvecs(2).^2 + vecrminusvecs(3).^2);     % ok
           %r = sum(vecrminusvecs.^2);
           vektorprodukt = [dl(2).*einheitsvecr(3) - dl(3).*einheitsvecr(2), ...
                dl(3).*einheitsvecr(1) - dl(1).*einheitsvecr(3), ...
                dl(1).*einheitsvecr(2) - dl(2).*einheitsvecr(1)];
            dB = konstante * I * vektorprodukt / (r.^2);
           
            dB = dB / skalierungsfaktor;                       % nur hier wird der Wert verändert bzw. skaliert
            
            B(x,y,z,1) = B(x,y,z,1) + dB(1);
            B(x,y,z,2) = B(x,y,z,2) + dB(2);
            B(x,y,z,3) = B(x,y,z,3) + dB(3);
            %B(x,y,z,1:3)=B(x,y,z,1:3)+dB(1:3);
            
            B(x,y,z,4) = B(x,y,z,4) + sqrt(dB(1).^2 + dB(2).^2 + dB(3).^2);
        end;
    end;
end;
close(hhh) 
end;
toc
n = 1:leiterelemente;
Lx = leiter(n,1);
Ly = leiter(n,2);
Lz = leiter(n,3);
%subplot(2,1,2), 
line(Lx,Ly,Lz,'Color','k','LineWidth',2); 
hold on
view(15,30);            % view(0,0) = Blickwinkel, 2D-Perspektive
grid on                 % Gitter anzeigen
xlim([0 matrix])
ylim([0 matrix])
zlim([0 5])
xlabel('x-Achse');
ylabel('y-Achse');
zlabel('z-Achse');
daspect([1 1 1])
[X,Y]=meshgrid(1:matrix);
U=(B(1:matrix,1:matrix,z,1))'; 
V=(B(1:matrix,1:matrix,z,2))';
streamslice(X,Y,U,V) % quiver, streamslice#

