Plotting second derivatives of magnetic fields.generated by Helmholtz Coils

Hello.

I'm relatively new to matlab, and struggling with proper indexing or formatting of loops. This code is supposed to plot the second derivative of the magnetic fields in the x direction of the central point between two coils. The nominal magnetic field strength seems to be correct, however the plot of the second deriavtive seems far from correct. I believe it comes mostly from my lack of understanding on MATLAB indexing.

If anyone has any reccomendations for better formatting or calculating those would be greatly appreciated.

% Constants
u0 = 4*pi*10e-12;% Permeability of free space in T*m/A
i = 1; % Current
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Magnetization vectors of Helmholtz coils
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

syms theta

% Initialization of arrays
bp = zeros(3, 3);
bp1 = zeros(3, 3);
bp2 = zeros(3, 3);
posdX = zeros(3, 41);
negdX = zeros(3, 41);


for X = 0:1:40    
    for dd = 1:3
        
                
            R = 1;   % Radius of coils
            D = X*.1;   % Separation Distance
            dX = D/2000; % Differential Element for second derivative
            W = .9*D;   % Arbitrary workspace

            normfactor = u0.*i./(2*R); % Normalization factor 

            % Position vectors
            p0(:,1) = W./2.*[0; 0; 1];
            % Deltas for Central Difference formula
            p0(:,2) = W./2.*[0+dX; 0; 1];
            p0(:,3) = W./2.*[0-dX; 0; 1];

            q1 = [R.*cos(theta); R.*sin(theta); D]; % Parameterization of top coil
            q2 = [R.*cos(theta); R.*sin(theta); -D]; % Parameterization of bottom coil

            dl = [-R.*sin(theta); R.*cos(theta); 0];

            % Top coil calculation
            q1p = q1 - p0(:,dd);
            p1q = p0(:,dd) - q1;
            divisor1 = norm(p1q)^3;
            sv1 = S(q1p); % Skew Symetric Matrix

            f1 = mtimes(sv1,dl./divisor1);

            fun1 = matlabFunction(f1);

            % Magnetic field vector generated by top coil
            bp1(:,X+1,dd) = u0.*i.*integral(fun1, 0,2*pi,'ArrayValued',true)./(4.*pi);



            % Bottom coil calculation
            q2p = q2 - p0(:,dd);
            p2q = p0(:,dd) - q2;
            divisor2 = norm(p2q)^3;

            sv2 = S(q2p); % Skew Symetric Matrix

            f2 = mtimes(sv2,dl./divisor2);

            fun2 = matlabFunction(f2);

            % Magnetic field vector generated by bottom coil 
            bp2(:,X+1,dd) = u0.*integral(fun2, 0,2*pi,'ArrayValued',true)./(4.*pi);

            bp = bp2+bp1;

            % Nominal magnetic field strength
            nombp(X+1) = norm(bp(:,X+1,1))./normfactor;

        end
        % Second derivative in the X direction;
        d2bpX(X+1) = norm((bp(:,:,2) + bp(:,:,3) - 2*bp(:,1,1))./(dX)^2);
    end


figure(1);
plot(0:.1:4,nombp(:))

figure(2);
plot(0:.1:4,d2bpX)

function S = S(v)
S = [0 -v(3) v(2);v(3) 0 -v(1); -v(2) v(1) 0];
end

function norm = n(v)
n = sqrt(v(1^2+v(2)^2+v(3)^2));
end

Answer 1

[Walter Roberson]

n = sqrt(v(1^2+v(2)^2+v(3)^2));

should be

n = sqrt(v(1)^2+v(2)^2+v(3)^2);

[Araxan]

Thank you, that did not fix the issue unfortunately.