clear all
clc
fprintf('\n\nDemonstration of Jiles-Atherton model - three hysteresis loops.');
Ms=1.6e6;
a=1100;
alpha=1.6e-3;
k=400;
c=0.2; % Parameters of Jiles-Atherton model
fprintf('Calculation for parameters: \na=%f(A/m), k=%f(A/m), c=%f, Ms=%e(A/m), alpha=%e \n\n',a,k,c,Ms,alpha);
H=[0:0.01:1 1:-0.01:-1 -1:0.01:1]'; % magnetizing field H - column vector
H=[H.*1000 H.*3000 H.*6000];
fprintf('Calculations started. Expected time of calculations: below 2 minutes.\n');
Bmodel = func_ja(a,k,c,Ms,alpha,H);
fprintf('Calculations finished.\n\n');
plot(H, Bmodel,'-o');
xlabel('H (A/m)');ylabel('B (T)');grid;
function BsimT = func_ja(a,k,c,Ms,alpha,HmeasT)
% INPUT:
% a - quantifies domain density, A/m (scalar)
% k - quantifies average energy required to break pinning side, A/m (scalar)
% c - magnetization reversibility, 0..1 (scalar)
% Ms - saturation magnetization, A/m (scalar)
% alpha - Bloch coefficient (scalar)
% H - magnetizing field, A/m (column vector)
% M0 - initial value of magnetization for ODE, A/m (scalar)
% OUTPUT:
% Hw - set of output values of magnetizing field H, A/m (vector)
% Bw - set of output values of flux density B, T (vector)
if ((a<0) || (k<0) || (c<0) || (c>1) || (alpha<0))
fprintf('\n\n*** ERROR in JAn_loops: unphysical parameters.\n\n');
BsimT=zeros(size(HmeasT));
return
end % check if parameters are physical
BsimT=zeros(size(HmeasT)); % Prepare output variable
InitialCurve=0;
if sum(HmeasT(1,:))==0
InitialCurve=1;
end
M0=0; % sample demagnetized at the beginning
for lT=1:size(HmeasT,2)
H=HmeasT(:,lT);
if InitialCurve==0
H=[0; H];
end % check if there is initial magnetization curve in H
[Hw,Bw] = func_jaloop(a,k,c,Ms,alpha,H,M0);
if InitialCurve==0
Hw(1)=[];
Bw(1)=[];
H(1)=[];
Bw=Bw+linspace(0,Bw(1)-Bw(numel(Bw)),numel(Bw))'; % Close hysteresis loop
Bw=Bw-(max(Bw)+min(Bw))./2 ; % Symmetrize Bw
end
BsimT(:,lT)=Bw; % Remember the result
end
end
function [Hw,Bw] = func_jaloop(a,k,c,Ms,alpha,H,M0)
mi0=4.*pi.*1e-7;
ip=1;
ik=2;
Hw=zeros(size(H));
Hw(1)=H(1);
Mw=zeros(size(H));
Mw(1)=M0;
while ik<numel(H)
if H(ik)>H(ip)
if H(ik+1)>=H(ik)
ik=ik+1;
else
[Hi, Mi]=func_jasolver(a,k,c,Ms,alpha,H(ip),H(ik),M0);
Hi(isnan(Hi))=1;
Mi(isnan(Mi))=1;
if numel(Hi)>2
Mw(ip+1:ik)=interp1(Hi,Mi,H(ip+1:ik),'pchip');
else
Mw(ip+1)=Mi(end);
end
Hw(ip+1:ik)=H(ip+1:ik);
M0=Mi(end);
ip=ik;
ik=ip+1;
end
end
if H(ik)<H(ip)
if H(ik+1)<=H(ik)
ik=ik+1;
else
[Hi, Mi]=func_jasolver(a,k,c,Ms,alpha,H(ip),H(ik),M0);
Hi(isnan(Hi))=1;
Mi(isnan(Mi))=1;
if numel(Hi)>2
Mw(ip+1:ik)=interp1(Hi,Mi,H(ip+1:ik),'pchip');
else
Mw(ip+1)=Mi(end);
end
Hw(ip+1:ik)=H(ip+1:ik);
M0=Mi(end);
ip=ik;
ik=ip+1;
end
end
if H(ik)==H(ip)
ik=ik+1;
end
end
[Hi, Mi]=func_jasolver(a,k,c,Ms,alpha,H(ip),H(ik),M0);
Hi(isnan(Hi))=1;
Mi(isnan(Mi))=1;
if numel(Hi)>2
Mw(ip+1:ik)=interp1(Hi,Mi,H(ip+1:ik),'pchip');
else
Mw(ip+1)=Mi(end);
end
Hw(ip+1:ik)=H(ip+1:ik);
M0=Mi(end);
Bw=Mw.*mi0+Hw.*mi0;
end
function [H,M] = func_jasolver(a,k,c,Ms,alpha,Hstart,Hend,M0)
options=odeset('RelTol',1e-4,'AbsTol',1e-6,'MaxStep',abs(Hend-Hstart)./10,'InitialStep',(Hend-Hstart)./10);
dMdH_=@(H,M) dMdH(a,k,c,Ms,alpha,M,H,Hstart,Hend);
[H,M] = ode45(dMdH_,[Hstart Hend],M0,options);
end
