/***************************************************************************
                          main.c  -  description
                             -------------------
    begin                : Don Jun 30 01:16:19 CEST 2005
    copyright            : (C) 2005 by 
    email                : 
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/

#include <math.h> // for pow

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

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

int main(int argc, char *argv[])
{

  FILE *outhandle1;
  FILE *outhandle2;

  long step;
  const char * const outfilename1="data1.crv";
  const char * const outfilename2="data2.crv";

  double r1, r1tag, r2, r3, r1deginhib, r1deg, r1degcomplexr1deginhib, delta;
  double K, L, M, genrate, genbasal, decrate, kon, koff, tagrate, untagrate;
  double r3actdeg, r3act, r3acttag, r3genrate2;
  char c;  

  r1deginhib=0.0;
  r1deg=0.0;
  r1degcomplexr1deginhib=0.0;
  delta=0.0;
  r3actdeg=0.0;
  r3act=0.0;

  K=1E-4;
  L=K*12;
  M=K*2000; //1200 17580
  
  genrate=K;
  r3genrate2=genrate/500;
  genbasal=K/15000; // leaky promoter
  decrate=K;
  kon=genrate*10000;
  koff=decrate;
  tagrate=genrate*100;
  untagrate=decrate;

  r1=0;
  r1tag=0; // tagged for destruction
  r2=0;
  r3=0;
  
  // open files for write

  if ((outhandle1=fopen(outfilename1,"w"))==0)
    {
      printf("Error opening file %s for write!\n", outfilename1);
      return EXIT_FAILURE;
    }
  printf("File %s opened for write!\n", outfilename1);
  if ((outhandle2=fopen(outfilename2,"w"))==0)
    {
      printf("Error opening file %s for write!\n", outfilename2);
      return EXIT_FAILURE;
    }
  printf("File %s opened for write!\n", outfilename2);

  printf("Multiple cycles (y/n)? ");
  scanf("%c",&c);

  printf("Beginn der Simulation, Start der Zeitschleife!\n");
  for (step=-200000; step < 200000; step=step++)         // long-time: 1500000
    {
      // protein synthesis
      r1deg=r1deg+genrate;
      if (step>10000)
          r1deginhib=r1deginhib+2.5*genrate;
      else
          r1deginhib=r1deginhib+genbasal;
       
      // chemical reaction of stoichiometric inhibitor binding
      delta=kon*r1deg*r1deginhib-koff*r1degcomplexr1deginhib;
      r1degcomplexr1deginhib=r1degcomplexr1deginhib+delta;
      r1deg=r1deg-delta;
      r1deginhib=r1deginhib-delta;

      // decay of proteins  
      r1deginhib=r1deginhib-r1deginhib*10*decrate;
      r1deg=r1deg-r1deg*10*decrate;
      r1degcomplexr1deginhib=r1degcomplexr1deginhib-r1degcomplexr1deginhib*10*decrate;      
 
      r1=r1+genbasal+genrate/(1+(pow(r3/K,2.1)));
      r2=r2+genbasal+genrate/(1+(pow((r1+r1tag)/K,2.1)));
      r3=r3+genbasal+genrate/(1+(pow(r2/K,2.1))); 
      if (c!='y') r3=r3+r3genrate2/(1+(pow(M/r3act,4))); 
      delta=tagrate*r1*r1deg-untagrate*r1tag;
      r1tag=r1tag+delta;
      r1=r1-delta;

      r1=r1-r1*decrate;  // v2 is .75*r1*decrate
      r1tag=r1tag-r1tag*100*decrate;// more rapid destruction of "tagged" form
      r2=r2-r2*decrate;
      r3=r3-r3*decrate;
      r3actdeg=r3actdeg+genbasal+genrate/(1+pow(r2/L,4))+genrate/(1+pow(r3/L,4));
      r3actdeg=r3actdeg-r3actdeg*decrate*10;

      r3act=r3act+genrate;
      delta=tagrate*r3act*r3actdeg-untagrate*r3acttag;
      r3acttag=r3acttag+delta;
      r3act=r3act-delta;

      r3act=r3act-r3act*decrate;
      r3acttag=r3acttag-r3acttag*100*decrate; //rapid degradation

      if ((step>-50000)&&(step%100==0))
       {
        fprintf(outhandle1, "%ld\t%ld\t%ld\t%ld\t%ld\t%ld\t%ld\t%ld\n",step,(long)((r1+r1tag)/K*3),(long)(r1tag/K*100),(long)(r2/K*3),(long)(r3/K*3),(long)(r1deginhib/K/10),(long)(r1deg/K/10),(long)(r1degcomplexr1deginhib/K/10));
        fprintf(outhandle2, "%ld\t%ld\t%ld\t%ld\t%ld\t%ld\t%ld\n",step,(long)((r1+r1tag)/K*3),(long)(r2/K*3),(long)(r3/K*3),(long)(r3actdeg/K/10),(long)(r3act/K/100),(long)(r3acttag/K/5));
       }
      
    }

  printf("Ende der Simulation!\n");

  fclose(outhandle1);
  printf("File %s closed!\n", outfilename1);

  fclose(outhandle2);
  printf("File %s closed!\n", outfilename2);

  system("gnuplot plot1.gp");
  system("gnuplot plot2.gp");

  return EXIT_SUCCESS;
}