function sigmaData  = threeSpinSim2(pdb,pdbNumCode,pdbreduced, sigmaData, diagDecay,conditions, user,work)

%%%%%%%%%%%%%%%%%%%% 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',...
        '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'])
    
diary('off')

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

% three 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
%}


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).threeSpin.RES1                     = []; 
   	sigmaData(i).threeSpin.ATOM1                    = [];
 	sigmaData(i).threeSpin.RES2                     = [];
  	sigmaData(i).threeSpin.ATOM2                    = [];
  	sigmaData(i).threeSpin.rPdb                     = [];
  	sigmaData(i).threeSpin.rho_H1                   = [];
  	sigmaData(i).threeSpin.rho_H2                   = [];
 	sigmaData(i).threeSpin.sigmaDir                 = [];
   	sigmaData(i).threeSpin.corrFactor               = [];
  	sigmaData(i).threeSpin.numOfPartners            = [];
  	sigmaData(i).threeSpin.sigmaCorr                = [];
   	sigmaData(i).threeSpin.corrFactors              = [];
 	sigmaData(i).threeSpin.build_up_H1H2            = []; 
 	sigmaData(i).threeSpin.two_spins_build_up_H1H2  = [];
end

autoRELAX           = rhoSim(pdbreduced,pdbNumCode,diagDecay, conditions,user, work);


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

    %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;            
            resSeq 	= [];
           	NumCode	= [];
           	resName	= [];
%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);
                resSeq(1)  	= sigmaData(i).numcode.ASSi;
                NumCode(1)  = sigmaData(i).numcode.ATOMi;
                resName(:,1)= sigmaData(i).ASSi_resName;
            else
                ser_num1 = sigmaData(i).numcode.ATOMi_pdbSerNo;
                resSeq(1)  	= sigmaData(i).numcode.ASSi;
                NumCode(1)  = sigmaData(i).numcode.ATOMi;
                resName(:,1)= sigmaData(i).ASSi_resName;
            end
            if flag_stereo2 == 1
                ser_num2 = sigmaData(i).numcode.ATOMj_HpdbSerNo(loop2);
                resSeq(2)  	= sigmaData(i).numcode.ASSj;
                NumCode(2)  = sigmaData(i).numcode.ATOMj;
                resName(:,2)= sigmaData(i).ASSj_resName;
            else
                ser_num2 = sigmaData(i).numcode.ATOMj_pdbSerNo;
                resSeq(2)  	= sigmaData(i).numcode.ASSj;
                NumCode(2)  = sigmaData(i).numcode.ATOMj;
                resName(:,2)= sigmaData(i).ASSj_resName;
            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
                H_not_used_1 = ser_num1;
                H_not_used_2 = ser_num2;
                
                % valine and leu pseudo methyl; get atomcode from ATOMi/j_HpdbSerNo
                if ismember( numcode_Atom1 , conditions.all.Q3H ) && flag_stereo1 == 1
                    numcode_Atom1 = sigmaData(i).numcode.ATOMi_H(loop1);
                end
                
                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 )
                    H_not_used_1 = sigmaData(i).numcode.ATOMi_HpdbSerNo;
                elseif ismember( numcode_Atom1 , conditions.all.Q3H )
                    H_not_used_1 = sigmaData(i).numcode.ATOMi_HpdbSerNo;
                    if flag_stereo1 == 1
                        if     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 )
                    H_not_used_1 = sigmaData(i).numcode.ATOMi_HpdbSerNo;
                elseif ismember( numcode_Atom1 , conditions.all.Q6H )
                    H_not_used_1 = sigmaData(i).numcode.ATOMi_HpdbSerNo;
                end


                if     ismember( numcode_Atom2 , conditions.all.Q2H )
                    H_not_used_2 = sigmaData(i).numcode.ATOMj_HpdbSerNo;
                elseif ismember( numcode_Atom2 , conditions.all.Q3H )
                    H_not_used_2 = sigmaData(i).numcode.ATOMj_HpdbSerNo;
                    if flag_stereo1 == 1
                        if     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  
                    end
                    
                elseif ismember( numcode_Atom2 , conditions.all.Q4H )
                    H_not_used_2 = sigmaData(i).numcode.ATOMj_HpdbSerNo;
                elseif ismember( numcode_Atom2 , conditions.all.Q6H )
                    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) ;

                prot = [] ;

%************************************************************
%        Get all protons involved in spin diffusion            
%************************************************************
 
                k           = 2;           
                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
                                k = k+1;
                                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];
                                    NumCode(k)	= pdbreduced.Model.Atom(j).NumCode;
                                    resName(:,k)= pdbreduced.Model.Atom(j).resName;
                                    resSeq(k)	= pdbreduced.Model.Atom(j).resSeq;
                                    ser_num(k)	= ser_numH;
                                end
                          end
                end
                Q_H_cnt_1 = 1;
                Q_H_cnt_2 = 1;
                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);
                        RELAX(j,k) = S(dist,J0,J2);

                    end
                end
 
%auto relaxation rates 
                for k = 1 : size_p
                    for j = 1 : size_p
                        if k~=j
                            for a = 1:size(autoRELAX,1)
                                if autoRELAX(a,3) == ser_num(k)
                                RELAX(k,k) 	= autoRELAX(a,4);
                                end
                            end

                        end
                    end
                end
                                
%***************************************************************
%              final calcul of the NOE build up
%***************************************************************
                
                two_spins_build_up_H1H2                 = [];
                two_spins_build_up_H1H2_scaled          = [];
                three_spins_build_up_H1H2_scaled        = [];
                s_two_spins_build_up_H1H2               = [];
                s_two_spins_build_up_H1H2_scaled_sum    = [];
                s_three_spins_build_up_H1H2_scaled_sum  = [];
                s_numPartners                           = [];
                s_H1H2                                  = [];
                
                %***************************************************************
                %          Isolated two spins theoretical build up NOE
                %***************************************************************
                rho_H1      = RELAX(1,1);
                rho_H2      = RELAX(2,2);
                sigma_H1H2  = RELAX(2,1);

% %                 sigmaEff_scaling = [1 1 1] ;
                % 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);
                
                %-----------------------------------------------
                RELAX3s      = zeros(3,3);
                RELAX3s(1,1) = RELAX(1,1);
                RELAX3s(1,2) = RELAX(1,2);
                RELAX3s(2,1) = RELAX(2,1);
                RELAX3s(2,2) = RELAX(2,2);
                
                ntime                           = max(size(time));
                for k = 3:size(RELAX,2);
                
                    resSeq_H3                   = resSeq(k);
                    NumCode_H3                  = NumCode(k);
                    resName_H3                  = resName(k);      	% convert char to double before comparing
                
                    RELAX3s(1,3)                = RELAX(1,k);
                    RELAX3s(3,1)                = RELAX(1,k);
                    RELAX3s(2,3)                = RELAX(2,k);
                    RELAX3s(3,2)                = RELAX(2,k);
                    RELAX3s(3,3)                = RELAX(k,k);
                    
                    %-----------------------------------------------
                    %expm approach
                    C                        	= [];
                    three_spins_build_up_H1H2   = [];
                    ntime                       = max(size(time));
                    for t = 1 : ntime
                        C = expm( -time(t)*RELAX3s );
                        three_spins_build_up_H1H2(t)= C(2,1);
                    end
                    
                    
% % %                     %-----------------------------------------------
% % %                     %The Solomon differential equation
% % %                     mixtimes= sigmaData(i).buildups.mixTimes;
% % %                     time    = (0:0.1:max(mixtimes) * 10^3) * 10^-3;
% % %                     
% % %                     three_spins_build_up_H1H2   = [];
% % %                     
% % %                     gamma                       = [RELAX3s(1,1) RELAX3s(1,3) RELAX3s(1,2); ...
% % %                                                    RELAX3s(1,3) RELAX3s(3,3) RELAX3s(2,3); ...
% % %                                                    RELAX3s(1,2) RELAX3s(2,3) RELAX3s(2,2)];
% % %                     options                     = [];   
% % %                     [t,buildup]                 = ode45(@Solomon,time,[1 0 0],options,gamma);
% % %                     three_spins_build_up_H1H2   = buildup(:,3);
% % % 
% % %                     three_spins_build_up_H1H2_scaled_sel  	= selBuildupAtMixtimes(three_spins_build_up_H1H2, time, mixtimes);
% % %                     three_spins_build_up_H1H2   = [];
% % %                     three_spins_build_up_H1H2   = three_spins_build_up_H1H2_scaled_sel;
% % %                     time = sigmaData(i).buildups.mixTimes;
% % %                     %-----------------------------------------------
                    
                    two_spins_build_up_H1H2_scaled(:,k)  	= scale2labSchemeFun(two_spins_build_up_H1H2,NumCode_H3,resName_H3,conditions,user);
                    three_spins_build_up_H1H2_scaled(:,k) 	= scale2labSchemeFun(three_spins_build_up_H1H2,NumCode_H3,resName_H3,conditions,user); 
                end
                s_two_spins_build_up_H1H2_scaled_sum        = sum(two_spins_build_up_H1H2_scaled,2);
                s_three_spins_build_up_H1H2_scaled_sum      = sum(three_spins_build_up_H1H2_scaled,2); 
                s_numPartners                               = size(path,2);
                s_H1H2                                      = sigma_H1H2;
                                


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

                exp_build_up_corrections   	= two_spins_build_up_H1H2' ./ (s_three_spins_build_up_H1H2_scaled_sum - s_two_spins_build_up_H1H2_scaled_sum + two_spins_build_up_H1H2');

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


% correct the NOE build up
                exp_rho_H1 = sigmaData(i).rhoi;
                exp_rho_H2 = sigmaData(i).rhoj;
                
                exp_build_up_H1H2 = sigmaData(i).buildups.exp_build_up_H1H2_norm;
                normalized_y = exp_build_up_H1H2 .* exp_build_up_corrections';

 % FIT the corrected build up
 
                start = -( exp_rho_H1 + exp_rho_H2 )/2;
                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).threeSpin.RES1                 = sigmaData(i).ASSi; 
                sigmaData(i).threeSpin.ATOM1                = sigmaData(i).numcode.ATOMi;
                sigmaData(i).threeSpin.RES2                 = sigmaData(i).ASSj;
                sigmaData(i).threeSpin.ATOM2                = sigmaData(i).numcode.ATOMj;
                
                sigmaData(i).threeSpin.rPdb                 = [ sigmaData(i).threeSpin.rPdb  norm( coor_H1 - coor_H2 ) ]; 
                
                sigmaData(i).threeSpin.rho_H1               = [ sigmaData(i).threeSpin.rho_H1 rho_H1];
                sigmaData(i).threeSpin.rho_H2               = [ sigmaData(i).threeSpin.rho_H2 rho_H2];
                sigmaData(i).threeSpin.sigmaDir             = [ sigmaData(i).threeSpin.sigmaDir sigma_H1H2 ]; 

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