#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>


#define dim 22


#define abt 1 /*sampling rate*/
#define N 70000 /*number of iterations*/ 


      /***********************************************************/
      /********************** PARAMETERS *************************/
      /***********************************************************/
      
    double mu                        = 0.45;
    double theta_HSL                 = 0.24;
    double theta_IPTG                = 0.0137;
    
    double rho_production            = 10;
    double rho_latch                 = 10;
    double rho_lysis                 = 10;
    double rho_antilysis             = 10;
    double rho_qs                    = 10;
    double rho_luxi                  = 10;
    
    
    double beta_production_max       = 0.44;
    double beta_production_min       = 0.013;
    double beta_production_max_luxi  = 0.002 ;
    double beta_production_min_luxi  = 0.00001;
    double beta_latch                = 0.28;
    double beta_lysis                = 0.0426;
    double beta_antilysis            = 0.0066;
    double beta_qs                   = 0.018;
    
    double K_laci                    = 700;  //guessed
    double K_ci                      = 7000; //guessed
    double K_tetr                    = 7000; //guessed
    double K_P                       = 700;  //guessed
    double K_IPTG_laci               = 1200; //guessed
    
    double n_laci                    = 2;
    double n_ci                      = 2;
    double n_tetr                    = 3;
    double n_P                       = 2; 
    
    double gamma_protein             = 0.1;
    
    double k_on_iptg_laci            = 0.001; 
    double k_off_iptg_laci           = 0.033; 
    
    double k_on_P                    = 0.00001;
    double k_off_P                   = 0.00333;
    
    double alpha_mrna                = 0.006;
    double alpha_X                   = 0.0000080225;
    double alpha_Y                   = 0.0000080225;
    double alpha_ci                  = 0.002888; // taged
    double alpha_laci                = 0.002888;// taged
    double alpha_tetr                = 0.002888; // taged
    double alpha_holin               = 0.0002 ; //guessed
    double alpha_antiholin           = 0.0002 ; //guessed
    double alpha_endolysin           = 0.0002; // guessed
    double alpha_luxi                = 0.00288811; // taged
    double alpha_luxr                = 0.0002;
    
    double alpha_hsl                 = 0.00017;
    
    double alpha_p                   = 0.0002;


const int PoissonRandomNumber(const double lambda);

main()
{

 srand(time(0));  //initialize random generator
 
  remove("stochastic.dat");
  double a;
  int i,j,k,b;
  float y[dim];
  double d1,d2,d3,d4,d5, d6;
  double lambda, tau, t;
 
 FILE *outfile;
 outfile=fopen("stochastic.dat","w");
 
// fprintf(outfile1,"prod1            latch2            lysis3          anti4        qs5                  X6            Y7                   ci8            laci9            tetr10            holin11         endo12             antiho13             luxi14               luxr15            hsl16              P17           iptg18         hsl_outside19        iptg_outside mrnaluxi\n");


  //initializing
for (b=0; b<1; b++) 
{ t=0;
  for(i=0; i<=21; i++)
  {
      y[i]=0.0;
  }  
   
  k=0;
  j=0;
  for(j=1;j<N-1;j++)
  {
         k++; 
    //INPUT SIGNALS
    
    // IPTG
   if (j>15000)
  {y[18]=50000; 
  }
    //HSL_outside
  if(j>25000)
  {y[19]=100;
  }
 if(j>35000)       //turns off IPTG
  {y[18=0;
  }

 
                  
    //definestau
    tau=0.5; //1.*rand()/((RAND_MAX));
 


    
/***********************************************************/
/******************** MRNA EQUATIONS ***********************/
/***********************************************************/
      
      
     
   //  mRNA_production - first term of equation
    lambda=rho_production*(1/(1+pow((y[8]/(K_laci*(1+pow(y[17]/K_IPTG_laci,n_laci)))),n_laci)))*((beta_production_max*pow(y[16],n_P)/(pow(y[16],n_P)+pow(K_P,n_P)))+beta_production_min)*tau;
    d1=PoissonRandomNumber(lambda);
     
     //mRNA_production - second term of equation
    lambda=(alpha_mrna+(log(2)/1800))*y[0]*tau;
     d2=PoissonRandomNumber(lambda);
    
     //mRNA_production - equation
    y[0]=y[0]+d1-d2;
     
      //mRNA_latch - first term of the equation
     lambda=rho_latch*beta_latch*(1/(1+pow((y[7]/K_ci),n_ci)))*tau;
     d1=PoissonRandomNumber(lambda);
      
      //mRNA_latch - second term of the euqation
     lambda=(alpha_mrna+(log(2)/1800))*y[1]*tau;
     d2=PoissonRandomNumber(lambda);
      
      //mRNA_latch - equation
     y[1]=y[1]+d1-d2;
      
      //mRNA_lysis - first term of the equation
    lambda=rho_lysis*beta_lysis*(1/(1+pow((y[9]/K_tetr),n_tetr)))*tau;
      d1=PoissonRandomNumber(lambda);
      
      //mRNA_lysis - second term of the euqation
      lambda=(alpha_mrna+(log(2)/1800))*y[2]*tau;
      d2=PoissonRandomNumber(lambda);
      
      //mRNA_lysis - equation
      y[2]=y[2]+d1-d2;
      
      //mRNA_antilysis - first term of the equation
      lambda=rho_antilysis*beta_antilysis*tau;
      d1=PoissonRandomNumber(lambda);
      
      //mRNA_antilysis - second term of the equation
      lambda=(alpha_mrna+(log(2)/1800))*y[3]*tau;
      d2=PoissonRandomNumber(lambda);
      
      //mRNA_antilysis - equation
      y[3]=y[3]+d1-d2;
      
      //mRNA_QS - first term of the equation
      lambda=rho_qs*beta_qs*tau;
      d1=PoissonRandomNumber(lambda);
      
      //mRNA_QS - second term of the equation
      lambda=(alpha_mrna+(log(2)/1800))*y[4]*tau;
      d2=PoissonRandomNumber(lambda); 
 
       //mRNA_QS - equation
      y[4]=y[4]+d1-d2;
      
      //mRNA_luxi - first term of the equation 
      lambda=rho_luxi*((beta_production_luxi*pow(y[16],n_P)/(pow(y[16],n_P)+pow(K_P,n_P)))+beta_production_min_luxi)*tau;
      d1=PoissonRandomNumber(lambda);
      
      //mRNA_luxi - second term of the equation
      lambda=(alpha_mrna+(log(2)/1800))*y[20]*tau;
      d2=PoissonRandomNumber(lambda);
      
      //mRNA_luxi - equation
      y[20]=y[20]+d1-d2;
      
 
      /***********************************************************/
      /****************** PROTEIN EQUATIONS **********************/
      /***********************************************************/
      
         
      //[X] - first term of the equation
      lambda=gamma_protein*y[0]*tau;
      d1=PoissonRandomNumber(lambda);
      
      //[X] - second term of the equation
      lambda=(alpha_X+(log(2)/1800))*y[5]*tau;
      d2=PoissonRandomNumber(lambda);
      
      //[X] - equation
      y[5]=y[5]+d1-d2;
    
      //[Y] - first term of the equation
      lambda=gamma_protein*y[0]*tau;
      d1=PoissonRandomNumber(lambda);
      
      //[Y] - second term of the equation
      lambda=(alpha_X+(log(2)/1800))*y[6]*tau;
      d2=PoissonRandomNumber(lambda);
      
      //[Y] - equation
      y[6]=y[6]+d1-d2;
   
       //[cI] - first term of the equation
      lambda=gamma_protein*y[0]*tau;
      d1=PoissonRandomNumber(lambda);
      
      //[cI] - second term of the equation
      lambda=(alpha_ci+(log(2)/1800))*y[7]*tau;
      d2=PoissonRandomNumber(lambda);
      
      //[cI] - equation
      y[7]=y[7]+d1-d2;
      
       //[LacI] - first term of the equation
      lambda=gamma_protein*y[1]*tau;
      d1=PoissonRandomNumber(lambda);
      
      //[LacI] - second term of the equation
      lambda=(alpha_laci+(log(2)/1800))*y[8]*tau;
      d2=PoissonRandomNumber(lambda);
      
      //[LacI] - equation
      y[8]=y[8]+d1-d2;
      
      //[TetR] - first term of the equation
      lambda=gamma_protein*y[1]*tau;
      d1=PoissonRandomNumber(lambda);
      
      //[TetR] - second term of the equation
      lambda=(alpha_tetr+(log(2)/1800))*y[9]*tau;
      d2=PoissonRandomNumber(lambda);
      
      //[TetR] - equation
      y[9]=y[9]+d1-d2;
      
      //[Holin] - first term of the equation
      lambda=gamma_protein*y[2]*tau;
      d1=PoissonRandomNumber(lambda);
      
      //[Holin] - second term of the equation
      lambda=(alpha_holin+(log(2)/1800))*y[10]*tau;
      d2=PoissonRandomNumber(lambda);
      
      //[Holin] - equation
      y[10]=y[10]+d1-d2;
      
       //[Endolysin] - first term of the equation
      lambda=gamma_protein*y[2]*tau;
      d1=PoissonRandomNumber(lambda);
      
      //[Endolysin] - second term of the equation
      lambda=(alpha_endolysin+(log(2)/1800))*y[11]*tau;
      d2=PoissonRandomNumber(lambda);
      
      //[Endolysin] - equation
      y[11]=y[11]+d1-d2;
      
      //[Antiholin] - first term of the equation
      lambda=gamma_protein*y[3]*tau;
      d1=PoissonRandomNumber(lambda);
      
      //[Antiholin] - second term of the equation
      lambda=(alpha_antiholin+(log(2)/1800))*y[12]*tau;
      d2=PoissonRandomNumber(lambda);
      
      //[Antiholin] - equation
      y[12]=y[12]+d1-d2;
      
      //[LuxI] - first term of the equation 0.6 value for the scheint orgaon sequence
      lambda=0.6*gamma_protein*y[20]*tau;
      d1=PoissonRandomNumber(lambda);
      
      //[LuxI] - second term of the equation
      lambda=(alpha_luxi+(log(2)/1800))*y[13]*tau;
      d2=PoissonRandomNumber(lambda);
      
      //[LuxI] - equation
      y[13]=y[13]+d1-d2;
      
      //[LuxR] - first term of the equation
      lambda=gamma_protein*y[4]*tau;
      d1=PoissonRandomNumber(lambda);
      
      //[LuxR] - second term of the equation
      lambda=k_off_P*y[16]*tau;
      d2=PoissonRandomNumber(lambda);
      
      //[LuxR] - third term of the equation
      lambda=k_on_P*(y[15]*y[14])*tau;
      d3=PoissonRandomNumber(lambda);
      
      //[LuxR] - fourth term of the equation
      lambda=(alpha_luxr+(log(2)/1800))*y[14]*tau;
      d4=PoissonRandomNumber(lambda);
      
      //[LuxR] - equation
      y[14]=y[14]+d1+d2-d3-d4;
      
      /***********************************************************/
      /****************** COMPLEX EQUATIONS **********************/
      /***********************************************************/
      
      //[HSL] - first term of the equation
      lambda=mu*y[13]*tau;
      d1=PoissonRandomNumber(lambda);
      
      //[HSL] - second term of the equation
      lambda=(k_off_P+alpha_luxr)*y[16]*tau;
      d2=PoissonRandomNumber(lambda);
      
      //[HSL] - third term of the equation
      lambda=k_on_P*y[15]*y[14]*tau;
      d3=PoissonRandomNumber(lambda);
      
      //[HSL] - fourth term of the equation
      lambda=theta_HSL*y[19]*tau;
      d4=PoissonRandomNumber(lambda);
      
       //[HSL] - fifth term of the equation
      lambda=theta_HSL*y[15]*tau;
      d5=PoissonRandomNumber(lambda);
      
      //[HSL] - sixth term of the equation
      lambda=(alpha_hsl+(log(2)/1800))*y[15]*tau;
      d6=PoissonRandomNumber(lambda);
      
      //[HSL] - equation
      y[15]=y[15]+d1+d2-d3+d4-d5-d6; 
      
      //[P] - first term of the equation
      lambda=k_on_P*y[15]*y[14]*tau;
      d1=PoissonRandomNumber(lambda);
      
      //[P] - second term of the equation
      lambda=k_off_P*y[16]*tau;
      d2=PoissonRandomNumber(lambda);
      
      //[P] - third term of the equation
      lambda=(alpha_p+(log(2)/1800))*y[16]*tau;
      d3=PoissonRandomNumber(lambda);
     
      //[P] - equation
      y[16]=y[16]+d1-d2-d3;
      
      //[IPTG] - first term of the equation
      lambda=theta_IPTG*y[18]*tau;
      d1=PoissonRandomNumber(lambda);
    
    //[IPTG] - second term of the equation
      lambda=theta_IPTG*y[17]*tau;
      d2=PoissonRandomNumber(lambda);
      
      //[IPTG] - equation
      y[17]=y[17]+d1-d2;
      
 
    
    
 
   fprintf(outfile,"%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",t,y[0],y[1],y[2],y[3],y[4],y[5],y[6],y[7],y[8],y[9],y[10],y[11],y[12],y[13],y[14],y[15],y[16],y[17],y[18],y[19],y[20]);    

}
}

 fclose(outfile);
 
  return 0;
}




const int PoissonRandomNumber(const double lambda)
{
  int k=0;                          //Counter
  const int max_k = 2*lambda;           //k upper limit
  double p = 1.*rand()/((RAND_MAX)); //uniform random number
  double P = exp(-lambda);          //probability
  double sum=P;                     //cumulant
  if (sum>=p) return 0;             //done allready
  for (k=1; k<max_k; ++k) {         //Loop over all k:s
    P*=lambda/(double)k;            //Calc next prob
    sum+=P;                         //Increase cumulant
    if (sum>=p) break;              //Leave loop
  }

  return k;                         //return random number
}