Self Modeling Curve Resolution (SMCR) Procedure

// Input: datamat is a 2D matrix whose columns are spectra obtained from two-component 
// solutions of different concentration. For example, the first spectrum can be that of a 
// pure solvent and the remaining spectra can have various (at least one) solute
// concentrations.
// There is no need to normalize the spectral measurements -- this is done
// inside the function.
// Output: component1 and component2, one of which is the minimum area solute correlated
// spectrum, while the other is the minimum area spectrum pertaining to the pure solvent.
// The pure solvent spectrum itself can be re-generated from the SMCR scores pertaining
// to the pure solvent spectrum. For example, if the first spectrum (in column 0) is the pure 
// solvent then the following command will re-generate the pure solvent spectrum (solvent).
// from the SMCR scores and PC loadings.
// solvent=score1[0]*PCLoadings1+score2[0]*PCLoadings2
// Requires: IP 6.22 or newer.
// Reference: Lawton and Sylvestre, Technometrics, volume 13, page 617, 1971
// Procedure written by: David Wilcox (wilcoxds at purdue dot edu)
 
 
Function SMCR(datamat)
          wave datamat
  
        // the following command normalizes and transposes the input matrix.
        MatrixOP/O/FREE inMat=scaleCols(datamat,(powr(sumcols(datamat),-1)^t))^t
        Variable i,numRows
        
           //Calculate scores and loadings from SVD 
          MatrixSVD inMat                              
          wave M_U, W_W, M_VT 
          MatrixTranspose M_VT
          MatrixOP/O W_W=diagonal(W_W)
          MatrixOP /FREE/O PCLoadings1=col(M_VT,0)
          MatrixOP/FREE/O PCLoadings2=col(M_VT,1)
          MatrixOP/FREE/O aug_scores = M_U x W_W
          MatrixOP/FREE/O score1=col(aug_scores ,0)     // PC1 scores for each input specxtrum
          MatrixOP/FREE/O score2=col(aug_scores ,1)     // PC2 scores for each input specxtrum
          
          //Find minimum positive and maximum negative of loadings ratio (slopes)
          MatrixOP/O/FREE Lratio = PCLoadings1/PCLoadings2
          numRows=DimSize(Lratio,0)
          variable m1,m2
          i=0
          m1=-10^100
          m2=10^100
          do
                if (Lratio[i] > 0)
                m1 = max(-Lratio[i],m1)
            elseif (Lratio[i] < 0)
                m2 = min(-Lratio[i],m2)
            endif
            i +=1
          while (i<numRows)
          
          // Linear regression of scores ratio
          CurveFit /Q line, score2/x=score1
          wave W_coef
          Variable m=W_coef[1]
          Variable b=W_coef[0]
          
          // Intersection of scores and loadings
          make/FREE/O/N=(2,2)  x1x2num={{b,1},{0,1}}
          make/FREE/O/N=(2,2)  x1y1denom={{-m,1},{-m1,1}}
          MatrixOP/FREE/O x1= Det(x1x2num)/Det(x1y1denom) // (x1,y1) is one end point of the PC score line
          make/FREE/O/N=(2,2)  y1num={{-m,b},{-m1,0}}
          MatrixOP/FREE/O y1= Det(y1num)/Det(x1y1denom)
          make/FREE/O/N=(2,2)  x2y2denom={{-m,1},{-m2,1}}
          MatrixOP/FREE/O x2= Det(x1x2num)/Det(x2y2denom) // (x2,y2) is one end point of the PC score line
          make/FREE/O/N=(2,2)  y2num={{-m,b},{-m2,0}}
          MatrixOP/FREE /O y2  = Det(y2num)/Det(x2y2denom) 
          
          // Make pure component spectra
          Variable x1v=x1[0],x2v=x2[0],y1v=y1[0],y2v=y2[0]
          MatrixOP/O    component1 = x1v*PCLoadings1 + y1v*PCLoadings2
          MatrixOP/O    component2 = x2v*PCLoadings1 + y2v*PCLoadings2
          KillWaves/Z M_U,W_W,M_VT,W_coef,W_sigma
end