%
%   vibfatigue_modified_dirlik.m  ver 1.1  by Tom Irvine
%
%   email: tom@irvinemail.org
%
function result = vibfatigue_modified_dirlik(psd_input, C, k, T, opts)
%% VIBFATIGUE_MODIFIED_DIRLIK - Spectral Fatigue (E_band Smooth Correction)
%
% ASCII-only version (FTP-safe, no encoding issues)
%
% Core model:
%   log_C = (1 - E_band) * [ b0*(1-alpha2^2) + b1*alpha2*(1-alpha2) ]
%

%% 0. Defaults
if nargin < 5, opts = struct; end
if ~isfield(opts,'useCorrection'), opts.useCorrection = true; end
if ~isfield(opts,'kInterp'),       opts.kInterp = true; end
if ~isfield(opts,'verbose'),       opts.verbose = true; end
if ~isfield(opts,'stressUnit'),    opts.stressUnit = 'units'; end

%% 1. Input checks
f   = psd_input(:,1);
Gxx = psd_input(:,2);

assert(all(diff(f)>0),'Frequency must be increasing');
assert(all(Gxx>=0),'PSD must be non-negative');

%% 2. Calibration table (LATEST)
k_cal = [3 4 5 6 7 8 9 10 11];
B_cal = [
  0.517938, -0.921463;  % 3
  0.839129, -1.606710;  % 4
  1.206716, -2.456515;  % 5
  1.573849, -3.338758;  % 6
  1.914259, -4.177807;  % 7 
  2.217859, -4.942144;  % 8
  2.482905, -5.623617;  % 9
  2.679377,   -6.165466;   % 10 
  2.864351,   -6.671826];  % 11

%% 3. Spectral moments
m0 = trapz(f,Gxx);
m1 = trapz(f,f.*Gxx);
m2 = trapz(f,f.^2.*Gxx);
m4 = trapz(f,f.^4.*Gxx);

assert(m0>0 && m2>0 && m4>0,'Invalid PSD moments');

%% 4. Derived parameters
alpha2 = m2/sqrt(m0*m4);
alpha2 = max(0,min(1,alpha2));

nu0_calc = sqrt(m2/m0);
nuP_calc = sqrt(m4/m2);

nu0_Hz = nu0_calc/(2*pi);
nuP_Hz = nuP_calc/(2*pi);

sigma = sqrt(m0);

%% 5. E_band
[~,ipk] = max(Gxx);
f_pk = f(ipk);

f_lo = f_pk/sqrt(2);
f_hi = f_pk*sqrt(2);

mask = (f>=f_lo & f<=f_hi);

if any(mask)
    E_band = trapz(f(mask),Gxx(mask))/m0;
else
    E_band = 0;
end

E_band = max(0,min(1,E_band));

%% 6. Dirlik parameters (robust)
xm = m1/m0*sqrt(m2/m4);

D1 = max(2*(xm-alpha2^2)/(1+alpha2^2),1e-6);

den = 1-alpha2-D1+D1^2;

if abs(den)<1e-12
    R = 0;
else
    R = (alpha2-xm-D1^2)/den;
end

if abs(1-R)<1e-12
    D2 = 0;
else
    D2 = max(min(den/(1-R),1),0);
end

D3 = max(1-D1-D2,1e-6);

Q = max(1.25*(alpha2-D3-D2*R)/D1,1e-6);

%% 7. Raw Dirlik damage
kk = k;
mh = kk/2;

T1 = D1*Q^kk*gamma(kk+1);
T2 = (sqrt(2)^kk)*gamma(mh+1)*(D2*abs(R)^kk + D3);

D_rate_raw = (nuP_calc/C)*(m0^mh)*(T1+T2);

%% 8. Correction
correction = 1.0;
logC = 0;

if opts.useCorrection

    % interpolate coefficients
    if opts.kInterp
        b0 = interp1(k_cal,B_cal(:,1),k,'pchip','extrap');
        b1 = interp1(k_cal,B_cal(:,2),k,'pchip','extrap');
    else
        [~,ik]=min(abs(k_cal-k));
        b0=B_cal(ik,1);
        b1=B_cal(ik,2);
    end

    phi0 = 1 - alpha2^2;
    phi1 = alpha2*(1-alpha2);

    scale = 1 - E_band;

    logC = scale*(b0*phi0 + b1*phi1);

    correction = exp(logC);
end

%% 9. Final outputs
D_rate = D_rate_raw * correction;
D_total = D_rate*T;
life = 1/max(D_rate,realmin);

%% 10. Output struct
result.D_rate = D_rate;
result.D_total = D_total;
result.D_rate_raw = D_rate_raw;
result.life_hours = life/3600;

result.alpha2 = alpha2;
result.E_band = E_band;
result.correction_factor = correction;
result.log_C_correction = logC;

result.nu0 = nu0_Hz;
result.nuP = nuP_Hz;
result.sigma_rms = sigma;

result.moments = [m0 m1 m2 m4];

result.dirlik_params = struct('D1',D1,'D2',D2,'D3',D3,'Q',Q,'R',R);

%% 11. Console output
if opts.verbose
    fprintf('\n--- Modified Dirlik ---\n');
    fprintf('alpha2 = %.4f\n',alpha2);
    fprintf('E_band = %.4f\n',E_band);
    fprintf('Correction = %.4f\n',correction);
    fprintf('Raw D_rate = %.4e\n',D_rate_raw);
    fprintf('Corr D_rate = %.4e\n',D_rate);
end

end