/*
 *    This program evaluates the j-factor from a distribution of 
 *    products of the multimerization-cyclization reaction. 
 *    The value of initial concentration of the monomers and 
 *    the symmetry of the sticky ends should also be provided.
 *    alexei@crab.chem.nyu.edu / version 2000-03
 */


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

void readdata(char*);
char* readstr(char*, char*);
void solution(void);

void tick(double);

int N=99;
double dtime=0.00001;
double c[100],ddc[100],o[100],ddo[100],jf[100],cg[100],og[100];

void main(int argc, char **argv)
{  
   if(argc==1) 
   {  printf("Data file is missing! \nUsage: %s data.file\n",argv[0]); 
      exit(1);
   }

   readdata(argv[1]); 

   solution();
}

void readdata(char *fname)
{  FILE *file;
   int i;
   char str[60], *tail;

   if((file=fopen(fname,"rt"))==NULL) exit(0);
   
   for(i=1;i<=N;i++) cg[i]=og[i]=0.;
   cg[0]=1.;
   
   while(fgets(str,sizeof(str),file)==str)
     {
       if     ((tail=readstr(str,"Concentration"     ))) cg[0]*=atof(tail);
       else if((tail=readstr(str,"Time step"         ))) dtime=atof(tail);
       else if((tail=readstr(str,"Identical ends"    )) && atoi(tail)) cg[0]*=2.;
       else if((tail=readstr(str,"Products:"         )))
          { for(i=1; i<=N; i++)
              { fscanf(file,"%le %le ",&cg[i],&og[i]);   }
          }
     }
   fclose(file);

   for(i=0; cg[i]>0. || og[i]>0.; i++);
   N=i;
   
   printf("%.1e %.1e JF:", cg[0], dtime);

   putchar('\n');       
   
}

char* readstr(char *c, char *s)
{  while(*s==*c++) s++;
   if(*s) return NULL;
   return c;
}

void solution()
{  int j,i,cond;
   double *chck,check[]={ 0.95,0.9,0.85,0.8,0.7,\
	                      0.6,0.5,0.4,0.3,0.2,0.15,0.1,\
						  0.05,0.03,0.02,0.01,0 };    
   double diffoc,diffoo,diffc,diffnc,diffo,diffno;
   double dt;
   double coef, coef1,coef2;

 cond=-1;
 coef1=0; coef2=8*(1/cg[1]-1);   

 while(2*(coef2-coef1)/(coef1+coef2) >0.005)
 {
   coef=(coef1+coef2)/2;
   for(j=1;j<N;j++) jf[j]=(og[j]==0.)?0.:coef*j*og[j]/cg[j];
   for(j=1;j<N;j++) c[j]=ddc[j]=o[j]=ddo[j]=0.;
   c[0]=c[1]=1.;
   chck=check;
   dt=dtime;
   i=0;
   for(diffc=0.,diffo=0.,j=1;j<N;j++)
   { if(cg[j]>0.) diffc+=(j*c[j]-cg[j])*(j*c[j]-cg[j])/(cg[j]);
     if(og[j]>0.) diffo+=(j*o[j]-og[j])*(j*o[j]-og[j])/(og[j]);
   }
   diffoc=diffc;
   diffoo=diffo;

   while(*chck>0)
   {   i++;
       tick(dt);

       for(diffnc=0.,diffno=0.,j=1;j<N;j++)
	   { if(cg[j]>0.) diffnc+=(j*c[j]-cg[j])*(j*c[j]-cg[j])/(cg[j]);
         if(og[j]>0.) diffno+=(j*o[j]-og[j])*(j*o[j]-og[j])/(og[j]);
	   }

       if(diffnc > diffc && diffc < diffoc) 
	   { if(cond==-1) cond=1;  putchar('C'); break;}
       if(diffno > diffo && diffo < diffoo) 
	   { if(cond==-1) cond=0;  putchar('O'); break;}

	   diffoc=diffc;
	   diffoo=diffo;
	   diffc=diffnc;
	   diffo=diffno;
	   
	   if(c[1] > *chck) continue;
	   

       if(*chck==0.1) { dt*=10; i/=10;}
       if(*chck==0.01) { dt*=10; i/=10;}
       if(*chck==0.001) { dt*=10; i/=10;}

       chck++;
   }
   printf(" Results: [%.3e %.3e]\n", coef, i*dt);
   if(cond==0) coef2=coef;
   else        coef1=coef;
   cond=-1;
 }

 printf("      Open                     Closed                   j-factor\n");
 for(j=1; j<=N; j++)
    printf("%2d    %.3e : %.3e    %.3e : %.3e    %.3e\n", j, j*c[j], cg[j], j*o[j], og[j], 0.5*jf[j]*cg[0]);

}


/* 
 *  The following is a step by Euler scheme to solve 
 *  the differential equations of the multimerization-cyclization reaction. 
 */

void tick(double dt)
{ int i, j;
  double cc;

   ddc[0]=-c[0]*c[0]*dt;

   for(i=1; i<=N; i++)
     {  ddo[i]= jf[i]*c[i]*dt;
        ddc[0]-=ddo[i];
           for(cc=0,j=1;j<i;j++)  cc+=c[j]*c[i-j];
        ddc[i]= ( cc - 2*c[i]*c[0] - jf[i]*c[i] )*dt;
     }

   for(i=0;i<=N;i++)
     {  c[i]+=ddc[i];
        o[i]+=ddo[i];
     }
}