/*
This code simulates the reptation model (eq.10) in Cates1987 without breakage and fusion
1. The reptation of a polymer within a tube
=> How does it change the viscosity with L0?
*/

#include<iostream>
#include<stdlib.h>
#include<cmath>
#include<sstream>
#include<fstream>
#include<string>
#include<iomanip>
#include<vector>
using namespace std;

/* CONSTANTS */
const int REPMAX=300;
const double dL=0.01; //
const double dt=0.01; //integration time
const double tau=1.0; //unit of time (100dt)
double L0=1.0; //L0 sets the units of length (100dL)
const double D0=0.01; //monomer diffusion. In units of L0^2/dt
                        //if L0=1; set D0=(1/dt) to get L0^2 monomer diffusion every
                        //tau=(1/dt) steps. [D0]=0.01[L0^2/dt]=1[L0^2/tau]
const int bin=10;//for histogram. Time units
const int Tmaxmax=10000000; // this is bc it wants something not variable for arrays

/* ARRAYS */
double pos_left[Tmaxmax];
double pos_right[Tmaxmax];
int RelaxTime[REPMAX];
int Histo[Tmaxmax];
int Survival[Tmaxmax];
int CountBreaks[REPMAX];
int CountFusions[REPMAX];

/* FUNCTIONS */
double Diff(double L);

/* MAIN */
int main(int argc, char* argv[]){
srand(time(NULL));
//cout << "Trep=" << Trep <<endl;

double L0=atof(argv[1]);
double Trep=pow(L0,3.0)/D0; //reptation time (units dt)
int Tmax=int(100*Trep); //true maximum time (units dt)

double xi=0;
double c1=(1.0/Trep)*(1/xi);//1/Trep*1/Xi
double Tbreak=1/(c1*L0);
//double xi=Tbreak/Trep;

/*Initialise arrays*/
for(int t=0;t<Tmax;t++) pos_left[t]=pos_right[t]=Survival[t]=0;
for(int i=0;i<REPMAX;i++)RelaxTime[i]=Tmax;

/* CREATE/OPEN FILES */
stringstream namepos;
namepos << "Position_"<<L0<<".dat";
ofstream writepos;
writepos.open(namepos.str().c_str());
//
stringstream namefi;
namefi << "Survival_"<<L0<<".dat";
ofstream writeSu;
writeSu.open(namefi.str().c_str());
//
stringstream namefieta;
namefieta << "eta_"<<L0<<".dat";
ofstream writeETA;
writeETA.open(namefieta.str().c_str());
//
stringstream namefitopo;
namefitopo << "topo_"<<L0<<".dat";
ofstream writetopo;
writetopo.open(namefitopo.str().c_str());
//

/*Start loop over attempts*/
for(int i=0; i<REPMAX; i++){
cout << "============= " <<endl;
cout << "== NEW REP == " <<endl;
cout << "============= " <<endl;

/* Define initial condition with Length exponentially distributed ~exp(-l/l0) */
//left end
pos_left[0]=-L0*log(rand()*1.0/RAND_MAX); //sample from exponential distro
//right end
pos_right[0]=-L0*log(rand()*1.0/RAND_MAX); //sampled from expoential ditro

/* Tube segment starts at 0 */
writepos << 0 << " " << -pos_left[0] <<" " << pos_right[0]<< " "<< i << endl;

double b1,b2,L;

/* =============== */
/*    START TIME   */
/* =============== */
/* Start time. The particle is at 0.*/
/* Pos_left is actually to negative value*/
for(int t=1;t<=Tmax; t++){

L=pos_right[t-1]+pos_left[t-1];
b1=pos_left[t-1];
b2=pos_right[t-1];

/* At each timestep there is a chance of breakage or fusion */
/* BREAKAGE at rate c1 per unit time per unit length */
/*
for(int l=0;l<int(b1/dL);l++){
//cout << "attemp break @ " << l << endl;
double pbreak_left=rand()*1.0/RAND_MAX;
    if(pbreak_left<c1*dL*dt){
    //cout << "pbreak_left " << pbreak_left << " " << c1*dL*dt <<endl;
    cout << "LEFT BREAK " << b1<<" -> " << l*dL << endl; //cin.get();
    b1=l*dL;
    CountBreaks[i]++;
    break;
    }
}

for(int l=0;l<int(b2/dL);l++){
double pbreak_right=rand()*1.0/RAND_MAX;
    if(pbreak_right<c1*dL*dt){
    //cout << "pbreak_right " << pbreak_right << " " << c1*dL*dt <<endl;
    cout << "RIGHT BREAK " << b2<<" -> " << l*dL << endl; //cin.get();
    b2=l*dL;
    CountBreaks[i]++;
    break;
    }
}
*/

/* FUSION at rate c2 per unit time per unit length */
/* in reality it becomes c2*N(l)=c2*c1/c2*exp(-l/L0)=c1*exp(-l/L0) */
/* so propto the concentration of chains of length l */
//double l_left=-L0*log(rand()*1.0/RAND_MAX);
/*
for(int l_left=1;l_left<1000;l_left++){
double pfus_left=rand()*1.0/RAND_MAX;
if(pfus_left<c1*dL*dt*exp(-l_left/L0)){
    cout << "LEFT FUSION " << b1<<" -> " << b1+l_left << endl;
    b1=b1+l_left;
    CountFusions[i]++;
    break;
}
}
//double l_right=-L0*log(rand()*1.0/RAND_MAX);
for(int l_right=1;l_right<1000;l_right++){
double pfus_right=rand()*1.0/RAND_MAX;
if(pfus_right<c1*dL*dt*exp(-l_right/L0)){
    cout << "RIGHT FUSION " << b2<<" -> " << b2+l_right << endl;
    b2=b2+l_right;
    CountFusions[i]++;
    break;
}
}
*/

/* NEW LENGTH */
L=b2+b1;

/* UPDATE */
//pos_left[t]=b1;
//pos_right[t]=b2;

/* DIFFUSE THE TWO ENDS */
/* USE b1 and b2 otherwise missing the breaks */
/* CAREFUL! They are not independent otherwise the chain changes length! */
int dir=(rand()%2)*2-1; //+1 or -1
pos_left[t]=b1+dir*sqrt(2*Diff(L)*dt);
pos_right[t]=b2-dir*sqrt(2*Diff(L)*dt);

cout << t*dt << " " << -pos_left[t] << " " << pos_right[t] << " "<< i << endl;
writepos << t*dt << " " << -pos_left[t] << " " << pos_right[t] << " "<< i << endl;

/* Survival function */
if(pos_left[t]>0 && pos_right[t]>0)Survival[t]++;

/* if one of them hit 0 break */
if(pos_left[t]<0 || pos_right[t]<0){
RelaxTime[i]=t;
break;
}

}//close loop over time

writepos << endl;
writepos << endl;

}//close loop over attempts

/* SURVIVAL */
writeSu <<"#t P_survived(t) " << endl;
for(int i=0;i<=Tmax;i++) writeSu << i*dt << " " << Survival[i]*1.0/REPMAX << endl;

/* VISCOSITY*/
/* this is the integral of the survival function */
double eta=0;
for(int i=0;i<Tmax;i++) eta+=Survival[i]*1.0/REPMAX*dt;
cout << L0 << " " << xi << " " <<  eta <<endl;
writeETA << L0 << " " << xi << " " <<  eta <<endl;


double meanb=0;
double meanfus=0;
for(int i=0;i<REPMAX;i++){
writetopo << i << " " << CountBreaks[i] << " " << CountFusions[i] << " " << RelaxTime[i] <<endl;
meanb+=CountBreaks[i]*1.0/RelaxTime[i];
meanfus+=CountFusions[i]*1.0/RelaxTime[i];
}
writetopo<<endl;
writetopo<<endl;
writetopo << "Xi c_1 <breaks/time> <fusion/time> " <<endl;
writetopo << xi << " " << c1 <<" " << meanb/REPMAX << " " << meanfus/REPMAX <<endl;

return 0;
}


/* curivlinear diffusion coeff and function of L */
double Diff(double L){
return D0/L; //L is in units of L0
}