function simul = simulations(user,conditions,threeSpin,pathway)


% Solving the Solomon equation numerically and fitting two_spin_build_up 
mixtimes                                    = threeSpin.buildups.mixTimes;
time                                        = mixtimes;                                        
% % time = (0:0.1:max(mixtimes) * 10^3) * 10^-3;


ATOMi_resSeq                                = pathway.ATOMi_resSeq;
ATOMi_NumCode                               = pathway.ATOMi_NumCode;
ATOMi_HpdbSerNo                             = pathway.ATOMi_HpdbSerNo;

ATOMj_resSeq                                = pathway.ATOMj_resSeq;
ATOMj_NumCode                               = pathway.ATOMj_NumCode;
ATOMj_HpdbSerNo                             = pathway.ATOMj_HpdbSerNo;

rho_H1                                      = pathway.ATOMi_rho;           % rho
rho_H2                                      = pathway.ATOMj_rho;           % rho
sigma_H1H2                                  = pathway.sigma_H1H2;          % sigma12
rPdb_H1H2                                   = pathway.rPdb_H1H2;

% two_spin_build_up

% we take the ratios of the two spin buildup simulation and three spin simulation
% therefore we set the scaling to one.           

sigmaEff_scaling                            = [1 1 1]; 

lambda_p                                    = lambdaPlusScaledFun(abs(rho_H1),abs(rho_H2),sigmaEff_scaling,sigma_H1H2);
lambda_m                                    = lambdaMinusScaledFun(abs(rho_H1),abs(rho_H2),sigmaEff_scaling,sigma_H1H2);
two_spins_build_up_H1H2                     = two_spins_build_up_Fun(lambda_p,lambda_m,sigma_H1H2,time,sigmaEff_scaling);

%-----------------------------------------------
partnerCount = length(pathway.spinDiffPartner);

for j = 1:partnerCount
    
    resSeq_H3(j)                            = pathway.spinDiffPartner(j).resSeq;
    NumCode_H3(j)                           = pathway.spinDiffPartner(j).NumCode;
    resName_H3{j}                       	= pathway.spinDiffPartner(j).resName;      	% convert char to double before comparing
    
    rho_H3                                  = pathway.spinDiffPartner(j).rho_H3;        % rho spindiff partner
    sigma_H1H3                              = pathway.spinDiffPartner(j).sigma_H1H3;    % sigma13
    sigma_H3H2                              = pathway.spinDiffPartner(j).sigma_H3H2;    % sigma32
           
    three_spins_build_up_H1H2               = [];
%-----------------------------------------------
%The Solomon differential equation
        
%   	gamma                                   = [rho_H1 sigma_H1H3 sigma_H1H2; sigma_H1H3 rho_H3 sigma_H3H2; sigma_H1H2 sigma_H3H2 rho_H2];
% %     options                                 = [];   
% %    	[t,three_spins_build_up_H1H2]           = ode45(@Solomon,time,[1 0 0],options,gamma);


    RELAX3s      	= [rho_H1 sigma_H1H2 sigma_H1H3; sigma_H1H2 rho_H2 sigma_H3H2; sigma_H1H3 sigma_H3H2 rho_H3];    
    C               = [];
    build_up_H1H2   = [];
  	nbr_mixing_time = max ( size(time) );

 	for k = 1 : nbr_mixing_time
        C = expm( -time(k)*RELAX3s );
       	build_up_H1H2 = [ build_up_H1H2;  C(2,1) ];
    end
    
    three_spins_build_up_H1H2 = build_up_H1H2';
    
    [two_spins_build_up_H1H2_scaled(:,j),scaleFactor(j)]      = scale2labSchemeFun(two_spins_build_up_H1H2',NumCode_H3(j),resName_H3{j},conditions,user);
    [three_spins_build_up_H1H2_scaled(:,j),scaleFactor(j)]    = scale2labSchemeFun(three_spins_build_up_H1H2(:,3),NumCode_H3(j),resName_H3{j},conditions,user);    

    
end

    two_spins_build_up_H1H2_sel                 =  two_spins_build_up_H1H2;  %  selBuildupAtMixtimes(two_spins_build_up_H1H2, time, mixtimes);

for j=1:partnerCount
    
    two_spins_build_up_H1H2_scaled_single       = two_spins_build_up_H1H2_scaled(:,j);
    two_spins_build_up_H1H2_scaled_sel(:,j)     = two_spins_build_up_H1H2_scaled(:,j); %selBuildupAtMixtimes(two_spins_build_up_H1H2_scaled_single, time, mixtimes);
    
    three_spins_build_up_H1H2_scaled_single     = three_spins_build_up_H1H2_scaled(:,j);
    three_spins_build_up_H1H2_scaled_sel(:,j)  	= three_spins_build_up_H1H2_scaled(:,j); %selBuildupAtMixtimes(three_spins_build_up_H1H2_scaled_single, time, mixtimes);
    
end
    
     sum(two_spins_build_up_H1H2_scaled,2);
     
     sum(three_spins_build_up_H1H2_scaled,2);
    
% % % if ATOMi_resSeq == 16 && ATOMi_NumCode == 1 ...
% % %     && ATOMj_resSeq == 22 && ATOMj_NumCode == 2
% % %         
% % %          two_spins_build_up_H1H2_sel
% % %          two_spins_build_up_H1H2_scaled_sel
% % %          three_spins_build_up_H1H2_scaled_sel
% % %          return
% % %         
% % % end
    
    
if length(pathway.spinDiffPartner) >= 1
        
    simul.ATOMi_resSeq                      = ATOMi_resSeq; 
    simul.ATOMi_NumCode                     = ATOMi_NumCode;      
    simul.ATOMi_HpdbSerNo                   = ATOMi_HpdbSerNo;     
    simul.ATOMj_resSeq                      = ATOMj_resSeq;    
    simul.ATOMj_NumCode                     = ATOMj_NumCode;   
    simul.ATOMj_HpdbSerNo                   = ATOMj_HpdbSerNo; 

    simul.mixtimes                          = mixtimes; 
    simul.two_spins_build_up_H1H2           = two_spins_build_up_H1H2_sel';
    simul.two_spins_build_up_H1H2_scaled_sum    = sum(two_spins_build_up_H1H2_scaled,2);
    simul.three_spins_build_up_H1H2_scaled_sum  = sum(three_spins_build_up_H1H2_scaled,2);
    simul.rho_H1                            = rho_H1;
    simul.rho_H2                            = rho_H2;
    simul.rPdb_H1H2                         = rPdb_H1H2;
    simul.sigma_H1H2                        = sigma_H1H2;
    simul.partnerCount                      = partnerCount;
    simul.resSeq_H3                         = resSeq_H3;
    simul.NumCode_H3                        = NumCode_H3;
    simul.resName_H3                        = resName_H3;
    simul.scaleFactor                       = scaleFactor;

    
else
    
    simul.ATOMi_resSeq                      = []; 
    simul.ATOMi_NumCode                     = [];      
    simul.ATOMi_HpdbSerNo                   = [];     
    simul.ATOMj_resSeq                      = [];    
    simul.ATOMj_NumCode                     = [];   
    simul.ATOMj_HpdbSerNo                   = []; 
    
    simul.mixtimes                          = [];
    simul.two_spins_build_up_H1H2           = [];
    simul.two_spins_build_up_H1H2_scaled    = [];
    simul.three_spins_build_up_H1H2_scaled  = [];
    simul.rho_H1                            = [];
    simul.rho_H2                            = [];
    simul.rPdb_H1H2                         = [];
    simul.sigma_H1H2                        = [];
    simul.partnerCount                      = [];
    simul.resSeq_H3                         = [];
    simul.NumCode_H3                        = [];
    simul.resName_H3                        = [];
    simul.scaleFactor                       = [];
    
end
          

    
  