%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%       Subroutine of sigmaTrue     %%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [X Y out work] = plotSexpStrueFun(sigmaData, out, work, user, conditions)

X = [];
Y = [];

if isfield(out, 'slopeSexpStrue')
    out = rmfield(out, 'slopeSexpStrue');
end
    

out.noTrueSigmaValue                = 0;
out.plottedPointsTrueSigmaValue     = 0;


for j=1:length(sigmaData)
    plotData(j).avgSigma = sigmaData(j).avgSigma;
end

if user.plotIsoAnisoSigma == 1
    subplot(1,1,1); hold; xlabel('sIsoExp'); ylabel('sTrue');
            
            w =  sprintf('sIsoExp vs. sTrue');
            title(w);
end

if user.plotIsoAnisoSigma == 2
    subplot(1,1,1); hold; xlabel('sAnIsoExp'); ylabel('sTrue');
            
            w =  sprintf('sAnIsoExp vs. sTrue');
            title(w);
end


for j = 1:length(sigmaData)
                
      if sigmaData(j).unique == 2               %The first one of the partners
          
          bothAtoms = [sigmaData(j).numcode.ATOMi, sigmaData(j).numcode.ATOMj];
            
          if ~isempty(sigmaData(j).trueSigma) %isnumeric(sigmaData(j).trueSigma) %&& isempty(sigmaData(j).trueSigma)
             
              if size(find(ismember(bothAtoms, [1, 2])), 2) == 2
                    X(j) = plotData(j).avgSigma;
                    Y(j) = sigmaData(j).trueSigma;
            
                    if strcmp(user.jindx, 'n') == 1   
                    
                        plot(X(j),Y(j),'ro','MarkerSize',5);
                            
                    elseif strcmp(user.jindx, 'y') == 1
                        text(X(j),Y(j),num2str(j),'color','r','FontSize',8,'HorizontalAlignment','center');
                    end
                    
                    out.plottedPointsTrueSigmaValue =   out.plottedPointsTrueSigmaValue + 1;
              end
          end
      end
end
  
  out.slopeSexpStrue = X / Y;
      
  if min(X) < min(Y)
        minvalue = min(X);
  else
        minvalue = min(Y);
  end
  
  if max(X) > max(Y)
        maxvalue = max(X);
  else
        maxvalue = max(Y);
  end
 
  axis([minvalue maxvalue  minvalue  maxvalue])
  
  hline = refline(1,0);
  set(hline,'Color','k','LineWidth',0.5)
      
  if user.plotIsoAnisoSigma == 1
      
      out.TAUCRATIO         = user.tauC/user.tauCsigmaSup;
      out.TAUCCORRFACTOR    = (out.TAUCRATIO - out.slopeSexpStrue)/ out.TAUCRATIO;
      out.PERCENT_ERROR     = abs(1 - out.slopeSexpStrue) * 100;
      out.CORRECT_TC_VALUE  = user.tauC * (1 - out.TAUCCORRFACTOR);
      
  diary('on')    
  
  fprintf(['\n\n\n\nFor sExp vs sTrue:\n\n', ...
         'Number of plotted points:             ', num2str(out(1).trueSigma),       ' \n\n', ...
         'TauC current:                         ', num2str(user.tauC),              ' \n\n',...
         'TauC sigma supplied:                  ', num2str(user.tauCsigmaSup),      ' \n\n',...
         'Ratio TauCs:                          ', num2str(out.TAUCRATIO),          ' \n\n',...
         'Fit slope:                            ', num2str(out.slopeSexpStrue),     ' \n\n',...
         'The percent error is:                 ', num2str(out.PERCENT_ERROR),      ' \n\n',...
         'The correction factor for TauC is:    ', num2str(out.TAUCCORRFACTOR),     ' \n\n',...
         'The correct TauC value should be:     ', num2str(out.CORRECT_TC_VALUE),   ' \n\n'])
     
  diary('off')    
  end