
#include "runsim.h"

double find_tau(
           struct  array_sizes          *arsz_ptr,
		   struct  rxn 					*rxn_ptr,
		   struct  cmpd_data 			*cmpds,
		   int						    **stoich_mat,
           double                       eps
		   )
{           	
    int const numrxns = arsz_ptr -> num_rxns;
    int const numcmpds = arsz_ptr -> num_cmpds;
    double mu, sigma, change, tau;
    int i,j;
    double tau1 = 1000, tau2 = 1000;
    double tau01, tau02;
    
    for (i=0; i<numcmpds; i++)
    {

        if (stoich_mat[numrxns][i] >0 )
        {
        mu = 0;
        sigma = 0;
            for (j=0; j<numrxns; j++)
            {
                mu += stoich_mat[j][i] * rxn_ptr[j].ai;
                sigma +=  stoich_mat[j][i] * stoich_mat[j][i] * rxn_ptr[j].ai;   
            }
                if (mu == 0) mu = 1.0e-06;
                if (sigma == 0) sigma = 1.0e-06;  
        change = eps * cmpds[i].value /stoich_mat[numrxns][i];
                if (change < 1) change = 1;
        tau01 = change / fabs(mu);
                if (tau01 < tau1) tau1 = tau01;
        tau02 = change * change / sigma;
                if (tau02 < tau2) tau2 = tau02;
        }
    }

    if (tau1 > tau2 ) 
        tau = tau2;
    else
        tau = tau1;
    
    return tau;
}
