function sigmaData  = spindiffusion(pdb,pdbNumCode,pdbreduced, sigmaData, conditions, user)

%%%%%%%%%%%%%%%%%%%% Display the settings used %%%%%%%%%%%%%%%%%%%%%%%%%%%%

diary('on')

fprintf(['\n\nCalculating spin diffusion with parameters:\n\n',...
        'spinDiffApp used:\t\t',        user.spinDiffApp,           '\n\n',...
       	'B0field set to:\t\t\t',        num2str(user.B0field),      '\n',...
        'noesyType set to:\t\t',        user.noesyType,             '\n',...
        'tauC set to:\t\t\t',           num2str(user.tauC),         '\n',...
       	'Amid deuteration set to:\t',   num2str(user.perCentD2O),	'\n',...
      	'pseudo methyl set to:\t\t',    num2str(user.pseudoMethyl), '\n',...
     	'upperlimit set to:\t\t',       num2str(user.upperlimit),   '\n\n',...
       	'isoAniso set to:\t\t',         num2str(user.isoAniso),     '\n\n'])
    
diary('off')

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

% Full relaxation matrix spin diffusion calculation
%{ 
 
inputs:
+ pdb with all atoms and pseudo atoms
+ pdbreduced with only the real protons and nothing else
    The two pdb should have the same SerNo for the same atoms
+ sigmaData structure for:
    the mixing time 
    the ATOMi,j_pdbSerNo
    the numcode of the atoms to check the multiplicity
    the experimental build up curve
    the two experimetal decay rates
    the Io for the initial magnetization
+ conditions structure to check the multiplicity of the pseudo atoms
+ i the index that loop over all cross peaks in the sigmaData  
+ user structure to read experimental setting like field, correlation time
    max_distance define the farthest proton involved in spin diffusion

Outputs are
            FM(j).RES1                      = []; 
            FM(j).ATOM1                     = [];
            FM(j).RES2                      = [];
            FM(j).ATOM2                     = [];

            FM(j).rPdb                      = [];  
            FM(j).sigmaDir                  = []; % 2 spins ?
            FM(j).corrFactor                = [];
            FM(j).sigma_spinDiff_corr       = [];
            FM(j).numOfPartners             = [];

            FM(j).PROCrPdb                  = []; 
            FM(j).PROCsigmaDir              = []; 
            FM(j).PROCsigma_spinDiff_corr   = [];
            FM(j).PROCcorrFactor            = []; 
            FM(j).PROCnumOfPartners         = [];     
%}


options         = optimset('display', 'off');
warning on
w               = user.B0field *1e6;
tauM            = user.tauC;
max_distance    = user.upperlimit;
size_pdbreduced = length(pdbreduced.Model.Atom);

for i = 1 : length(sigmaData)
  	sigmaData(i).FM.RES1                    = []; 
   	sigmaData(i).FM.ATOM1                   = [];
 	sigmaData(i).FM.RES2                    = [];
  	sigmaData(i).FM.ATOM2                   = [];
  	sigmaData(i).FM.rPdb                    = [];
  	sigmaData(i).FM.rho_H1                  = [];
  	sigmaData(i).FM.rho_H2                  = [];
 	sigmaData(i).FM.sigmaDir                = [];
   	sigmaData(i).FM.corrFactor              = [];
  	sigmaData(i).FM.numOfPartners           = [];
  	sigmaData(i).FM.sigmaCorr               = [];
   	sigmaData(i).FM.corrFactors             = [];
 	sigmaData(i).FM.build_up_H1H2           = []; 
 	sigmaData(i).FM.two_spins_build_up_H1H2 = [];
end

deleteScreenMsg                                  = [];


for i = 1 : length(sigmaData)
    
    if (mod(i,2) == 0) || i == length(sigmaData) || i == 1
        screenMsg = sprintf(['\n', num2str(length(sigmaData) + 1 - i), ' NOEs remaining.\n']);
        fprintf([deleteScreenMsg, screenMsg])
        deleteScreenMsg = repmat(sprintf('\b'), 1, length(screenMsg));
    
        if i == length(sigmaData)
            fprintf(deleteScreenMsg)
            fprintf('\nDONE\n')
        end
    end
    
% % %     i
% % %     sigmaData(i)
    
    %loop1 and loop2 are indices for the stereopartners if any
    
    loop1=1;
    loop2=1;
    
    loop1_end=1;
    loop2_end=1;
    
    flag_stereo1 = sigmaData(i).numcode.stereo(1);
    flag_stereo2 = sigmaData(i).numcode.stereo(2);
        
    if flag_stereo1 == 1
        loop1_end = 2;
    end
    
    if flag_stereo2 == 1
        loop2_end = 2;
    end
    
    numcode_Atom1 = sigmaData(i).numcode.ATOMi;
    %Methyl stereo partners: val(101/102)/leu(104/105)
    if flag_stereo1 == 1
        H_numCode = [];
        H_serNo   = [];
        if     numcode_Atom1 == 101
            H_numCode = [101 102];
      	elseif numcode_Atom1 == 102
            H_numCode = [102 101];
       	elseif numcode_Atom1 == 104
            H_numCode = [104 105];
        elseif numcode_Atom1 == 105
            H_numCode = [105 104];
        end
        
        ASSi = sigmaData(i).numcode.ASSi;
        for s = 1:length(pdbNumCode.Model.Atom)
        	if ~isempty(H_numCode)
            	for b = 1:length(H_numCode)
                  	if ASSi == pdbNumCode.Model.Atom(s).resSeq && H_numCode(b) == pdbNumCode.Model.Atom(s).NumCode
                        H_serNo(b)    = pdbNumCode.Model.Atom(s).AtomSerNo;
                   	end
                end
          	end
        end 
        if ~isempty(H_numCode)
            sigmaData(i).numcode.ATOMi_H            = H_numCode;
            sigmaData(i).numcode.ATOMi_HpdbSerNo    = H_serNo;
        end
    end
    
    
    numcode_Atom2 = sigmaData(i).numcode.ATOMj;
    %Methyl stereo partners: val(101/102)/leu(104/105)
    if flag_stereo2 == 1
        H_numCode = [];
        H_serNo   = [];
        if     numcode_Atom2 == 101
            H_numCode = [101 102];
      	elseif numcode_Atom2 == 102
            H_numCode = [102 101];
       	elseif numcode_Atom2 == 104
            H_numCode = [104 105];
        elseif numcode_Atom2 == 105
            H_numCode = [105 104];
        end
        
        ASSj = sigmaData(i).numcode.ASSj;
        for s = 1:length(pdbNumCode.Model.Atom)
        	if ~isempty(H_numCode)
            	for b = 1:length(H_numCode)
                  	if ASSj == pdbNumCode.Model.Atom(s).resSeq && H_numCode(b) == pdbNumCode.Model.Atom(s).NumCode
                        H_serNo(b)    = pdbNumCode.Model.Atom(s).AtomSerNo;
                   	end
                end
          	end
        end 
        if ~isempty(H_numCode)
            sigmaData(i).numcode.ATOMj_H            = H_numCode;
            sigmaData(i).numcode.ATOMj_HpdbSerNo    = H_serNo;
        end
    end
    
    for loop1 = 1 : loop1_end
        for loop2 = 1 : loop2_end

    
            time = sigmaData(i).buildups.mixTimes;

            build_up_H1     = [];
            build_up_H2     = [];
            build_up_H1H2   = [];

%Serial number for H1 and H2 (Q1 and Q2)
% and Vector coordinates [X Y Z] of H1 and H2 (Q1 and Q2)

            if flag_stereo1 == 1
                ser_num1 = sigmaData(i).numcode.ATOMi_HpdbSerNo(loop1);
            else
                ser_num1 = sigmaData(i).numcode.ATOMi_pdbSerNo;
            end
            if flag_stereo2 == 1
                ser_num2 = sigmaData(i).numcode.ATOMj_HpdbSerNo(loop2);
            else
                ser_num2 = sigmaData(i).numcode.ATOMj_pdbSerNo;
            end
        
            if ser_num1 ~= ser_num2 % check for NOE between stereo partners
                
                ser_numH = 0;

                coor_H1 = AtomCoor(pdb, ser_num1);
                coor_H2 = AtomCoor(pdb, ser_num2);
                coor_H  = [];

                % ==> Look if the protons involved in the NOE are pseudo atoms
                Q_H_cnt_1 = 1;
                Q_H_cnt_2 = 1;
                H_not_used_1 = ser_num1;
                H_not_used_2 = ser_num2;
                
                % valine and leu pseudo methyl; get atomcode from ATOMi/j_HpdbSerNo
% %                 numcode_Atom1 = sigmaData(i).numcode.ATOMi;
                if ismember( numcode_Atom1 , conditions.all.Q3H ) && flag_stereo1 == 1
                    numcode_Atom1 = sigmaData(i).numcode.ATOMi_H(loop1);
                end
                
% %                 numcode_Atom2 = sigmaData(i).numcode.ATOMj;
                if ismember( numcode_Atom2 , conditions.all.Q3H ) && flag_stereo2 == 1
                    numcode_Atom2 = sigmaData(i).numcode.ATOMj_H(loop2);
                end

                if     ismember( numcode_Atom1 , conditions.all.Q2H )
                    Q_H_cnt_1 = 2;
                    H_not_used_1 = sigmaData(i).numcode.ATOMi_HpdbSerNo;
                elseif ismember( numcode_Atom1 , conditions.all.Q3H )
                    Q_H_cnt_1 = 3;
                    H_not_used_1 = sigmaData(i).numcode.ATOMi_HpdbSerNo;
                    if flag_stereo1 == 1
                        if     numcode_Atom1 == 100
                        H_not_used_1nc = [305 315 325];
                        elseif numcode_Atom1 == 103
                        H_not_used_1nc = [335 345 355];    
                        elseif numcode_Atom1 == 106
                        H_not_used_1nc = [365 375 385];
                        elseif numcode_Atom1 == 107
                        H_not_used_1nc = [250 255 260]; 
                        elseif numcode_Atom1 == 101
                        H_not_used_1nc = [120 130 140];
                        elseif numcode_Atom1 == 102
                        H_not_used_1nc = [170 180 190];
                        elseif numcode_Atom1 == 104
                        H_not_used_1nc = [220 230 240];
                        elseif numcode_Atom1 == 105
                        H_not_used_1nc = [270 280 290];
                        end
                        
                        ASSi = sigmaData(i).numcode.ASSi;
                        for s = 1:length(pdbreduced.Model.Atom)
                            if ~isempty(H_not_used_1nc)
                                for b = 1:length(H_not_used_1nc)
                                    if ASSi == pdbreduced.Model.Atom(s).resSeq && H_not_used_1nc(b) == pdbreduced.Model.Atom(s).NumCode
                                    H_not_used_1(b)    = pdbreduced.Model.Atom(s).AtomSerNo;
                                    end
                            	end
                            end
                        end
                    end
                elseif ismember( numcode_Atom1 , conditions.all.Q4H )
                    Q_H_cnt_1 = 4;
                    H_not_used_1 = sigmaData(i).numcode.ATOMi_HpdbSerNo;
                elseif ismember( numcode_Atom1 , conditions.all.Q6H )
                    Q_H_cnt_1= 6;
                    H_not_used_1 = sigmaData(i).numcode.ATOMi_HpdbSerNo;
                end


                if     ismember( numcode_Atom2 , conditions.all.Q2H )
                    Q_H_cnt_2 = 2;
                    H_not_used_2 = sigmaData(i).numcode.ATOMj_HpdbSerNo;
                elseif ismember( numcode_Atom2 , conditions.all.Q3H )
                    Q_H_cnt_2 = 3;
                    H_not_used_2 = sigmaData(i).numcode.ATOMj_HpdbSerNo;
                    if     numcode_Atom2 == 100
                        H_not_used_2_nc = [305 315 325];
                    elseif numcode_Atom2 == 103
                        H_not_used_2_nc = [335 345 355];    
                    elseif numcode_Atom2 == 106
                        H_not_used_2_nc = [365 375 385];
                    elseif numcode_Atom2 == 107
                        H_not_used_2_nc = [250 255 260];
                    elseif numcode_Atom2 == 101
                        H_not_used_2_nc = [120 130 140];
                    elseif numcode_Atom2 == 102
                        H_not_used_2_nc = [170 180 190];
                    elseif numcode_Atom2 == 104
                        H_not_used_2_nc = [220 230 240];
                    elseif numcode_Atom2 == 105
                        H_not_used_2_nc = [270 280 290];
                    end
                    
                    ASSj = sigmaData(i).numcode.ASSj;
                    for s = 1:length(pdbreduced.Model.Atom)
                       	if ~isempty(H_not_used_2_nc)
                          	for b = 1:length(H_not_used_2_nc)
                            	if ASSj == pdbreduced.Model.Atom(s).resSeq && H_not_used_2_nc(b) == pdbreduced.Model.Atom(s).NumCode
                             	H_not_used_2(b)    = pdbreduced.Model.Atom(s).AtomSerNo;
                                end
                          	end
                      	end
                  	end   
                    
                    
                elseif ismember( numcode_Atom2 , conditions.all.Q4H )
                    Q_H_cnt_2 = 4;
                    H_not_used_2 = sigmaData(i).numcode.ATOMj_HpdbSerNo;
                elseif ismember( numcode_Atom2 , conditions.all.Q6H )
                    Q_H_cnt_2 = 6;
                    H_not_used_2 = sigmaData(i).numcode.ATOMj_HpdbSerNo;
                end

% List of all protons we do not want as partner for spin diffusion
%   all H belonging to the pseudo  
% OR
%   the two protons involved in the NOE
                H_not_used_SerNo = [ H_not_used_1 H_not_used_2 ] ;

%************************************************************
%************************************************************

                J0 = J(tauM, w, 0) ;
                J1 = J(tauM, w, 1) ;
                J2 = J(tauM, w, 2) ;

%declare
                prot = [] ;

 
%************************************************************
%        Get all protons involved in spin diffusion            
%************************************************************
 
 
                for j = 1 : size_pdbreduced % go through the pdb atoms

                    % select only protons or pseudo atoms present in library
                    % for spin diffusion simulation
                    % if pdbreduced.Model.Atom(j).H_lib == 1 

                     % Remove (from the list) protons belonging to the pseudo atom involved
                     % in the NOE

                     ser_numH = pdbreduced.Model.Atom(j).AtomSerNo;
                     Q_H_cnt  = 1; % All members from the pdb are now single protons only


                          % the proton does NOT belong to one of the two pseudo atoms
                          if ismember( ser_numH , H_not_used_SerNo ) == 0

                                coor_H = [ pdbreduced.Model.Atom(j).X pdbreduced.Model.Atom(j).Y pdbreduced.Model.Atom(j).Z ]; % get the coordinates of the H

                                if norm( coor_H - coor_H1 ) <= max_distance || norm( coor_H - coor_H2 ) <= max_distance % is H close enough to H1 or H2
                                    prot = [prot ; coor_H Q_H_cnt];

                                end
                          end

                    %end
                end




                prot = [ coor_H1 Q_H_cnt_1 ; coor_H2 Q_H_cnt_2 ; prot ]; % introduce the coordinates of H1 and H2 at the begining of the array "prot"
 
  
   
%************************************************************
%************************************************************
 
%define size of the bloc matrix 
                
                size_p = size(prot,1);
                RELAX = zeros( size_p , size_p );
 
%************************************************************
%************************************************************
 
% For pseudo atoms we introduce a factor that takes into account 
% the protons that are not treated explicitely 
 
%cross relaxation rates
                for k = 1 : size_p-1
                    for j = k+1 : size_p

                        u = [ prot(k,1) prot(k,2) prot(k,3) ] ;
                        v = [ prot(j,1) prot(j,2) prot(j,3) ] ;
                        dist = norm(u-v); 

                        RELAX(k,j) = S(dist,J0,J2) * prot(k,4) * prot(j,4) ;
                        RELAX(j,k) = S(dist,J0,J2) * prot(j,4) * prot(k,4) ;

                    end
                end
 
%auto relaxation rates 
                for k = 1 : size_p
                    for j = 1 : size_p
                        if k~=j

                            u = [ prot(k,1) prot(k,2) prot(k,3) ] ;
                            v = [ prot(j,1) prot(j,2) prot(j,3) ] ;
                            dist = norm(u-v);
                            RELAX(k,k) = RELAX(k,k) + R(dist,J0,J1,J2) * prot(j,4) * prot(k,4) ;

                        end
                    end
                end
                
%***************************************************************
%              final calcul of the NOE build up
%***************************************************************

                C = [];

                nbr_mixing_time = max ( size(time) );

                for j = 1 : nbr_mixing_time

                    C = expm( -time(j)*RELAX );
                    build_up_H1   = [ build_up_H1  ;  C(1,1) ];
                    build_up_H2   = [ build_up_H2  ;  C(2,2) ];
                    build_up_H1H2 = [ build_up_H1H2;  C(2,1) ];

                end

%***************************************************************
%          Isolated two spins theoretical build up NOE
%***************************************************************

                % Save the sigma_H1H2 for the 2 spins system
                sigma_H1H2 = RELAX(2,1);

                rho_H1 = RELAX(1,1);
                rho_H2 = RELAX(2,2);
                
                % For Q atoms we modified the rate and kept Io = 1
                sigmaEff_scaling= [sigmaData(i).numcode.sigmaEffScaling(1) 1 sigmaData(i).numcode.sigmaEffScaling(3)] ; %
                
                lambda_p        = lambdaPlusScaledFun(rho_H1 ,rho_H2,sigmaEff_scaling,sigma_H1H2);
                lambda_m        = lambdaMinusScaledFun(rho_H1 ,rho_H2,sigmaEff_scaling,sigma_H1H2);

                % Estimated two_spins_build_up_H1H2 
                two_spins_build_up_H1H2 = two_spins_build_up_Fun(lambda_p, lambda_m, sigma_H1H2, time, sigmaEff_scaling);
                                


%***************************************************************
%        Spin diffusion correction factors for the build_up
%***************************************************************

                exp_build_up_corrections = two_spins_build_up_H1H2./build_up_H1H2';

%***************************************************************
%            Fit sigma on the corrected NOE build up
%***************************************************************


                % Get experimental data
                exp_build_up_H1H2   = sigmaData(i).buildups.exp_build_up_H1H2_norm;
                
                exp_rho_H1          = sigmaData(i).rhoi;
                exp_rho_H2          = sigmaData(i).rhoj;


                % correct the NOE build up
                y                   = exp_build_up_H1H2 .* exp_build_up_corrections;




% Correcting factor due to degenerated group of protons "Q"
% If the noe is between two atoms groups with different number of protons
% then we need to correct the delta_I_zero by a factor that is the number 
% of protons, "Q_H_cnt", represented by the pseudo atom "Q" of 
% the non-normalizing peak
 

%if the 2nd protons group is a pseudo atom "Q" we correct by the implicit
%number of protons
       

%                 exp_delta_I_zero_H1 = abs( sigmaData(i).Izeroij )* Q_H_cnt_2;
%                 normalized_y = y / exp_delta_I_zero_H1;
                normalized_y = y;

                
 % FIT the corrected build up
 
                start = -( exp_rho_H1 + exp_rho_H2 )/2;
% %                 sigmaEff_scaling = [sigmaData(i).numcode.sigmaEffScaling(1) 1 sigmaData(i).numcode.sigmaEffScaling(3)] ; %
                sigmaEff_scaling = sigmaData(i).numcode.sigmaEffScaling;
                
                 [corrected_estimated_exp_sigma, sse_target_fun_sigma, exitflag_sigma] = fminsearch(@(x)fitsigma(x, exp_rho_H1, exp_rho_H2, time ,normalized_y, sigmaEff_scaling), start, options);





%***************************************************************
%                   Compiling output
%***************************************************************

                sigmaData(i).FM.RES1                 = sigmaData(i).ASSi; 
                sigmaData(i).FM.ATOM1                = sigmaData(i).numcode.ATOMi;
                sigmaData(i).FM.RES2                 = sigmaData(i).ASSj;
                sigmaData(i).FM.ATOM2                = sigmaData(i).numcode.ATOMj;
                
                sigmaData(i).FM.rPdb                 = [ sigmaData(i).FM.rPdb  norm( coor_H1 - coor_H2 ) ]; 
                
                sigmaData(i).FM.rho_H1               = [ sigmaData(i).FM.rho_H1 rho_H1];
                sigmaData(i).FM.rho_H2               = [ sigmaData(i).FM.rho_H2 rho_H2];
                sigmaData(i).FM.sigmaDir             = [ sigmaData(i).FM.sigmaDir sigma_H1H2 ]; % 2 spins ?

                sigmaData(i).FM.corrFactor           = [ sigmaData(i).FM.corrFactor abs( corrected_estimated_exp_sigma / sigmaData(i).sigmaij ) ];
                sigmaData(i).FM.sigmaCorr            = [ sigmaData(i).FM.sigmaCorr corrected_estimated_exp_sigma ]; %
                sigmaData(i).FM.numOfPartners        = [ sigmaData(i).FM.numOfPartners   max( size(RELAX) ) - 2 ]; % 2 for the two atoms involved in the NOE
            	
                sigmaData(i).FM.corrFactors          = [ sigmaData(i).FM.corrFactors exp_build_up_corrections ]; 
                sigmaData(i).FM.build_up_H1H2        = build_up_H1H2; 
                sigmaData(i).FM.two_spins_build_up_H1H2 = two_spins_build_up_H1H2;
            end
            
        end
    end
  
    sigmaData(i).FM.PROCrPdb                        = mean( sigmaData(i).FM.rPdb );  
    sigmaData(i).FM.PROCsigmaDir                    = mean( sigmaData(i).FM.sigmaDir ); 
    sigmaData(i).FM.PROC_sigma_spinDiff_corr        = mean( sigmaData(i).FM.sigmaCorr );
    sigmaData(i).FM.PROCcorrFactor                  = mean( sigmaData(i).FM.corrFactor );
            
            
% Double assignment for comparison with 3 spins.            
            
    sigmaData(i).spinDiff.RES1                      = sigmaData(i).FM.RES1; 
    sigmaData(i).spinDiff.ATOM1                     = sigmaData(i).FM.ATOM1;
    sigmaData(i).spinDiff.RES2                      = sigmaData(i).FM.RES2;
    sigmaData(i).spinDiff.ATOM2                     = sigmaData(i).FM.ATOM2;
    
    sigmaData(i).spinDiff.sigmaDir                  = sigmaData(i).FM.sigmaDir;       
    sigmaData(i).spinDiff.corrFactor                = sigmaData(i).FM.corrFactor;               
    sigmaData(i).spinDiff.sigma_spinDiff_corr       = sigmaData(i).FM.sigmaCorr; 
    sigmaData(i).spinDiff.numOfPartners             = sigmaData(i).FM.numOfPartners; 
      
    sigmaData(i).spinDiff.PROCrPdb                  = sigmaData(i).FM.PROCrPdb;        
    sigmaData(i).spinDiff.PROCsigmaDir              = sigmaData(i).FM.PROCsigmaDir; 
    sigmaData(i).spinDiff.PROC_sigma_spinDiff_corr  = sigmaData(i).FM.PROC_sigma_spinDiff_corr;
    sigmaData(i).spinDiff.PROCcorrFactor            = sigmaData(i).FM.PROCcorrFactor; 
    
end


 



%{ 
@nd method for sigma correction
    
 
    start = 1;
    % Fit rho1
    y = build_up_H1;
    [estimated_rho1, sse_target_fun_rho1, exitflag_rho1] = fminsearch(@(x)fitrho(x,t,y),start);

    % Fit rho2
    y = build_up_H2;
    [estimated_rho2, sse_target_fun_rho2, exitflag_rho2] = fminsearch(@(x)fitrho(x,t,y),start);

    % Fit sigma
    y = build_up_H1H2;
    delta_I_zero_H1 = 1;
    rho_H1 = estimated_rho1;
    rho_H2 = estimated_rho2;
    [estimated_sigma_spin_diffusion, sse_target_fun_sigma, exitflag_sigma] = fminsearch(@(x)fitsigma(x, delta_I_zero_H1, rho_H1, rho_H2, t ,y), start);

    
    
    % Make correction on the experimental fitted sigma
    correction = sigma_H1H2/estimated_sigma_spin_diffusion;
    corrected_exp_sigma = exp_sigma * correction ;

    
    % Convert the corrected_sigma to a corrected distance
     sigma2distance(corrected_sigma, J2, J0);
    

        
    %}
