Implementation of the analytical expression for the magnetic field of a circular current loop and interpreta​tion/repre​sentation of the results

I want to analytically approximate the magnetic field of a few coil arrangements. For this purpose i found a very helpful paper: 20140002333.pdf (nasa.gov). On page 8 of the PDF document are the analytic expressions for the field components of the magnetic field in spherical coordinates:

This is my implementation:

function [B_r,B_theta] = magneticField_circularCoil(I,N,a,r,theta)
%MAGNETICFIELDCOMPONENTS Calculates the magnetic field components B_r and
%B_theta (spherical coordinates)
%   B_r: B component in r direction
%   B_theta: B component in theta direction
%   I: current through conductor
%   N: number of coil windings
%   a: radius of the coil
%   r: distance from the origin (spherical coordinates)
%   theta: angle to z-axis (spherical coordinates) IN DEGREES
%
%   Source for used analytic formula: 
%   https://ntrs.nasa.gov/api/citations/20140002333/downloads/20140002333.pdf

    mu0 = 4.*pi.*1e-7;
    
    alpha2 = a.^2 + r.^2 - 2.*a.*r.*sind(theta);
    beta2  = a.^2 + r.^2 + 2.*a.*r.*sind(theta);
    k2     = 1 - alpha2./beta2;
    C      = mu0 * I./pi;
    
    [K_k2,E_k2] = ellipke(k2);
    
    B_r = N.*(C.*a.^2.*cosd(theta))./(alpha2.*sqrt(beta2)) .* E_k2;
    B_theta = N.*C./(2.*alpha2.*sqrt(beta2).*sind(theta)) .* ((r.^2+a.^2.*cosd(2.*theta)).*E_k2 - alpha2.*K_k2);
    B_phi = 0;

end

To test the function, I wrote the following code:

%% Analytical calculation of the magnetic field of the Helmholtz coil arrangement %%
% Approximation: The coil diameter is neglected. All windings "in one
% place"
% approximation: The magnetic table top is assumed to act as a perfect 
% magnetic "mirror" is assumed.
%

format compact;

% Radius of the Coil in meters:
a = 0.2;
% Current through Coil in amperes:
I = 5.0;
% Number of Coil windings:
N = 154; % source: datasheet Helmholtz coils

r_test = sqrt(0.2.^2+0.2.^2);

[B_r1,B_theta1] = magneticField_circularCoil(I,N,a,r_test,45.0)
[B_r2,B_theta2] = magneticField_circularCoil(I,N,a,r_test,135.0)

This leads to the following expected results (the magnitude of the resulting field is the same but it's in different directions):

>> Magnetfeld_Helmholtzspule_analytisch
B_r1 =
5.2273e-04
B_theta1 =
-4.2355e-06
B_r2 =
-5.2273e-04
B_theta2 =
-4.2355e-06

Is my implementation of the field components correct?

And how could I represent the superimposed field of two (or more) coils? I would appreciate any ideas!

[Alan Stevens]

You define beta2 as a.^2 + r.^2 + 2.*a.*r, but the printed text has the 2ar term multiplied by sin(theta) in a similar manner to that for alpha2. (I've no idea if that will solve your problem though!).

[Florian G]

Thanks, I updated the code accordingly!

Answer 1

[Alan Stevens]

When I run it I get different signs as well as slightly different magnitudes:

%% Analytical calculation of the magnetic field of the Helmholtz coil arrangement %%
% Approximation: The coil diameter is neglected. All windings "in one
% place"
% approximation: The magnetic table top is assumed to act as a perfect 
% magnetic "mirror" is assumed.
%

format compact;

% Radius of the Coil in meters:
a = 0.2;
% Current through Coil in amperes:
I = 5.0;
% Number of Coil windings:
N = 154; % source: datasheet Helmholtz coils

r_test = sqrt(0.2.^2+0.2.^2);

[B_r1,B_theta1] = magneticField_circularCoil(I,N,a,r_test,45.0)

B_r1 = 5.7391e-04

B_theta1 = 4.8589e-05

[B_r2,B_theta2] = magneticField_circularCoil(I,N,a,r_test,135.0)

B_r2 = -5.7391e-04

B_theta2 = 4.8589e-05

function [B_r,B_theta] = magneticField_circularCoil(I,N,a,r,theta)
%MAGNETICFIELDCOMPONENTS Calculates the magnetic field components B_r and
%B_theta (spherical coordinates)
%   B_r: B component in r direction
%   B_theta: B component in theta direction
%   I: current through conductor
%   N: number of coil windings
%   a: radius of the coil
%   r: distance from the origin (spherical coordinates)
%   theta: angle to z-axis (spherical coordinates) IN DEGREES
%
%   Source for used analytic formula: 
%   https://ntrs.nasa.gov/api/citations/20140002333/downloads/20140002333.pdf

    mu0 = 4.*pi.*1e-7;
    
    alpha2 = a.^2 + r.^2 - 2.*a.*r.*sind(theta);
    beta2  = a.^2 + r.^2 + 2.*a.*r.*sind(theta);
    k2     = 1 - alpha2./beta2;
    C      = mu0 * I./pi;
    
    [K_k2,E_k2] = ellipke(k2);
    
    B_r = N.*(C.*a.^2.*cosd(theta))./(alpha2.*sqrt(beta2)) .* E_k2;
    B_theta = N.*C./(2.*alpha2.*sqrt(beta2).*sind(theta)) .* ((r.^2+a.^2.*cosd(2.*theta)).*E_k2 - alpha2.*K_k2);
    B_phi = 0;

end

[Florian G]

Yes, you are right. I resolved the issue with your previous comment. Thank you very much!