Number of magnetic nulls in random fields

nMagneticNullsInRandomFields.m

%% Number of magnetic nulls in random fields applying conditions in paper
N     = 100000;
iReal = 0;
iNull = 0;
for i=1:N
	M        = rand(3)-0.5;
	traceM   = M(1,1)+M(2,2)+M(3,3);
	maxM     = max(abs(M(:)));
	E        = eig(M);
	maxRealE = max(abs(real(E)));
	
	if (abs(traceM/maxRealE) <= 0.4) ...
			&& (abs(traceM/maxM) <= 0.4)
		iNull = iNull + 1;
		if isreal(eig(M))
			iReal=iReal+1;
		end
	end
	
end

disp([num2str((iNull-iReal)/iNull*100,3) '% of all nulls are spiral type when analyzing fully random magnetic field']);

iReal = 0;
iNull = 0;
for i=1:N
	M        = rand(3)-0.5;
	traceM   = M(1,1)+M(2,2)+M(3,3);
	M        = M-traceM/3;            % to make B divergence free
	maxM     = max(abs(M(:)));
	E        = eig(M);
	maxRealE = max(abs(real(E)));
	
	if (abs(traceM/maxRealE) < 0.4) ...
			&& (abs(traceM/maxM) < 0.4)
		iNull = iNull + 1;
		if isreal(eig(M))
			iReal=iReal+1;
		end
	end
	
end

disp([num2str((iNull-iReal)/iNull*100,3) '% of all nulls are spiral type when analyzing fully random magnetic field satisfying divB=0']);

%% Number of magnetic nulls in random fields applying conditions in paper
N     = 100000;
iReal = 0;
iNull = 0;
for i=1:N
	M        = rand(3)-0.5;
	traceM   = M(1,1)+M(2,2)+M(3,3);
	maxM     = max(abs(M(:)));
	E        = eig(M);
	maxRealE = max(abs(real(E)));
	
	if (abs(traceM/maxRealE) <= 0.4) ...
			&& (abs(traceM/maxM) <= 0.4)
		bAtOrigo = [rand(1) rand(1) rand(1)];
		positionNull = -bAtOrigo/(M');
		if max(abs(positionNull(:))) < 1
			iNull = iNull + 1;
			if isreal(eig(M))
				iReal=iReal+1;
			end
		end
	end
	
end

disp([num2str((iNull-iReal)/iNull*100,3) '% of all nulls are spiral type when analyzing fully random magnetic field within box x,y,z=+-1']);