
%  rainflow_stress_damage_nasgro.m  ver 1.0  by Tom Irvine

%  rainflow_stress_damage_nasgro.m  ver 1.1  by Tom Irvine
function[D]=rainflow_stress_damage_nasgro(amp_mean_cycles,A,B,C,P)
sz=size(amp_mean_cycles);
damage=0;

% DIAGNOSTIC COUNTERS
n_total=0;
n_skipped_R=0;
n_skipped_Seq=0;
n_skipped_compressive=0;
n_counted=0;
Sa_skipped_sum=0;
Sa_counted_sum=0;
Seq_counted_sum=0;

for i=1:sz(1)
    stress=amp_mean_cycles(i,1);
    S_mean=amp_mean_cycles(i,2);
    n=amp_mean_cycles(i,3);
    Smax=stress+S_mean;
    Smin=S_mean-stress;
    R=Smin/Smax;
    
    n_total=n_total+n;
    
    if(Smax<=0)
        n_skipped_compressive=n_skipped_compressive+n;
        Sa_skipped_sum=Sa_skipped_sum+(stress*n);
        continue;
    end
    
    if(R<-1)
        R=-1;
    end
    
    if(R>1)
        n_skipped_R=n_skipped_R+n;
        continue;
    end
    
    Seq=Smax*(1-R)^P;
    
    if(Seq<=C)
        n_skipped_Seq=n_skipped_Seq+n;
        continue;
    end
    
    log_Nf = A - B*log10(Seq-C);
    Nf=10^(log_Nf);
    damage=damage+n/Nf;
    n_counted=n_counted+n;
    Sa_counted_sum=Sa_counted_sum+(stress*n);
    Seq_counted_sum=Seq_counted_sum+(Seq*n);
end

% DIAGNOSTIC SUMMARY
fprintf('--- NASGRO Damage Diagnostics ---\n');
fprintf('Total cycles:                %g\n', n_total);
fprintf('Counted cycles:              %g\n', n_counted);
fprintf('Skipped (compressive):       %g\n', n_skipped_compressive);
fprintf('Skipped (R>1):               %g\n', n_skipped_R);
fprintf('Skipped (Seq<=C):            %g\n', n_skipped_Seq);
fprintf('Raw damage:                  %g\n', damage);

if(n_skipped_compressive>0)
    fprintf('Mean amplitude (skipped):    %g ksi\n', Sa_skipped_sum/n_skipped_compressive);
end
if(n_counted>0)
    fprintf('Mean amplitude (counted):    %g ksi\n', Sa_counted_sum/n_counted);
    fprintf('Mean Seq (counted):          %g ksi\n', Seq_counted_sum/n_counted);
end

fprintf('Skipped fraction:            %g %%\n', 100*n_skipped_compressive/n_total);

D=damage;