title: Matlab: Magnetic fields calculated with biot-savarts-law
body: <p>This is the code for the computation of magnetic fields using the biot-savart-law. I hope to get some tipps for the optimization of this code. Regrettably I use german language :( I never will do this again. :)</p>
  
  <pre><code>tic
  clear all; clc; clf
  skalierungsfaktor = 10^-6; % vom m-Bereich zum mm-Bereich wg. T = Vs / m^2
  I = 1; % in A
  konstante = 10^-10; % mu0 / (4 pi) in Vs / (A mm) % Faktor 10^3 kleiner wg. mm statt m
  matrix = 30;
  B = zeros(200,200,200,5);
  zaehler = 0; % für Zeitmessung
  xmax = matrix;
  ymax = matrix;
  zmax = 1;
  radius = 9; % Spulenradius in mm
  genauigkeit = 5; % 1 = 6.28 Elemente pro Kreis; 2 = 12.56 Elemente pro Kreis; 4 bis 5 scheint gut zu sein
  windungen = 10;
  leiterelemente = ceil(2 * 3.14152 * genauigkeit * windungen) % 2 * Pi * genauigkeit für eine Umrundung
  leiter = repmat(leiterelemente+1,3);
  windungen = leiterelemente / genauigkeit / 2 / 3.1415927;
  spulenlaenge = 20; % Spulenlaenge in mm
  steigung = spulenlaenge / windungen
  for i = 1:leiterelemente+1;
      leiter(i,1) = i * steigung / (genauigkeit * 2 * 3.1415927) + 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
  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
              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,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
  </code></pre>