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

typedef double matrix[3][3];
typedef double vector[3];
typedef vector triplet[3];
typedef double converter[32768];

double dotprod(vector, vector);

struct Test {
       char   name[12];                     /* test name               */
       int    fr;                           /* frq. of testing & # of data */
       long   step;                         /* # of current testing    */
       long   time;                         /* time of current testing */
       int    init;                          /* 1 - linear, 0 - circle */
       double mean,mnsq;
       double *res,*sqres;                /* mean value & sq value   */
       double (*value)(struct Test *);    /* returning curr. value */
       };

struct Distribution {
        int flag;
        double xmin, xmax;
        double (*probability)(double);
        converter m;            
                     };

double jjfactor();

double endtoend(struct Test *);
double isclosed(struct Test *);
double angledist(struct Test *);

#define INITT   { {"endtoend",0,0,0,0,0,0,0,0,endtoend} }

#define NORES         1.79E+308
#define TIME          2.23E-308
                      
#define M_PI          3.1415926535897932384626433832795
#define M_1_16384     6.103515625E-5
#define M_PI_16384    1.91747598485705153714760948686E-4
#define M_15_LN2      10.397207708399179641258481821873
#define M_32768_SQRTPI  18487.36427369287801077760
#define M_SQRTPI_32768 5.40909988679661873563273E-5
#define randnorm(d)    sqrt(2*d*(M_15_LN2-log(rand()+1.)))*sin(M_PI_16384*rand())

char* readstr(char*, char*);
double distance(int, int);
void chain4view();
double tottwist(void);
double totenergy(void);


struct Data { double alpha,beta,gamma,ddac,ddb; 
			  double *dac, *db;
			  triplet *xh;
            } set;

triplet *xx;
vector *rr;
double *tw,*erg; 
int NS,NT;
struct Test test[]= INITT;
int step=0;
double rdata[3];
double restrict[3];
double condit[3];
double current[3];
double wolume;

double dotprod(vector a, vector b)
{ return a[0]*b[0]+a[1]*b[1]+a[2]*b[2]; }

double* crossprod(vector res, vector a, vector b)
{  res[0]=a[1]*b[2] - a[2]*b[1];
   res[1]=a[2]*b[0] - a[0]*b[2];
   res[2]=a[0]*b[1] - a[1]*b[0];
   return res;
}

/* The next function is approximate and good for normalization of a vector
   with the length close to 1 */
void norm(vector a)
{ double l;
   l=1.5-0.5*dotprod(a,a);
   a[0]*=l; a[1]*=l; a[2]*=l;
}

void out(triplet x)
{  printf(" (%f %f %f) %f\n", x[0][0],x[0][1],x[0][2], dotprod(x[0],x[0])); 
   printf(" (%f %f %f) %f\n", x[1][0],x[1][1],x[1][2], dotprod(x[1],x[1])); 
   printf(" (%f %f %f) %f\n", x[2][0],x[2][1],x[2][2], dotprod(x[2],x[2])); 
   printf("      %f %f %f\n", dotprod(x[0],x[1]), dotprod(x[1],x[2]), dotprod(x[2],x[0])); 
   getchar();
}

/* t=trans(x)*/
void setmatrix(matrix t, triplet x)
{  t[0][0]=x[0][0]; t[0][1]=x[1][0]; t[0][2]=x[2][0];
   t[1][0]=x[0][1]; t[1][1]=x[1][1]; t[1][2]=x[2][1];
   t[2][0]=x[0][2]; t[2][1]=x[1][2]; t[2][2]=x[2][2];
}

void settriplet(triplet x, double alpha, double beta, double gamma)
{ double ca,sa, cb,sb, cc,sc;

   ca=cos(alpha); sa=sin(alpha);
   cb=cos(beta);  sb=sin(beta);
   cc=cos(gamma); sc=sin(gamma);

   x[0][0]= ca*cb*cc-sa*sc;
   x[0][1]= ca*cb*sc+sa*cc;
   x[0][2]=-ca*sb;

   x[1][0]=-sa*cb*cc-ca*sc;
   x[1][1]=-sa*cb*sc+ca*cc;
   x[1][2]= sa*sb;

   x[2][0]= sb*cc;
   x[2][1]= sb*sc;
   x[2][2]= cb;
}

void rotmatrix(matrix t, vector x, double alpha)
{  double co, si, oc, x0, x1, x2;

   co=cos(alpha);   si=sin(alpha);

   oc=1.-co;
   x0=x[0]; x1=x[1]; x2=x[2];

    t[0][0]=x0*x0*oc+co; 
    t[1][1]=x1*x1*oc+co; 
    t[2][2]=x2*x2*oc+co;
    t[0][1]=t[1][0]=x0*x1*oc; 
    t[0][2]=t[2][0]=x0*x2*oc; 
    t[1][2]=t[2][1]=x1*x2*oc; 
    x0*=si; t[2][1]+=x0; t[1][2]-=x0; 
    x1*=si; t[0][2]+=x1; t[2][0]-=x1; 
    x2*=si; t[1][0]+=x2; t[0][1]-=x2;
}

void randvector(vector z)
{ double cb,sb,cc,sc;
   cb=((double)(rand()-16384))/16384;
   sb=sqrt(1.-cb*cb);
   cc=M_PI_16384*rand();
   sc=sin(cc);   cc=cos(cc);
   z[0]= sb*cc;
   z[1]= sb*sc;
   z[2]= cb;

/* or, alternatively, 
  int a,b,c,ll
  double l;
   do {
   a=(rand()-16384);
   b=(rand()-16384);
   c=(rand()-16384);
   ll=a*a+b*b+c*c;
   }  while(ll>805306368 || ll==0)
   l=1./sqrt(ll);
   z[0]=a*l;
   z[1]=b*l;
   z[2]=c*l;
*/

}

void matrixvector(vector x, matrix t, vector y)
{  x[0]=dotprod(t[0],y);
   x[1]=dotprod(t[1],y);
   x[2]=dotprod(t[2],y);
}

/* The next function also orthonormalizes the resulting triplet.
   It's not a simple rotation of each vector in the triplet!     */
void rotation(triplet x, matrix t, triplet x0)
{ double c1;
   matrixvector(x[2],t,x0[2]);
       norm(x[2]);
   matrixvector(x[1],t,x0[1]);
       c1=dotprod(x[1],x[2]);
       x[1][0]-=c1*x[2][0];
       x[1][1]-=c1*x[2][1];
       x[1][2]-=c1*x[2][2];
       norm(x[1]);
   crossprod(x[0],x[1],x[2]);   
} 

void casttriplet(triplet x, triplet x0)
{  x[0][0]=x0[0][0]; x[0][1]=x0[0][1]; x[0][2]=x0[0][2];
   x[1][0]=x0[1][0]; x[1][1]=x0[1][1]; x[1][2]=x0[1][2];
   x[2][0]=x0[2][0]; x[2][1]=x0[2][1]; x[2][2]=x0[2][2];
}

/* The next function is a relic after having been merged into 'rotation'.
   It corrects orthonormal triplet. */
void correction(triplet a)
{ double c1; 
   norm(a[2]);
   c1=dotprod(a[1],a[2]);
   a[1][0]-=c1*a[2][0];
   a[1][1]-=c1*a[2][1];
   a[1][2]-=c1*a[2][2];
   norm(a[1]);
   crossprod(a[0],a[1],a[2]);
}

void cmptriplet(double *bend, double *twist, triplet x,  triplet x0)
{ double cosb,pcosac,psinac;
   cosb=dotprod(x[2],x0[2]);
     *bend=acos(cosb);
   pcosac=(dotprod(x[0],x0[0])+dotprod(x[1],x0[1]));
   psinac=(dotprod(x[0],x0[1])-dotprod(x[1],x0[0]));
     *twist=atan2(psinac,pcosac);
}


double energy(triplet x, triplet x0, triplet xs, double dac, double db) 
{  triplet xe;
   matrix t;
   double euler_ac, euler_b;
   
   setmatrix(t,x0);

   rotation(xe, t, xs);

   cmptriplet(&euler_b, &euler_ac, x, xe);
   
   return 0.5*(euler_b*euler_b/db + euler_ac*euler_ac/dac);
}

double energyadd(vector a)
{ return sqrt(dotprod(a,a)); 
}

int notallowed(int k1, int k2, matrix rot)
{ double erg1, erg2, bend, twist, sqdist;
  triplet xs, xe;
  vector e1,e2;

   erg1=erg2=0.;
   sqdist=current[2];

   rotation(xs,rot,xx[k1]);
   if(k1!=0)  erg1=energy(xs,xx[k1-1], set.xh[k1-1],set.dac[k1-1],set.db[k1-1]);
   rotation(xe,rot,xx[k2-1]);
   if(k2!=NS+1) erg2=energy(xx[k2],xe, set.xh[k2-1],set.dac[k2-1],set.db[k2-1]);

   if( exp(erg[k1]+erg[k2]-erg1-erg2) < rand()/32768.) return 1; 

   if(k1==0 || k2==NS+1)
   {   e1[0]=rr[k2][0]-rr[k1][0];
       e1[1]=rr[k2][1]-rr[k1][1];
       e1[2]=rr[k2][2]-rr[k1][2];
       matrixvector(e2,rot,e1);
       e1[0]=rr[NS][0]+e2[0]-e1[0]-rr[0][0];
       e1[1]=rr[NS][1]+e2[1]-e1[1]-rr[0][1];
       e1[2]=rr[NS][2]+e2[2]-e1[2]-rr[0][2];
       sqdist=dotprod(e1,e1);
   }
   
   if( condit[2] < sqdist ) return 1;

   cmptriplet(&bend, &twist, (k1==0)?xs:xx[0], (k2==NS+1)?xe:xx[NS]);
   twist=fabs(twist);

   if(bend > condit[1] || twist > condit[0]) return 1;

    current[2]=sqdist;
	current[1]=bend;
	current[0]=twist;
    erg[k1]=erg1;
    erg[k2]=erg2;
   return 0;
}

void rotcontour(int start, int end, matrix t)
{ int i;
  triplet x;

  for(i=start; i<end; i++)
  {    rotation(x,t,xx[i]);
       casttriplet(xx[i],x);
  }

  if(start!=0)
  for(i=start; i<end || i<NS; i++)
  {    rr[i+1][0]=rr[i][0]+xx[i][2][0];
       rr[i+1][1]=rr[i][1]+xx[i][2][1];
       rr[i+1][2]=rr[i][2]+xx[i][2][2];      
  }
  else
  for(i=end-1; i>=0; i--)
  {    rr[i][0]=rr[i+1][0]-xx[i][2][0];
       rr[i][1]=rr[i+1][1]-xx[i][2][1];
       rr[i][2]=rr[i+1][2]-xx[i][2][2];      
  }
	rr[NS+1][0]=rr[NS][0];
    rr[NS+1][1]=rr[NS][1];
    rr[NS+1][2]=rr[NS][2];
}

void move()
{ int k1,k2,*count;
  vector a;
  matrix rot;
  double l,*ampl;
  static double ampl1=0.1,ampl2=0.1;
  static int   count1=0, count2=0;

   k1=rand()%NS;
   k2=rand()%NS + 1;

  if(k2>k1) 
   { a[0]=rr[k2][0] - rr[k1][0];
	 a[1]=rr[k2][1] - rr[k1][1];
	 a[2]=rr[k2][2] - rr[k1][2];
	 l=1/sqrt(dotprod(a,a));
	 a[0]*=l; a[1]*=l; a[2]*=l;
	 ampl=&ampl2;	 count=&count2;
   }
  else
   { k1=rand()%NS + 1;
	 if(2*k1>NS) k2=NS+1;
	 else 
	 { k2=k1; k1=0;	 }

	 randvector(a);
     ampl=&ampl1;	 count=&count1;
   }
             
   rotmatrix(rot,a,(*ampl)*(rand()-16384)/16384.);

   if(notallowed(k1,k2,rot))
   {   if(--*count <-30) {*ampl*=0.9; *count=0; }
	   return;
   }
       if(++*count > 30) {*ampl*=1.1; *count=0; }
	   
   rotcontour(k1,k2,rot);
}

void generate()
{ int i;
  double angle,v,mv,dv;
  matrix t;

			condit[0]=2*(restrict[0]=rdata[0]);
			condit[1]=restrict[1]=rdata[1];
			condit[2]=restrict[2]=rdata[2];

    rr[0][0]=rr[0][1]=rr[0][2]=0;
    xx[0][0][0]=1; xx[0][0][1]=0; xx[0][0][2]=0;
    xx[0][1][0]=0; xx[0][1][1]=1; xx[0][1][2]=0;
    xx[0][2][0]=0; xx[0][2][1]=0; xx[0][2][2]=1;
    rr[0][0]=1; rr[0][1]=rr[0][2]=0;
	erg[0]=0;
	angle=2*M_PI/NS;
    rotmatrix(t,xx[0][0],angle);    

	for(i=0; i<NS;i++)    
	{   rotation(xx[i+1],t,xx[i]);
	    erg[i+1]=energy(xx[i+1],xx[i],set.xh[i],set.dac[i],set.db[i]);
	    rr[i+1][0]=rr[i][0]+xx[i][2][0];
        rr[i+1][1]=rr[i][1]+xx[i][2][1];
        rr[i+1][2]=rr[i][2]+xx[i][2][2];
	}
	rr[NS+1][0]=rr[NS][0];
	rr[NS+1][1]=rr[NS][1];
	rr[NS+1][2]=rr[NS][2];
	erg[NS+1]=0;
/*
	for(i=1; i<=NS; i++)    
	{   rotmatrix(t,xx[i][2],2*M_PI*i/NS);    
        rotation(x,t,xx[i]);
		casttriplet(xx[i],x);
	    erg[i]=energy(xx[i],xx[i-1],set.xh[i-1],set.dac[i-1],set.db[i-1]);
	}
*/
    cmptriplet(&current[1],&current[0],xx[0],xx[NS]);
    current[0]=fabs(current[0]);
	current[2]=distance(0,NS);

	for(mv=0,dv=0,i=0;i<400000;i++) 
	{  move();
	   v=tottwist();
       mv += (   v-mv )/(i+1);
       dv += ( v*v-dv )/(i+1);
    
	}

/*	printf("writhe %e %e\n",mv,sqrt(dv-mv*mv));*/
}

void allocmemory()
{ xx=(triplet*)malloc(sizeof(triplet)*(NS+1));
  rr=(vector*)malloc(sizeof(vector)*(NS+2));
  erg=(double*)malloc(sizeof(double)*(NS+2));
  set.xh=(triplet*)malloc(sizeof(triplet)*(NS+1));
  set.dac=(double*)malloc(sizeof(double)*(NS+1));
  set.db=(double*)malloc(sizeof(double)*(NS+1));
  if(xx==NULL || rr==NULL || erg==NULL || set.xh==NULL || set.dac==NULL || set.db==0)
  {   printf("Memory problem\n"); exit(0);
  }
}

void setparameters()
{  int i=0;
	 settriplet(set.xh[i], set.alpha, set.beta, set.gamma);
         set.dac[i]=set.ddac; 
         set.db[i]=set.ddb;
         i++;
    	  for(;i<=NS;i++)
	  {   
		   set.dac[i]=set.ddac; 
		   set.db[i]=set.ddb;
	           casttriplet(set.xh[i], set.xh[i-1]);
	  }
}


void readdata(char *fname)
{  FILE *file;
   int nf,bp;
   double bnd,p,rpt,C;
   char c[80], *tail;

   if((file=fopen(fname,"rt"))==NULL) 
   {   printf("Can't open %s\n", fname); exit(0);
   }

   while(fgets(c,sizeof(c),file)==c)
   {
       if((tail=readstr(c,"Number of fragments"))) nf=atoi(tail);
       if((tail=readstr(c,"Fragment length, bp"))) bp=atoi(tail);
       if((tail=readstr(c,"Static bend"))) bnd=atof(tail);
       if((tail=readstr(c,"DNA Kuhn length, bp"))) p=atof(tail);
       if((tail=readstr(c,"Helical repeat, bp"))) rpt=atof(tail);
       if((tail=readstr(c,"Torsional rigidity, ergcm"))) C=atof(tail);
		
   }
   fclose(file);

   NS=nf*bp/rpt;
   set.alpha=0.;
   set.gamma=2*M_PI*(bp/rpt-(int)(bp/rpt+0.5))*nf/NS;
   set.beta=bnd*nf/NS;
   set.ddac=(3.4e-8*bp*nf/NS)*1.38e-16*293/C;
   p=2/(1+p*NS/(bp*nf)); /* reusing p */
   set.ddb = p + 2*p*p/3 + 26*p*p*p/45 + 526*p*p*p*p/945;
   allocmemory();
   setparameters();

   rdata[0]=0.1; rdata[1]=0.1; rdata[2]=0.1;
   wolume=1.5/( rdata[0]*(1-cos(rdata[1]))*sqrt(rdata[2])*rdata[2]);

   p=3.4e-9*bp*nf/NS;      /* reusing p again */
   wolume/=6.022e23*p*p*p;    /* convertion to moles per liter */

}

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

void tests()
{struct Test *tst;
 double v;
   for(tst=test; tst<test+NT; tst++)
     if(tst->step && step%tst->fr==0)
       { v=(*(tst->value))(tst);
	    if(v==NORES) continue;
	 tst->mean += (  v-tst->mean )/tst->step;
	 tst->mnsq += (v*v-tst->mnsq )/tst->step;
         if(tst->step%1024==0)
         printf("%10s %5d n=%-11ld x=%-12.4e Mx=%-12.4e Sx=%-12.4e\n",\
         tst->name,tst->fr,tst->step,v,tst->mean,sqrt(tst->mnsq - tst->mean*tst->mean));
         tst->step++;
         fflush(stdout);
       }
}

void simulation(double *mmv, double *eev)
{ double v, mv=0., dv=0.;
  int i=0;
  do
  {
    generate();
    do
    {  step++;
       move();
       v=jjfactor();
    } while(v==NORES);
    i++;
    mv += (   v-mv )/i;
    dv += ( v*v-dv )/i;
  } while(i<9);
  *mmv=mv; 
  *eev=sqrt(dv-mv*mv);
}


int main(int argc,char **argv)
{ double mv,ev;
	 srand((unsigned)time(NULL));

     NT=sizeof(test)/sizeof(test[0]);
     
	 if(argc>1)  readdata(argv[1]); 
	 else        exit(0);  
     
     simulation(&mv,&ev);

     printf("%e %e\n",mv,ev);

   return 1;
}

double endtoend(struct Test *test)
{  return distance(0,NS);
}

double isclosed(struct Test *test)
{ return ( endtoend(test)< 1.);
}

double jjfactor()
{  int freq=400000;
   double pc;

   static int dt=0, count=0, dtall=0;
   static double v=1.;
  
   if(  restrict[count] > current[count] )   dt++;
   
   dtall++;
   if(dtall%freq==0)
   {  pc=(double)dt/freq;
      v*=pc;
      /*printf("%.2e (%f) ",pc, tottwist()); */
	  restrict[count]=condit[count]; 
	  if(dt==freq)  
	  {   restrict[count]=condit[count]=999;
		  count++; 
		  if(count==3) 
		  {   pc=v*wolume;
			  /*printf(" J-factor: %e %e  Done!\n", v, pc);*/
			  fflush(stdout);
			  dt=0; count=0; dtall=0; v=1.;                           
		      return pc;
		  }
	  }
      condit[count]*=2;    
      dt=0;
      dtall=0;
   }

   return NORES;
}


double distance(int i, int j)
{ vector h;
     h[0]=rr[j][0]-rr[i][0];
     h[1]=rr[j][1]-rr[i][1];
     h[2]=rr[j][2]-rr[i][2];
  return dotprod(h,h);
}

double tottwist()
{ int i;
  double twist,bend,total=0;

  for(i=0;i<NS;i++)
  {
    cmptriplet(&bend, &twist, xx[i+1], xx[i]);
	total+=twist;
  }
    cmptriplet(&bend, &twist, xx[NS], xx[0]);

  return total-twist;
}

double totenergy()
{ int i;
  double toten=0;

  for(i=1;i<=NS;i++) toten+=erg[i];

  return toten;
}