#include <iostream>
#include <cmath>
#include <iomanip>

using namespace std;
/* This program, "friction.cc" calculating the velocity of the domino effect for 
   dominoes sliding with friction coeffcient mu over each other while rotating
*/
enum {Nt = 40}; //number of separations s/h.
//aspect parameters 
const double  d=0.17863, mu=0.1; // d-> d/h, mu=friction, 
double sep, dth, duration; //sep=s/h, dth=dtheta, duration=time
double A,B,C; //coefficients of differential equation for omega(theta)
int Ni, Nc;

void coefficients(int n, double x);//computes A,B,C
double derivative(int n , int i, double oma);//calculates d omega/ d theta
double collision(int n, double oma);//collision
double time(int n, double oma);//integrates time

int main() {
  int i, n;
  double thetac, oma, omp, vas;
  cout << setprecision(10);
  cout << setiosflags(ios :: showpoint);
  cout << setiosflags(ios :: fixed);
  for(i=1; i<Nt; i++) {
    Ni = 20*i;   //estimate for integration mesh  
    Nc = 8*(Nt-i+5); //estimate of number of collisions for convergence
    oma = 1.0;
    sep = i; sep /= Nt;
    thetac = asin(sep);  dth = thetac/Ni; 
    for(n=1; n<Nc; n++) {  
      omp = time(n,oma); //time of interval n
      oma = collision(n,omp);  //new starting rotation velocity from collision
    }   
    vas = (sep+d)/duration;  vas *= sqrt(3/(1+d*d)); //asymptotic velocity
    cout << sep << '\t' << vas << endl;
  }
  return 0;
} 

void coefficients(int n, double x) {
  int j;
  double  X, Y, y, xn, yn, w, wn, r, v, T, ici, ai, bi, ti, zi;
  y = sqrt(1 - x*x); 
  r= w = 1.0; v = 0.0; //initial values w=psi', v=psi"
  A = 1.0; B = 0.0; C = 0.5*(x - d*y); //initial values coefficients
  for(j=1; j<n; j++) {
    X = (sep+d)* y - d;  Y = sqrt(1 - X*X);
    xn= Y*x + X*y;  yn = -X*x + Y*y; //next angle
    T = 0.5*(xn - d*yn);  //new torque
    ici = 1/Y; ti = X*ici;
    zi = 1 - (sep+d)*x*ici; 
    ai = zi - mu*d*ici; bi = 1 + mu*ti; 
    if(ai<0) break;  else r *= ai/bi; //if ai<0 end of recursion
    if(zi<0) {wn = 0; v = 0;} 
    else {wn = w*zi; v *= zi;  v += ti*(wn-w)*(wn-w)-(sep+d)*y*w*w*ici;}
    A += r*wn; B += r*v; C += r*T; //formation of coefficients
    x=xn; y=yn; w=wn;  //input for next step in loop
  }  
}

double derivative(int n, int i, double oma) {
  int j;
  coefficients(n,sin((i+0.5)*dth)); //coefficients A,B,C at theta=(i+0.5)*dth
  return (C/oma - B*oma)/(A+0.5*dth*(C/(oma*oma)+B)); //3rd order scheme
}
 
double collision(int n, double oma) {
  int j;
  coefficients(n+1,0); //coefficients at theta=0, computes A(0)
  return oma*(A-1)/A; //new starting value
  
}

double time(int n, double oma) {
  double sum = 1/oma; //start interval
  for(int i=0,s=1; i<=Ni; ++i,s=-s)  {
    oma += dth*derivative(n,i,oma); //new omega
    sum += (i<Ni ? s+3: 1)/oma; //implementing Simpson weights
  }
  duration = sum*dth/3; //duration interval n
  return oma;
}