clear; close all; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % FinIrrSDA is a code that calculates shape and distribution anisotropy % % for irregular assemblies of strongly magnetic particles (ferrofluid- % % impregnated pores or magnetite) in a non-magnetic matrix. Input data % % for the code consists of a table with columns for particle ID, x coor-% % dinate of the particle center, y-coordinate of the particle center and% % z-coordinate of the particle center, dimension of the longest axis, % % dimension of the intermediate axis, dimension of the shortest axis of % % best-fit ellipsoid for the particle, declination of the longest axis, % % inclination of the longest axis, declination of the intermediate axis,% % inclination of the intermediate axis, declination of the short axis, % % and inclination of the short axis in this order. Units can be chosen % % freely for dimesnsional properties, as long as the same unit is used % % throughout the table. Units for angular properties are degrees. The % % unit for intrinsic susceptibility is volume-normalized SI. To switch % % between elliptical and cylindrical particles, comment/uncomment the % % calculation for self-demagnetization tensors on L48-66 or L68-98, % % respectively. % For a full description of the model refer to Biedermann, A.R., 2020, % FinIrrSDA: A 3D model for magnetic shape and distribution anisotropy % of finite irregular arrangements of particles with different sizes, % geometries, and orientations, DOI... % Please cite DOI... when using this code. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% data = input('id, coordinates, dimension, orientation: '); % particle properties pID = data(:,1); % pore/particle ID xcoord = data(:,2); % x coordinate of pore centre ycoord = data(:,3); % y coordinate of pore centre zcoord = data(:,4); % z coordinate of pore centre a = data(:,5); % longest pore axis b = data(:,6); % intermediate pore axis c = data(:,7); % shortest pore axis Da = data(:,8); % Declination of long axis Ia = data(:,9); % Inclination of long axis Db = data(:,10); % Declination of int axis Ib = data(:,11); % Inclination of int axis Dc = data(:,12); % Declination of short axis Ic = data(:,13); % Inclination of short axis tic %% assuming elliptical or cylindrical particles, calculate the demagnetization factors Na = zeros(length(pID),1); Nb = Na; Nc = Na; V = Na; % cylindrical particles (Sato & Ishii 1989 approximation) for i = 1:length(pID) if a(i) == b(i) nn = c(i)/a(i); % dimensional ratio defined by Sato & Ishii, for oblate cylinders Na(i) = (2*nn/sqrt(pi))/(2*(2*nn/sqrt(pi))+1); Nb(i) = Na(i); Nc(i) = 1/(2*(2*nn/sqrt(pi))+1); V(i) = pi/4*nn*a(i)^3; %volume of a oblate cylinder elseif b(i) == c(i) nn = a(i)/c(i); % dimensional ratio defined by Sato & Ishii, for prolate cylinders V(i) = pi/4*nn*c(i)^3; %volume of a prolate cylinder Nc(i) = (2*nn/sqrt(pi))/(2*(2*nn/sqrt(pi))+1); Nb(i) = Nc(i); Na(i) = 1/(2*(2*nn/sqrt(pi))+1); end end % % elliptical particles (Osborn, 1945) % for i = 1:length(pID) % theta = acos(c./a); % 0 <= theta <= pi/2; bc c/a is between 0 and 1 % phi = acos(b./a); % 0 <= phi <= pi/2 % k = sin(phi)./sin(theta); % alpha = asin(k); % 0<=alpha<=pi/2 % % El_F = ellipticF(theta,sqrt(k)); % elliptic integral of the first kind % El_E = ellipticE(theta,sqrt(k)); % elliptic integral of the second kind % % if a(i) == b(i) && b(i) == c(i) % Na(i) = 1/3; % Nb(i) = 1/3; % Nc(i) = 1/3; % elseif a(i) == b(i) % Nc(i) = (cos(phi(i))*cos(theta(i)))/((sin(theta(i)))^3*(cos(alpha(i)))^2) * ((sin(theta(i))*cos(phi(i)))/cos(theta(i)) - El_E(i)); % Na(i) = (1-Nc(i))/2; % Nb(i) = Na(i); % elseif b(i) == c(i) % Na(i) = (cos(phi(i))*cos(theta(i)))/((sin(theta(i)))^3*(sin(alpha(i)))^2) * (El_F(i) - El_E(i)); % Nb(i) = (1-Na(i))/2; % Nc(i) = Nb(i); % else % Na(i) = (cos(phi(i))*cos(theta(i)))/((sin(theta(i)))^3*(sin(alpha(i)))^2) * (El_F(i) - El_E(i)); % Nb(i) = (cos(phi(i))*cos(theta(i)))/((sin(theta(i)))^3*(sin(alpha(i)))^2*(cos(alpha(i)))^2) * (El_E(i) - (cos(alpha(i)))^2*El_F(i) - ((sin(alpha))^2*sin(theta)*cos(theta))/cos(phi)); % Nc(i) = (cos(phi(i))*cos(theta(i)))/((sin(theta(i)))^3*(cos(alpha(i)))^2) * ((sin(theta(i))*cos(phi(i)))/cos(theta(i)) - El_E(i)); % end % % V(i) = 1/6*pi*a(i)*b(i)*c(i); %volume of an ellipsoid (4/3*pi*abc if a,b,c are semi-axes, our a,b,c are full axes) % % end %% define intrinsic susceptibility ki = 0.367272727; % input('intrinsic susceptibility: '); %% compute apparent susceptibility kobs for each particle k_a = ki./(1+Na*ki); % in particle coordinate system k_b = ki./(1+Nb*ki); k_c = ki./(1+Nc*ki); % transform to sample coordinate system for i = 1:length(pID) vect_a = vectorfromDI(Da(i), Ia(i)); % direction of long axis vect_b = vectorfromDI(Db(i), Ib(i)); % direction of int axis vect_c = vectorfromDI(Dc(i), Ic(i)); % direction of short axis ev = [vect_a vect_b vect_c]; % eigenvector matrix kev = diag([k_a(i) k_b(i) k_c(i)]); % eigenvalue matrix k_xyz = ev*kev/ev; % susceptibility matrix in xyz coordinates k_xyz_elements(i,:) = [k_xyz(1,1) k_xyz(2,2) k_xyz(3,3) k_xyz(1,2) k_xyz(2,3) k_xyz(1,3)]; end %% compute total shape anisotropy % % equal weight for each particle % k_xyz_total_elements = sum(k_xyz_elements)/length(pID); % % disp('Total shape anisotropy, equal weight for each particle: ') % k_xyz_total = [k_xyz_total_elements(1) k_xyz_total_elements(4) k_xyz_total_elements(6) % k_xyz_total_elements(4) k_xyz_total_elements(2) k_xyz_total_elements(5) % k_xyz_total_elements(6) k_xyz_total_elements(5) k_xyz_total_elements(3)] % weighing by total pore or particle volume for i = 1:length(pID) k_xyz_elements_volume(i,:) = V(i) * k_xyz_elements(i,:); end k_xyz_total_elements_volume = sum(k_xyz_elements_volume,1)/sum(V); % disp('Total shape anisotropy, weighted by particle volume: ') k_xyz_shape = [k_xyz_total_elements_volume(1) k_xyz_total_elements_volume(4) k_xyz_total_elements_volume(6) k_xyz_total_elements_volume(4) k_xyz_total_elements_volume(2) k_xyz_total_elements_volume(5) k_xyz_total_elements_volume(6) k_xyz_total_elements_volume(5) k_xyz_total_elements_volume(3)]; kDI_xyz_shape = kDIfromtensor(k_xyz_total_elements_volume); %% shape and distribution anisotropy % calc M and secondary fields for 6 independent field directions % calc effective field at each particle and effective magnetization % add magnetizations of all particles % calc field-parallel component of magnetization % get anisotropy ellipsoid from parallel-components Hx = [1 0 0]'; Mx = calcshapedistaniso_manyiter(Hx, pID, xcoord, ycoord, zcoord, k_xyz_elements, V); Hy = [0 1 0]'; My = calcshapedistaniso_manyiter(Hy, pID, xcoord, ycoord, zcoord, k_xyz_elements, V); Hz = [0 0 1]'; Mz = calcshapedistaniso_manyiter(Hz, pID, xcoord, ycoord, zcoord, k_xyz_elements, V); Hxy = [1 1 0]'/sqrt(2); Mxy = calcshapedistaniso_manyiter(Hxy, pID, xcoord, ycoord, zcoord, k_xyz_elements, V); Hyz = [0 1 1]'/sqrt(2); Myz = calcshapedistaniso_manyiter(Hyz, pID, xcoord, ycoord, zcoord, k_xyz_elements, V); Hzx = [1 0 1]'/sqrt(2); Mzx = calcshapedistaniso_manyiter(Hzx, pID, xcoord, ycoord, zcoord, k_xyz_elements, V); % Mx, My, Mz, Mxy, Myz and Mzx are the parallel-components of magnetization if % field was applied parallel to x, y, z, xy, yz and zx % -> calculate total anisotropy tensor Mparallel = [Mx My Mz Mxy Myz Mzx]'; % design matrix for Hx Hy Hz Hxy Hyz Hzx A = [1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0.5 0.5 0 1 0 0 0 0.5 0.5 0 1 0 0.5 0 0.5 0 0 1]; k_shapedist_elements = inv(A'*A)*A'* Mparallel; % disp('Total shape and distribution anisotropy, equal weight per particle: ') k_xyz_total = [k_shapedist_elements(1) k_shapedist_elements(4) k_shapedist_elements(6) k_shapedist_elements(4) k_shapedist_elements(2) k_shapedist_elements(5) k_shapedist_elements(6) k_shapedist_elements(5) k_shapedist_elements(3)]; kDI_total = kDIfromtensor(k_shapedist_elements'); kxvect = k_xyz_total* [1 0 0]'; kyvect = k_xyz_total* [0 1 0]'; kzvect = k_xyz_total* [0 0 1]'; kx = kxvect(1); ky = kyvect(2); kz = kzvect(3); toc