% % Octave / MATLAB script for plotting magnetic field lines, generated by coils. % It discretizes coils in a number of discrete elements (and plot them). % Then, each element is used to compute the magnetic field DIRECTION (not the actual amplitude) at various locations. % Finally, the field lines are obtained from the magnetic field. % % NOTE: the amplitude of the magnetic field is not computed, but you can easily modify the script to compute it. % Field lines are computed using the midpoint euler integration scheme (which is 2nd order in accuracy). % To speed-up the computation, reduce Nx_probe, Ny_probe, Nele (for each cell) and . % % COPYLEFT: Stefano Boccelli - BoccelliEngineering % close all clear clc try page_screen_output(0); end % +++++++++++++++ Coils details +++++++++++++++++++ % For each coil, put one entry in the cell array coils{1}.I = 1; % [A] coil current coils{1}.Nturns = 150; % [-] Number of turns coils{1}.Nele = 50; % [-] Number of elements for discretizing the coil coils{1}.Rc = 0.012; % [m] coil radius coils{1}.x_coil = 0.0; % [m] x-position of the coil coils{2}.I = 1; coils{2}.Nturns = 150; coils{2}.Nele = 50; coils{2}.Rc = 0.012; coils{2}.x_coil = 0.1; % Create coils fprintf('Creating coils...\n') figure hold on for cID= 1:numel(coils) th_vect = linspace(0,2*pi, coils{cID}.Nele+1); xc = zeros(size(th_vect)) + coils{cID}.x_coil; yc = coils{cID}.Rc*cos(th_vect); zc = coils{cID}.Rc*sin(th_vect); coils{cID}.xc = xc; coils{cID}.yc = yc; coils{cID}.zc = zc; coils{cID}.ri_vec = [(xc(2:end)+xc(1:end-1))/2; (yc(2:end) + yc(1:end-1))/2; (zc(2:end) + zc(1:end-1))/2]; % Wire elements center coils{cID}.Li_vec = [xc(2:end)-xc(1:end-1); yc(2:end) - yc(1:end-1); zc(2:end) - zc(1:end-1)]; % Wire elements vectors plot3(xc, yc, zc, 'k', 'linewidth',2) xlabel('x [m]') ylabel('y [m]') zlabel('z [m]') axis equal view([60,30]) end % +++++++ Create domain for computing the magnetic field (x,y plane) ++++++++++++++++ % Domain limits x_min_probe = -0.05; x_max_probe = 0.15; y_min_probe = -0.075; y_max_probe = 0.075; % Number of elements Nx_probe = 50; Ny_probe = 50; x_probe = linspace(x_min_probe, x_max_probe, Nx_probe); y_probe = linspace(y_min_probe, y_max_probe, Ny_probe); [XX, YY] = meshgrid(x_probe, y_probe); % Compute magnetic field BBx = zeros(size(XX)); BBy = zeros(size(XX)); BBz = zeros(size(XX)); % Loop on positions and compute the magnetic field for ii = 1:Nx_probe fprintf('Computing magnetic field - step %d of %d\n', ii, Nx_probe) for jj = 1:Ny_probe % Current position r_probe = [x_probe(ii); y_probe(jj); 0]; % Loop on coils for cID = 1:numel(coils) % Loop on coil elements for kk = 1:coils{cID}.Nele r_k = coils{cID}.ri_vec(:,kk); % Position vector of wire element % The constant in front of the expression is mu_0 / 4 pi. TOLER = 1.0e-5; % Tolerance (see dist_k = ...) dist_k = max(norm(r_probe - r_k), TOLER); % Avoid computing elements that are closer than tolerance TOLER, or computation blows up. % Inaccurate very close to the wires. B_vec_k = (4*pi*1e-7)/(4*pi)*coils{cID}.I*coils{cID}.Nturns*cross(coils{cID}.Li_vec(:,kk), (r_probe - r_k)/( dist_k^3 )); % Update field BBx(jj,ii) = BBx(jj,ii) + B_vec_k(1); BBy(jj,ii) = BBy(jj,ii) + B_vec_k(2); BBz(jj,ii) = BBz(jj,ii) + B_vec_k(3); end end end end figure quiver(XX, YY, BBx, BBy) dx = 0.0005; % Step for streamline integration % field lines % y0_vect = [y_min_probe+0.001:0.005:y_max_probe-0.001]; step_y0 = 0.005; y0_vect = [ [y_min_probe+0.001:step_y0:0], 0, [y_max_probe-0.001:-step_y0:0] ]; for f_ID = 1:numel(y0_vect) % Loop on field lines fprintf('Plotting field line %d of %d\n', f_ID, numel(y0_vect)) x0 = (x_max_probe + x_min_probe)/2; y0 = y0_vect(f_ID); xx = x0; yy = y0; % -------- Forward integration ------------ for tt = 2:200 % Current position x_now = xx(tt-1); y_now = yy(tt-1); if ( (x_now > x_max_probe) || (x_now < x_min_probe) || (y_now > y_max_probe) || (y_now < y_min_probe) ) break end %%% % FORWARD EULER %%% % Interpolate fields at current position %%% Bx_now = interp2(XX, YY, BBx, x_now, y_now); %%% By_now = interp2(XX, YY, BBy, x_now, y_now); %%% B_now = sqrt(Bx_now^2 + By_now^2); %%% % Forward euler integration %%% xx(tt) = x_now + dx*Bx_now/B_now; %%% yy(tt) = y_now + dx*By_now/B_now; % MIDPOINT EULER (2nd order) % Interpolate fields at current position Bx_now = interp2(XX, YY, BBx, x_now, y_now); By_now = interp2(XX, YY, BBy, x_now, y_now); B_now = sqrt(Bx_now^2 + By_now^2); % Do half step x_half = x_now + dx/2*Bx_now/B_now; y_half = y_now + dx/2*By_now/B_now; % Re-compute field Bx_now = interp2(XX, YY, BBx, x_half, y_half); By_now = interp2(XX, YY, BBy, x_half, y_half); B_now = sqrt(Bx_now^2 + By_now^2); % Full step xx(tt) = x_now + dx*Bx_now/B_now; yy(tt) = y_now + dx*By_now/B_now; end hold on plot(xx, yy, 'b', 'linewidth',2) % -------- Backward integration ------------ dx = - dx; % Reverse integration direction for tt = 2:200 % Current position x_now = xx(tt-1); y_now = yy(tt-1); if ( (x_now > x_max_probe) || (x_now < x_min_probe) || (y_now > y_max_probe) || (y_now < y_min_probe) ) break end %%% % FORWARD EULER %%% % Interpolate fields at current position %%% Bx_now = interp2(XX, YY, BBx, x_now, y_now); %%% By_now = interp2(XX, YY, BBy, x_now, y_now); %%% B_now = sqrt(Bx_now^2 + By_now^2); %%% % Forward euler integration %%% xx(tt) = x_now + dx*Bx_now/B_now; %%% yy(tt) = y_now + dx*By_now/B_now; % MIDPOINT EULER (2nd order) % Interpolate fields at current position Bx_now = interp2(XX, YY, BBx, x_now, y_now); By_now = interp2(XX, YY, BBy, x_now, y_now); B_now = sqrt(Bx_now^2 + By_now^2); % Do half step x_half = x_now + dx/2*Bx_now/B_now; y_half = y_now + dx/2*By_now/B_now; % Re-compute field Bx_now = interp2(XX, YY, BBx, x_half, y_half); By_now = interp2(XX, YY, BBy, x_half, y_half); B_now = sqrt(Bx_now^2 + By_now^2); % Full step xx(tt) = x_now + dx*Bx_now/B_now; yy(tt) = y_now + dx*By_now/B_now; end hold on plot(xx, yy, 'b', 'linewidth',2) pause(0.001) end