#include "runsim.h"
	
void create_stoich_matrix(
     struct array_sizes   *arsz_ptr, 
     struct rxn           *rxn_ptr, 
     int                  **stoich_mat )
{
     struct rxn       arxn;
     int numrxns = arsz_ptr -> num_rxns;
     int numcmpds = arsz_ptr -> num_cmpds;
     int i, mol_id=0, hor;
     int numsubs, numprods;
     int is, ip;
     
     for (i=0; i<numrxns; i++)
       {
          arxn = rxn_ptr[i];
       
          switch (arxn.rxn_index) 
          {
           case UNI_UNI:
            { 
           mol_id = arxn.substrates[0];
           stoich_mat[i][mol_id] = -1;
           mol_id = arxn.products[0];    
           stoich_mat[i][mol_id] = 1;
            break;
            }
          case UNI_MULTI: 
           {
           mol_id = arxn.substrates[0];
           stoich_mat[i][mol_id] = -1;
           numprods = arsz_ptr -> num_prods[i];
	    	if (numprods != 0) 
		    	for (ip = 0; ip < numprods; ip++) 
                {
                 mol_id = arxn.products[ip];
                 stoich_mat[i][mol_id] += arxn.numprods_made[ip];
                }
             break;
            }
              case BINDING:
              case HOMO_LTOR: 
             {
              mol_id = arxn.substrates[0];
              stoich_mat[i][mol_id] += -1;
              mol_id = arxn.substrates[1];
              stoich_mat[i][mol_id] += -1;
              mol_id = arxn.products[0];
              stoich_mat[i][mol_id] = 1;
              break;
            }
             case UNBINDING:
            {
              mol_id = arxn.substrates[0];
              stoich_mat[i][mol_id] = -1;
              mol_id = arxn.products[0];
              stoich_mat[i][mol_id] += 1;
              mol_id = arxn.products[1];
              stoich_mat[i][mol_id] += 1;
              break; 
            }
              case NTH_ORDER:
              case GENM:
              case TRANSV1:
              case UNDEFINED:
	
              {
		numsubs = arsz_ptr -> num_subs[i];
		for (is = 0; is < numsubs; is++) 
                  {
			mol_id = arxn.substrates[is];
			stoich_mat[i][mol_id] -= arxn.numsubs_used[is];
		  }
		numprods = arsz_ptr -> num_prods[i];
		if (numprods != 0) 
                    {
			for (ip = 0; ip < numprods; ip++)  
                         {
				mol_id = arxn.products[ip];
				stoich_mat[i][mol_id] += arxn.numprods_made[ip];
			}
		  }
              break;
	     }
                                                 default:
                 printf("a specified reaction type is not defined - calc_stoich_mat\n");
     }
     }
     /*  Estimate the Highest Order of Reaction for each of the compounds.   
         It is needed for tau-step selection.  
                 See Cao, Gillespie, Petzold, J.Chem.Phys, 124,044109 (2006)*/
        for (ip = 0; ip < numcmpds; ip++) 
        {
               hor = 0;                 
		   for (is = 0; is < numrxns; is++)
           {
                    if (stoich_mat[is][ip] < 0)
                    {  
                        mol_id = 0;
                               for (i = 0; i < numcmpds; i++)
                                   if (stoich_mat[is][i] < 0)
                                   {
                                          mol_id += 1;
                                   }
                    }
                    if (mol_id > hor)
                    {
                         hor = mol_id;
                         mol_id = 0;
                    }
           }
           stoich_mat[numrxns][ip] = hor;   
        }
     


/*     for (is = 0; is < numrxns+1; is++)  
      		for (ip = 0; ip < numcmpds; ip++) 
              if (stoich_mat[is][ip] != 0)
                     	printf("\t %d - %d, STH = %d\n", is, ip, stoich_mat[is][ip]);
*/
     
}
     
