manual_lau_holtzman_fig2.m

function manual_lau_holtzman_fig2()
    % recreate fig 2 of lau and holtzman 2019 grl
    % https://doi.org/10.1029/2019GL083529
    beta = 1e-4
    alf = 1/3

    Mu_in = 60*1e9;
    eta_ss =  1.8922e+21;
    tau_M = eta_ss / Mu_in;

    f = logspace(-15,1, 100);
    w = 2 * pi * f;

    % andrade
    M_real = 1 + beta * gamma(1+alf) * cos(alf * pi /2) ./ (w.^alf);
    M_imag = 1. ./ (w .* tau_M) + beta * gamma(1+alf)*sin(alf*pi/2)./(w.^alf);
    M_fac = M_real - i * M_imag;
    M_complex = Mu_in ./ M_fac;
    [and_Qinv, and_eta_app, and_eta_normalized] = get_complex_visc(M_complex, tau_M, eta_ss, w,0);

    % burgers
    eta_t = eta_ss/10;
    M_t = Mu_in * 0.8;
    M_complex = i * w * eta_ss .* (1 + i * w * eta_t/M_t);
    eta_fac = eta_ss / M_t + eta_ss/Mu_in + eta_t / M_t + i * w * eta_ss * eta_t / (Mu_in * M_t);
    M_complex = M_complex ./ (1 + i * w .* eta_fac);
    [burg_Qinv, burg_eta_app, burg_eta_normalized] = get_complex_visc(M_complex, tau_M, eta_ss, w,1);

    % maxwell
    M_complex = i * w * eta_ss ./ (1 + i * w * tau_M);
    [mxwl_Qinv, mxwl_eta_app, mxwl_eta_normalized] = get_complex_visc(M_complex, tau_M, eta_ss, w,1);

    % Zener
    M_complex = Mu_in * M_t / (Mu_in + M_t) * (1 + i*w*eta_t/M_t);
    M_complex = M_complex ./ (1 + i * w * eta_t / (Mu_in + M_t));
    [z_Qinv, z_eta_app, z_eta_normalized] = get_complex_visc(M_complex, tau_M, eta_ss, w,1);



    tau_f = 1./ tau_M;

    figure()
    subplot(3,1,1)
    f_axis = f;
    loglog(f_axis, mxwl_eta_app, 'linewidth', 1.0,'displayname', 'maxwell')
    hold on
    loglog(f_axis, and_eta_app,'r', 'linewidth', 1.0,'displayname', 'andrade')
    loglog(f_axis, burg_eta_app,'m', 'linewidth', 1.0,'displayname', 'burger')
    loglog(f_axis, z_eta_app,'--k', 'linewidth', 1.0,'displayname', 'zener')
    loglog([tau_f, tau_f], [1e12,1e24],'--k', 'displayname', 'f_M')
    ylim([1e12,1e24])
    ylabel('||\eta*||')
    xlim([1e-13,10])
    legend()

    subplot(3,1,2)
    loglog(f_axis, mxwl_Qinv, 'linewidth', 1.0,'displayname', 'maxwell')
    hold on
    loglog(f_axis, and_Qinv, 'r', 'linewidth', 1.0,'displayname', 'andrade')
    loglog(f_axis, burg_Qinv,'m', 'linewidth', 1.0,'displayname', 'burger')
    loglog(f_axis, z_Qinv,'--k', 'linewidth', 1.0,'displayname', 'zener')
    loglog([tau_f, tau_f], [1e-8,1e4],'--k')
    ylabel('Q^{-1}')
    xlim([1e-13,10])
    ylim([1e-8,1e4])

    subplot(3,1,3)
    semilogx(f_axis, mxwl_eta_normalized, 'linewidth', 1.0,'displayname', 'maxwell')
    hold on
    semilogx(f_axis, and_eta_normalized, 'r', 'linewidth', 1.0,'displayname', 'andrade')
    semilogx(f_axis, burg_eta_normalized,'m', 'linewidth', 1.0,'displayname', 'burger')
    semilogx(f_axis, z_eta_normalized, '--k', 'linewidth', 1.0,'displayname', 'zener')
    semilogx([tau_f, tau_f], [0, 1.5],'--k')
    semilogx([f_axis(1), f_axis(end)], [1,1],'--k')
    ylabel('normalized ||{\eta}*||')
    xlabel('f [Hz]')
    xlim([1e-13,10])
    ylim([0, 1.5])

end

function [Qinv, eta_app, eta_normalized] = get_complex_visc(M_complex, tau_M, eta_ss, w, use_Q_fac)
    M1 = real(M_complex);
    M2 = imag(M_complex);

    J1 = real(1. ./ M_complex);
    J2 = imag(1. ./ M_complex);
    if use_Q_fac == 1
        M1_M2_fac = (1 + sqrt(1+(M2./M1).^2)) / 2;
    else
        M1_M2_fac = 1.0;
    end
    Qinv = (M2./M1) ./ M1_M2_fac;

    % complex compliance
    J = J1 + i * J2;

    % complex modulus
    M = 1./J;

    % complex viscosity
    eta_star= -i .* (1 ./ w) .* M;  % complex viscosity

    % apparent viscosity
    eta_app = abs(eta_star);

    % complex maxwell viscosity
    eta_maxwell = ( eta_ss ./ (1 + i * w * tau_M));

    % maxwell-normalized apparent viscosity
    eta_normalized = abs(eta_star) ./ abs(eta_maxwell);

end