%{
A. J. Bergqvist, "A simple vector generalization of the Jiles-Atherton model of hysteresis," 
in IEEE Transactions on Magnetics, vol. 32, no. 5, pp. 4213-4215, Sept. 1996, doi: 10.1109/20.539337
%}

function [H, M, B, Man] = func_Bergqvist(JA, H_amp, f, nCycles, nSteps)
    % Inputs:
    % JA: Struct with fields Ms, a, k, c, alpha
    % H_amp: Amplitude of external field
    % f: Frequency [Hz]
    % nCycles: Number of cycles
    % nSteps: Steps per cycle
    
    mu0 = 4*pi*1e-7;

    delta_H = H_amp/nSteps;
    H1 = 0 : delta_H : H_amp;
    H2 = (H_amp-delta_H) : -delta_H : -H_amp;
    H3 = (-H_amp+delta_H):  delta_H : -delta_H;

    H_cycle = [H1, H2, H3];

    H = repmat(H_cycle, 1, nCycles);


    % Initialization
    M   = zeros(size(H));
    Man = zeros(size(H));
    
    % Initial value for M
    M(1) = 0;

    for i = 2:length(H)-1
        % RK4-method
        dt = H(i)-H(i-1);
        [k1, He1] = dMdt(JA, H(i),   M(i), dt);
        [k2, ~]   = dMdt(JA, H(i)+dt/2, M(i)+dt/2*k1, dt);
        [k3, ~]   = dMdt(JA, H(i)+dt/2, M(i)+dt/2*k2, dt);
        [k4, ~]   = dMdt(JA, H(i)+dt,   M(i)+dt*k3, dt);
        
        M(i+1) = M(i) + dt/6 * (k1 + 2*k2 + 2*k3 + k4);
        Man(i+1) = func_Man(JA, He1); 
    end

    B = mu0 * (H + M); 
end


function [dM, Heff] = dMdt(JA, H, M,dt)
    Heff = H + JA.alpha * M;
    
     % Take care of singularity at Heff~=0 and M~=0
    if abs(Heff) < 1e-6 && abs(M) < 1e-6
        dM = JA.Ms / (3 * JA.a) * (1 + JA.c);  % linear approximation
        return;
    end

    Man   = func_Man(JA, Heff);
    dMandH = func_dMan(JA, Heff);

    delta    = sign(dt);
    delta_m  = double((Man - M) * dt > 0);
    phi      = (Man - M) * delta_m;

    numerator   = JA.k * delta * JA.c * dMandH + phi;
    denominator = JA.k * delta - phi * JA.alpha - JA.k * delta * JA.alpha * JA.c * dMandH;

    if abs(denominator) < 1e-12
        dM = 0; 
    else
        dM = numerator / denominator + dMandH * JA.c;
    end
end


function Man = func_Man(JA, H)
    x = H / JA.a;
    if abs(x) < 0.1
        Man = JA.Ms * x / 3; % Taylor-approximation
    else
        Man = JA.Ms * (coth(x) - 1/x);
    end
end

function dMandH = func_dMan(JA, H)
    x = H / JA.a;
    if abs(x) < 0.1
        dMandH = JA.Ms / (3 * JA.a); % Taylor-approximation
    else
        dMandH = JA.Ms * (1 - coth(x)^2 + (1/x)^2) / JA.a;
    end
end
