%{ 
Xue, G., Bai, H., Li, T., Ren, Z., Liu, X., & Lu, C. (2022). 
Numerical Solving Method for Jiles-Atherton Model and Influence Analysis of the Initial Magnetic Field on Hysteresis. 
Mathematics, 10(23), 4431. https://doi.org/10.3390/math10234431
%}

%%% edited: 21.11.2025

function [H, M, B, Man] = func_Jiles86(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));
    dMan = zeros(size(H));
    dMan(1,1) = 1/(3*JA.a/JA.Ms - JA.alpha);
    
    % Startvalue 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, H);
    dMandHe = func_dMan(JA, H,Man);

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

    phi = delta_m*(1-JA.c)*(Man-M);
    numerator = delta*JA.k*JA.c*dMandHe + phi;
    denominator = delta*JA.k - JA.alpha*phi;
    
    
    if abs(denominator) < 1e-12
        dM = 0;  
    else
        dM = (delta*JA.k*JA.c*dMandHe + phi)/...
        (delta*JA.k - JA.alpha*phi);
    end
end

function [dMan] = func_dMan(JA,H,Man)
    if H == 0
            dMan = 1/(3*JA.a/JA.Ms - JA.alpha);
    else

        Hf = (H + JA.alpha*Man)/JA.a;

        num = sinh(Hf)^2 -Hf^2;
        denum = JA.a*Hf^2*sinh(Hf)^2/JA.Ms - JA.alpha*(sinh(Hf)^2-Hf^2);

        dMan = num/denum;
    end
end

function [Man] = func_Man(JA,H)
    if H == 0
        Man = 0;
    else
        error = 1;
        Man_old_1 = JA.Ms*H*0.15/JA.a;
        Man_old_2 = JA.Ms*H*0.21/JA.a;
        while error > 1e-6
            num = func_f_Man(JA,H,Man_old_1)*(Man_old_1 - Man_old_2);
            denum = func_f_Man(JA,H,Man_old_1)-func_f_Man(JA,H,Man_old_2);
            Man_new = Man_old_1 - num/denum;

            Man_old_2 = Man_old_1;
            Man_old_1 = Man_new;

            error = abs(Man_old_1-Man_old_2)/abs(Man_old_1);
        end
        Man = Man_new;
    end
end
function [Man] = func_Man_init(JA,H)
    error = 1;
    Man_old_1 = 0.25;
    Man_old_2 = 0.0635;
    while error > 1e-6
        num = func_f_Man(JA,H,Man_old_1)*(Man_old_1 - Man_old_2);
        denum = func_f_Man(JA,H,Man_old_1)-func_f_Man(JA,H,Man_old_2);
        Man_new = Man_old_1 - num/denum;
        
        Man_old_2 = Man_old_1;
        Man_old_1 = Man_new;
        
        error = abs(Man_old_1-Man_old_2)/abs(Man_old_1);
    end
    Man = Man_new;
end

function f_Man = func_f_Man(JA,H,Man)
    He = (H + JA.alpha*Man)/JA.a;
    f_Man = Man - JA.Ms*( coth(He) - 1/He );
end 