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

using namespace std;

/* This program, "ideal.cc", calculates the velocity, for frictionless  
   sliding of the dominoes, of the fully developed domino effect. 
   It is faster and more accurate than "friction.cc" with mu=0.0.
*/ 

enum {Nt = 40};//total number of separations
double  sep, dth, H, I, J, H0, I0, J0, P, oma;
const double d = 0.17863; //aspect ratio of the dominoes 
int Ni, Nc; //Ni=# if points in theta integration, Nc=# of steps in recursion

void coefficients(double x);   //computes H,I,J for x=sin(theta)
double iomega(double x);   //yields invers omega(x=sin(theta))
double time();   //integrates over theta

int main() {
  double alpha, beta;
  cout << setprecision(10);
  cout << setiosflags(ios :: showpoint);
  cout << setiosflags(ios :: fixed);  
  for(int i=1; i<Nt; i++) {
    Ni = 20*i; Nc=10*(Nt-i+5); //estimates for integration mesh  and number of collisions
    sep = i;  sep /= Nt;  dth = asin(sep)/Ni; //sets separation and increment theta
    coefficients(0);  H0=H;  I0=I;  J0=J; //initial values are set
    //cout << H0 << "  " << I0 << "  " << J0 << endl;
    P=(sep-d*sqrt(sep*sep+2.0*d*sep))/(sep+d); //fuel
    //For the theory of Banks the next line should be substituted
    //P = 1- sqrt(1-sep*sep) -d*sep;//Banks
    if(P<0) {cout << sep << " no domino effect possible" << endl; continue;}
    //make a choice out of the following three possobilities (default: present theory):
    alpha=I0; beta=I0-1;//present theory
    //alpha=J0; beta=J0-1;//Shaw
    //alpha=1; beta=sqrt(1-sep*sep); //Banks
    oma=beta*sqrt(P/(alpha*alpha*(I0-1)-beta*beta*I0)); //initial omega
    //cout << P << " init oma = " << oma << " fin oma = " << omf << "   ";
    double tim=time(); //integrates the time for the present separation
    double vas = (sep+d)*sqrt(3/(1+d*d))/tim; //asymptotic velocity
    cout << sep << '\t' << vas <<  endl;
  }
  return 0;
} 

void coefficients(double x) {//potential and effective moments
  double  X, Y, Z, y, xn,yn, w,wn;
  y= sqrt(1 - x*x); w = 1; I=J=1.0;
  H = y + d * x;  
  for(int j=1; j<Nc; j++) {
    X = (sep+d)* y - d; Y = sqrt(1 - X*X);
    xn = Y *x + X *y;    yn = -X *x + Y *y;
    H += yn + d * xn;
    Z = 1.0-(sep+d)*x /Y;
    wn = Z*w; I +=wn*wn; J +=wn;
    x=xn; y=yn; w=wn;
  }  
}

double iomega(double x) {
  coefficients(x);
  return  sqrt(I/(I0*oma*oma+H0-H));
}
 

double time() {//time per segment
  double sum = 1/oma; 
  for(int i=1,s=1; i<=Ni; i++,s=-s) {
    double the = i*dth; double x = sin(the);
    sum += (i<Ni ? (s+3) :1)*iomega(x); //implementing Simpson's rule
  }
  //cout << "final oma= " << 1/iomega(sep) << "   ";
  return  sum*dth/3;
}