#include "runsim.h"

double calc_rxn_prob_erg(struct rxn *arxn, struct cmpd_data *cmpds, struct run_constants *r_c, struct array_sizes *arsz_ptr, int rxnnum)
{
	int i, id, numsubs;
	    double ai=0; 
        double order=0;
        double num_needed=0;
        double tot_num_avail=0;
        double num_sets_avail; 
        double tempval; 
	
		switch (arxn->rxn_index) {
			case BINDING:
			case UNIQUE_BI_BI:
			{
				ai = arxn->rates[0].value / r_c->volume;
				id = arxn->substrates[0]; 
				ai = ai * cmpds[id].value;
				id = arxn->substrates[1];       
				ai = ai * cmpds[id].value;
				break;
			}
			case UNBINDING:
			case UNI_UNI:
			case UNI_MULTI:
			{
				ai = arxn->rates[0].value;
				id = arxn->substrates[0];
				ai = ai * cmpds[id].value;
   				break;
			}
			case HOMO_LTOR:
			{
				ai = arxn->rates[0].value / r_c->volume;
				id = arxn->substrates[0];
				tot_num_avail = cmpds[id].value;
				ai = ai * tot_num_avail * (tot_num_avail - 1) / 2;
				break;
			}
			case TRANSV1:
			{
				numsubs = arsz_ptr -> num_subs[rxnnum];
				ai = arxn->rates[0].value / r_c->volume;
				num_needed = arxn->numsubs_used[0];
				id = arxn->substrates[0];
				tot_num_avail = cmpds[id].value;
				num_sets_avail = floor(tot_num_avail / num_needed);         
				if (num_sets_avail == 0) {
					ai = 0;
				}
				else {
					ai = ai * tot_num_avail;
				}
		
				tempval = 0;
				for (i = 1; i < numsubs; i++) {
					num_needed = arxn->numsubs_used[i];
					id = arxn->substrates[i];
					tot_num_avail = cmpds[id].value; 
					num_sets_avail = floor(tot_num_avail / num_needed);
					if (num_sets_avail == 0) {
						ai = 0;
						tempval = DBL_MAX;
						break;
					}
					else {
						tempval = tempval + (1 / tot_num_avail);
					}
				}
				ai = ai / tempval;
				break;
			}

			case GENM:
			{
				numsubs = arsz_ptr -> num_subs[rxnnum];
				ai = arxn->rates[0].value;
				num_needed = arxn->numsubs_used[0];
				id = arxn->substrates[0];
				tot_num_avail =  cmpds[id].value;
				num_sets_avail = floor(tot_num_avail / num_needed);
				if (num_sets_avail == 0) {
					ai = 0;
				}
				else {
					ai = ai * tot_num_avail / r_c->volume;
				}
		
				tempval = 1;
				for (i = 1; i < numsubs; i++) {
					num_needed = arxn->numsubs_used[i];
					id = arxn->substrates[i];
					tot_num_avail =  cmpds[id].value;
					num_sets_avail = floor(tot_num_avail / num_needed);
					if (num_sets_avail == 0) {
						tempval = 0;
						break;
					}
					else {
						tempval = tempval * tot_num_avail / (tot_num_avail + (arxn->rates[i].value * r_c->volume));
					}
				}
				ai = ai * tempval;
				break;
			}

			case NTH_ORDER:
			{
				numsubs = arsz_ptr -> num_subs[rxnnum];
				order = arxn->rates[1].value;
				ai = arxn->rates[0].value / pow(r_c->volume, (order - 1));
				if (order != 0) 
                {
					for (i = 0; i < numsubs; i++) 
                    {
						num_needed = arxn->numsubs_used[i];
						id = arxn->substrates[i];
						tot_num_avail = cmpds[id].value;                
						num_sets_avail = floor(tot_num_avail / num_needed);
    		
						if (num_sets_avail == 0) 
						                	ai = 0;
						else 
							ai = ai * pow(tot_num_avail, arxn->misc_subval[i]);
					}
				}
				break;
			}

			case UNDEFINED:
			{
				numsubs = arsz_ptr -> num_subs[rxnnum];
				ai = 1;
                                   
				for (i = 0; i < numsubs; i++) 
                {
					num_needed = arxn->numsubs_used[i];
					id = arxn->substrates[i];
					tot_num_avail = cmpds[id].value;                    
					num_sets_avail = floor(tot_num_avail / num_needed);
					if (num_sets_avail == 0) 
                    {
						ai = 0;	
						break;
					}
				}
				if (ai != 0)
				  ai = rxn_prob_undef(arxn, cmpds);
                              break;
			  }

                      default:            
                 printf("a specified reaction type is not defined - calc_rxn_prob\n");
		}  
        
	return ai;
}
