function diagDecay = diagDecayFit8Stats(diagDecay)

fprintf('\nExponential fitting of rho to diagonal peak intensities\n')

% Fit rho to diagonal peak intensities
for i = 1:length(diagDecay)
    
    if ~isempty(diagDecay(i).mixTimes)

    
        time                = diagDecay(i).mixTimes;
        exp_decay_H1        = diagDecay(i).DPrho;
    
        [estimates, model] 	= fitfunRhoSingle(time, exp_decay_H1);
        exp_rho_H1          = estimates(1);
        exp_const_H1        = estimates(2);
    
    
        diagDecay(i).rho   	= exp_rho_H1;
        diagDecay(i).Const 	= exp_const_H1;
    
    end
    
end

% Statistics for diagonal peaks
for i = 1:length(diagDecay)
    
    if ~isempty(diagDecay(i).mixTimes)
    
        for k = 1:length(diagDecay(i).mixTimes)
        
            diagDecay(i).PointtheoRho(k) 	= diagDecay(i).Const * exp(-diagDecay(i).rho * diagDecay(i).mixTimes(k));
            diagDecay(i).diffRho(k)     	= (diagDecay(i).PointtheoRho(k) - diagDecay(i).DPrho(k));
        
        end
    
        diagDecay(i).chisqRho               = sum(diagDecay(i).diffRho.^ 2) / (length(diagDecay(i).mixTimes) - 0);
        diagDecay(i).chisq05Rho         	= diagDecay(i).chisqRho^0.5;
        diagDecay(i).chiNRho            	= diagDecay(i).chisq05Rho/max(abs(diagDecay(i).DPrho));
        
    end
end