/***************************************************************************
                          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;

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

  double r1, r2, r3, r1act, r1acttag, r1actdeg, endofcycle, delta;
  double K, L, M, N, genrate, r1genrate2, genbasal, decrate, kon, koff, tagrate, untagrate;
  char c='n';  

  r1=0.0;
  r2=0.0;
  r3=0.0;
  r1act=0.0;
  r1acttag=0.0;
  r1actdeg=0.0;
  endofcycle=0.0;

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

  // 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);

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

  printf("\nBeginn der Simulation, Start der Zeitschleife!\n");
  for (step=-200000; step < 200000; step=step++)         // long-time: 1500000
    {
      // protein r1actdeg synthesis
      if (step>10000)
         r1actdeg=r1actdeg+2.5*genrate;
      else
         r1actdeg=r1actdeg+genbasal;
      // and degradation
      r1actdeg=r1actdeg-r1actdeg*10*decrate;

      // r1act synthesis
      r1act=r1act+genrate/10;
      // protein r1act tagging for rapid degradation
      delta=tagrate*r1act*r1actdeg-untagrate*r1acttag;
      r1acttag=r1acttag+delta;
      r1act=r1act-delta;
      // and degradation
      r1act=r1act-r1act*decrate;
      r1acttag=r1acttag-r1acttag*100*decrate;// more rapid destruction of "tagged" form

      // cyclin synthesis - reciprocal inhibition, 1/(1+pow(x/Km,Hill coeff))
      r1=r1+r1genrate2/(1+pow(L/r1act,4)); 
      if ((c=='y'))
        {
         r1=r1+genbasal+genrate/(1+pow(r3/K,2.1));
         r2=r2+genbasal+genrate/(1+pow(r1/K,2.1));
         r3=r3+genbasal+genrate/(1+pow(r2/K,2.1)); 
        }
      else
        {
         r1=r1+genbasal+genrate/(1+pow(r3/K,2.1))/(1+pow(endofcycle/N,1.4));
         r2=r2+genbasal+genrate/(1+pow(r1/K,2.1))/(1+pow(endofcycle/N,1.4));
         r3=r3+genbasal+genrate/(1+pow(r2/K,2.1))/(1+pow(endofcycle/N,1.4));
        }
      // endofcycle generation - logical NOR(r2,r3)
      endofcycle=endofcycle+genbasal+10*genrate/(1+pow(r2/M,3))/(1+pow(r3/M,3));
      // and degradation
      endofcycle=endofcycle-endofcycle*50*decrate;

      // cyclin degradation
      r1=r1-r1*decrate;  // v2 is .75*r1*decrate
      r2=r2-r2*decrate;
      r3=r3-r3*decrate;

      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/K*3),(long)(r2/K*3),(long)(r3/K*3),(long)(r1actdeg/K/30),(long)(r1act/K/30),(long)(r1acttag/K),(long)(endofcycle/K/30));
       }

    }

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

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

  system("gnuplot plot1.gp");

  return EXIT_SUCCESS;

}