function [diagDecay] = readDiagDecaySup(diagDecay,pdbHQnumCode,user)

if ~strcmp(user.buildupRhoFileName, 'n') 
    fprintf('\nUsing user supplied rho and izero values\n')

    
    file                                            = fopen (user.buildupRhoFileName, 'r');
    dataRhoInt                                      = fscanf (file, '%f  %f  %f  %f  %f', [5 inf])';
    fclose (file);

    for i = 1:length(dataRhoInt)
        
        diagDecaySup(i).numcode.ASSi                = dataRhoInt(i,1);
    	diagDecaySup(i).numcode.ATOMi               = dataRhoInt(i,2);
        diagDecaySup(i).fullSetDPrho                = zeros(length(user.mixing_times));
        diagDecaySup(i).rho                         = dataRhoInt(i,3);
        diagDecaySup(i).Const                       = dataRhoInt(i,4);
        
        ASS                                         = diagDecaySup(i).numcode.ASSi;
    	ATOM                                        = diagDecaySup(i).numcode.ATOMi;
      	ATOM_pdbSerNo                               = pdbSerNoFun(ASS,ATOM,pdbHQnumCode);
    
      	diagDecaySup(i).numcode.ATOMi_pdbSerNo      = ATOM_pdbSerNo;
    
        if isempty(ATOM_pdbSerNo)
            fprintf(['\nNOT all diagonal peak assignments have been matched to the pdb file'...
             '\nMissing: ' num2str(ASS) 'and' num2str(ATOM) 'in your ' user.buildupRhoFileName 'file'...
             '\nThe spindiffusion calculation is likely to fail!\n'])
        end
        
    end
    


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
% Overwrite fitted data with supplied if available
   
    for i=1:length(diagDecay)
       
        for k = 1:length(diagDecaySup)
            
            if diagDecaySup(k).numcode.ASSi     == diagDecay(i).numcode.ASSi...
            && diagDecaySup(k).numcode.ATOMi   	== diagDecay(i).numcode.ATOMi
                
                diagDecay(i).numcode            = diagDecaySup(k).numcode;
                diagDecay(i).fullSetmixTimes  	= [];
                diagDecay(i).fullSetDPrho    	= diagDecaySup(k).fullSetDPrho;
                diagDecay(i).mixTimes         	= [];
                diagDecay(i).DPrho            	= [];
                diagDecay(i).rho               	= [];
                diagDecay(i).Const            	= [];
                diagDecay(i).chiNRho         	= 999;
                diagDecay(i).PointtheoRho     	= [];
                diagDecay(i).diffRho          	= [];
                diagDecay(i).chisqRho          	= [];
                diagDecay(i).chisq05Rho      	= [];
                
                diagDecay(i).rho                = diagDecaySup(k).rho;
                diagDecay(i).Const            	= diagDecaySup(k).Const;
                
            end
            
        end
        
    end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% add a diagDecay Structure entry for supplied from file.

    addCount                                = length(diagDecay) + 1;

    for k = 1:length(diagDecaySup)
    
        foundFlag                               = 0;
    
        for i=1:length(diagDecay)
        
            if diagDecaySup(k).numcode.ASSi         == diagDecay(i).numcode.ASSi...
            && diagDecaySup(k).numcode.ATOMi        == diagDecay(i).numcode.ATOMi
            
                foundFlag                           = 1;
            
            end
                
        end
    
        if foundFlag == 0
        
      	diagDecay(addCount).numcode                 = diagDecaySup(k).numcode;
% % %        	diagDecay(addCount).resRho                  = diagDecay(addCount).numcode.ASSi;
% % %        	diagDecay(addCount).atomRho                 = diagDecay(addCount).numcode.ATOMi;
       	diagDecay(addCount).secondStruct            = [];
       	diagDecay(addCount).fullSetmixTimes         = [];
      	diagDecay(addCount).fullSetDPrho            = diagDecaySup(k).fullSetDPrho;
      	diagDecay(addCount).mixTimes                = [];
      	diagDecay(addCount).DPrho                   = [];
     	diagDecay(addCount).rho                     = diagDecaySup(k).rho;
    	diagDecay(addCount).Const                   = diagDecaySup(k).Const;
     	diagDecay(addCount).chiNRho                 = 999;
      	diagDecay(addCount).PointtheoRho            = [];
      	diagDecay(addCount).diffRho                 = [];
      	diagDecay(addCount).chisqRho                = [];
      	diagDecay(addCount).chisq05Rho              = [];
        
        addCount                                    = addCount + 1;
        
        end
            
        end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%