#include<iostream>
//#include<fstream>
#include<iomanip>
#include<cstdlib>
#include<stdlib.h>

using namespace std;

extern  void spline_init(double*, double*, int, double*);
extern  void spline_interpolation(double*, double*, double*, int, double, double*, double* );

extern "C"{void dsyev_( char* JOBZ, char* UPLO, int* N, double* A, int* LDA, double* W, double* WORK, int* LWORK, int* INFO );}


class QC_MOM 
{

  private: static const int Ndvr = 128;
           static const int LLeg = 30;

 
           double* Xval; double* Yval; double* MOMS; int Npoints;
	   double* VEIGEN; double* Matrix;
           int Ntraj; 
           double pr; double pt; double pp;
	   double X,Y,R, THETA , PHI, BETA;  
         
           double mu;
           void diagonalisation(void);
           int INFO2;
           double* VMIN; double* REQ;
           double t_kin(int i, int j, double dx, double mass);
	   void prop();
//	   double* MD_Y; double* MD_X; double* MOMS_MD; int N_MD; double MD_TIME;

	   double* RMIN; double* RMAX;
           int QVIB; int QL; int QM; //QUAMTUM NUMBERS
           constexpr static double L_surf = 2.89769791; // MS.RPBE X vector length for 1 unit cell, with 170 K lattice expansion 
           //FOR THE MOLECULAR BEAM
           double BEAM_FUNC(double x, double a, double x0);
           void set_vbeam(double x){VBEAM = x;};
           double ALPHA; double VBEAM;
           double NORM_BEAM;
           void MAKE_STREAM_SETUP();
           double INT_A; double INT_B; //Inteval for the stream  v_stream_min and v_stream_max
           double DX_STREAM; //INTEGRATION INTERVAL FOR THE BEAM        
           int N_STREAM;     //NUMBER OF INTEGRATION INTERVALS FOR THE BEAM
           double RAND_VBEAM; //randomly selected volecity of the beam
           double* VEIGEN_ql;
           double TSTREAM;
           void SEARCH_MIN_POT(double dx, int Npoints, double* OUT_REQ, double* OUT_VMIN, int L);
           double PARTITION_FUNCTION; double kbT; double KbTrot;
           void DETERMINE_R(int qvib, int lrot, double* R_out, double rand_nr);
           int* QVIB_LROT_SUM;
           void  DETERMINE_THETA_PHI(int lrot, int mJ);
           void DETERMINE_RMIN_RMAX(double* xvals, double* yvals, double* moms, int qvib, int lrot);
           int BEAM_SWITCH;
           int STATE_SWITCH;
           int* QVIBLM;
           double KbT;
           void add_qvib(int qvib, int lrot){QVIB_LROT_SUM[qvib*LLeg + lrot] += 1;};
           void GET_EROT(double x, double L, double* erot, double* derot);
  public: 
      void INITPOS(char* file, double mass, int TON, double alpha, double vstream, double Tstream, int state, int qvib, int ql, int qm);
      int check_status(){return(INFO2);}; 
      void getposr(double* x, double* px, int i);
      void init_pos_mom();
      double get_eigenval(int i){return(VEIGEN[i]);};

      double get_rmax(int i ){cout << "RMAX[" << i << "] = " << RMAX[i] << "\n\n"; return(RMAX[i]);};
      double get_rmin(int i ){return(RMIN[i]);};
      double get_vstream(){return(RAND_VBEAM);};

      double GET_EXPECT_INT();
      double GET_EXPECT_INT_EXPECTED();
      double GET_POT(double R);

      int get_qm(){return(QM);};
      int get_qvib(){return(QVIB);};
      int get_ql(){return(QL);};

};
void QC_MOM::GET_EROT(double x, double L, double* erot, double* derot){
   *erot  = (double)(L*(L+1.0))/(2.*mu*x*x);
   *derot = -(double)(L*(L+1.0))/(mu*x*x*x);
}
double QC_MOM::GET_POT(double rin){
  double vpot = 0.; double dvdx = 0.;
  spline_interpolation(Xval, Yval, MOMS, Npoints, rin, &vpot, &dvdx );
  return(vpot);
}
double QC_MOM::GET_EXPECT_INT(){
  double y = 0.;
  for(int i = 0; i < 4; i++){
    for(int j = 0; j < LLeg; j++){
      y += QVIB_LROT_SUM[i*LLeg + j]*VEIGEN_ql[i*LLeg + j];
//      cout << i << " " << j << " " << QVIB_LROT_SUM[i*LLeg + j] << " " << VEIGEN_ql[i*LLeg + j] << "\n"; 
   }
  }
 
  return(y);
};
double QC_MOM::GET_EXPECT_INT_EXPECTED(){
   double sum = 0.; double fac = 0.; double fac2 = 0.;
   for(int i = 0; i < 4; i++){
     for(int j = 0; j < LLeg; j++){
       fac  = (VEIGEN_ql[i*LLeg] - VEIGEN_ql[0])/KbT;
       fac2 = (VEIGEN_ql[i*LLeg + j] - VEIGEN_ql[i*LLeg])/KbTrot;
       sum += (2.*j+1.)*VEIGEN_ql[i*LLeg + j] * exp(-fac) * exp(-fac2);     
     }
   }
   sum /= PARTITION_FUNCTION;
   return(sum);  
}

double QC_MOM::BEAM_FUNC(double x, double a, double x0){
  double z = (x-x0)*(x-x0)/(a*a);
  return(exp(-z)*x*x*x);
};
void QC_MOM::MAKE_STREAM_SETUP(){
  DX_STREAM = VBEAM*1e-04;
  INT_A = 0.; INT_B = 0.;
  int yes = 1; 
  double beam_min_val = 2e-07; // minimum beam value to determine the valid beam interval
  double ymax = 0.; double xmax = 0.;
  int count = 0; double v = 0.; double fac = 0.;
  int count2 = 0;
  //FIRST FIND MAXIMUM VALUE OF THE BEAM DISTRIBUTION

  while(yes == 1){
    v = DX_STREAM*count;
    fac = BEAM_FUNC(v, ALPHA, VBEAM) - BEAM_FUNC(v+DX_STREAM, ALPHA, VBEAM);
    if(fac > 0.){ yes = 0; ymax = BEAM_FUNC(v, ALPHA, VBEAM); xmax = v;} 
    count++;
  }
  count = 0; yes = 1;  NORM_BEAM = ymax*1.1;

  //DETERMINE THE INTERVAL RANGE
  while(yes == 1){
    v = DX_STREAM*count;
    fac = BEAM_FUNC(v, ALPHA, VBEAM)/NORM_BEAM;
    if(fac > beam_min_val  ){ INT_A = v; yes = 0;}
    count++;
  } 
  yes = 1;
   while(yes == 1){
    v = DX_STREAM*count;
    fac = BEAM_FUNC(v, ALPHA, VBEAM)/NORM_BEAM;
    if(fac < beam_min_val  ){ INT_B = v; yes = 0;}
    count2++; 
    count++;
  }

  N_STREAM = count2;

  ofstream out_stream; out_stream.open("mol_beam_info.dat");
  out_stream << "# vstream                            " << VBEAM << " m/s n";
  out_stream << "# alpha                              " << ALPHA << " m/s n";
  out_stream << "# BEAM INTERVAL:               INT_A " << INT_A << " INT_B " << INT_B << "\n";
  out_stream << "# BEAM VALUE AT INTERVAL END POINTS  " << BEAM_FUNC(INT_A, ALPHA, VBEAM)/NORM_BEAM  << "   " << BEAM_FUNC(INT_B, ALPHA, VBEAM)/NORM_BEAM << "\n";
  out_stream << "# BEAM INTERGRATION INCREMENT        " << DX_STREAM  << "\n";
  out_stream << "# NUMBER OF INTERGRATION INCREMENTS  " << N_STREAM   << "\n";
  out_stream << "# MAXIMUM VALUE OF BEAM DISTRIBUTION " << NORM_BEAM  << "\n";

  for(int i = 0; i < N_STREAM; i++){
    v = DX_STREAM*i + INT_A;
    out_stream << v << " " << BEAM_FUNC(v, ALPHA, VBEAM)/NORM_BEAM << "\n";
  }


};
void QC_MOM::diagonalisation(void ){
     char JOBZ = 'V';
     char UPLO = 'U';
     int  LDA  = Ndvr;
     int LWORK = 4*Ndvr; int INFO;
     double* WORK = new double[LWORK];
     int N = Ndvr;
     dsyev_(&JOBZ, &UPLO, &N, Matrix, &LDA, VEIGEN, WORK, &LWORK, &INFO );
     INFO2 = INFO;
//     for(int i = 0; i < N; i++) VEIGEN[i] *= 0.95555556;
//     cout << INFO << "\n";
     delete WORK;
};
double QC_MOM::t_kin(int i, int j, double dx, double mass){
       double val = 1./(2.*mass*dx*dx) * pow((-1), (i-j))   ; double y;
       if(i != j)  y = 2./ ((i-j)*(i-j)) - 2./ ((i+j)*(i+j)) ;   
       if(i == j)  y = M_PI*M_PI/3. - 1./(2*i*i) ;
       return(y*val);
};

void QC_MOM::SEARCH_MIN_POT(double dx, int Npoints, double* OUT_REQ, double* OUT_VMIN, int L){ 
  //SEARCH FOR LOWEST POTENTIAL ENERGY AND EQUILIBRIUM GEOMETRY
  double dt = 1e-02; double conv = 0.; double xold = Xval[0]; double xnew;
  int count = 0; double vpot; double dvdx; double Erot = 0.; double dErot = 0.;
  int yes = 1; double xsave; double vsave;
  do{
    spline_interpolation(Xval, Yval, MOMS, Npoints, xold, &vpot, &dvdx );
    GET_EROT(xold, L, &Erot, &dErot);
    vpot += Erot;
    dvdx += dErot;
    xnew  = xold - dt*dvdx;
    conv = fabs(dvdx);//fabs(xnew-xold);
    if(conv < 1e-10 && yes == 1){
      xsave = xold;
      vsave = vpot;
      yes = 0;
    }
    xold = xnew;
    count++;
  }while(yes == 1);
  *OUT_REQ = xold; *OUT_VMIN = vpot;
};

void QC_MOM::INITPOS(char* file, double mass, int TON, double alpha, double vstream, double Tstream, int state, int qvib, int ql, int qm){
//  time_t seed_time;
//  time(&seed_time);
//  srandom((int)seed_time);
 
  int NStates = LLeg*4; //number of ro-vibrational states l = 0 ... LLeg - 1 and v = 0,1,2
  Matrix = new double[Ndvr*Ndvr];
  VEIGEN = new double[Ndvr];
  VEIGEN_ql = new double[NStates]; //ro-vibrational eigenstates Ev,l wheer v = 0,1,2 and L = 0, ... LLeg - 1
  ALPHA = alpha; VBEAM = vstream;  MAKE_STREAM_SETUP(); 
  TSTREAM = Tstream; //beam temperature
  mu = mass;
  BEAM_SWITCH = TON;
  STATE_SWITCH = state;
  QVIBLM = new int[3] { qvib, ql, qm };
  RMIN = new double[NStates];
  RMAX = new double[NStates]; 
  VMIN = new double[LLeg];  
  REQ = new double[LLeg];
  QVIB_LROT_SUM = new int[NStates];
  for(int i = 0; i < NStates; i++) QVIB_LROT_SUM[i] = 0;

  ifstream inp; inp.open(file);
  inp >> Npoints; double dum;
  Xval = new double[Npoints];
  Yval = new double[Npoints];
  MOMS = new double[Npoints];
  for(int i = 0; i < Npoints; i++){
    inp >> Xval[i] >> dum;
    Xval[i] *= 1.8897261;
    //DETERMINE THE 
    Yval[i] = dum/27.211383; 
  }
  
  spline_init(Xval, Yval, Npoints, MOMS);
 

  double* V = new double[Ndvr];
  double dx = (Xval[Npoints-1] - Xval[0])/(Ndvr);
  double dvdx, vpot, Erot, dErot;


  //SEARCH FOR LOWEST POTENTIAL ENERGY AND EQUILIBRIUM GEOMETRY
  SEARCH_MIN_POT(dx, Npoints, &REQ[0], &VMIN[0], 0);
  for(int i = 0; i < Npoints; i++) Yval[i] = (Yval[i] - VMIN[0]);
  spline_init(Xval, Yval, Npoints, MOMS);

  for(int L = 0; L < LLeg; L++){
    SEARCH_MIN_POT(dx, Npoints, &REQ[L], &VMIN[L], L);
    GET_EROT(REQ[L], L, &Erot, &dErot);
//    spline_interpolation(Xval, Yval, MOMS, Npoints, REQ[L], &vpot, &dvdx );
//    cout << "SEARCH " << " " << L << " " << REQ[L] << " " << (VMIN[L]-VMIN[0]) - Erot - vpot<< "\n";
  }
  for(int i = 0; i < Npoints; i++) Yval[i] = (Yval[i] - VMIN[0]);
  spline_init(Xval, Yval, Npoints, MOMS);  
  

  for(int L = 0; L < LLeg; L++){
//    int valL = L*Npoints;
    for(int i = 0; i < Ndvr; i++){
      double x = Xval[0] + dx*i;
      spline_interpolation(Xval, Yval, MOMS, Npoints, x, &vpot, &dvdx );
      GET_EROT(x, L, &Erot, &dErot);
      V[i] = vpot + Erot;
    }
    for(int i = 0; i < Ndvr; i++){
      for(int j = 0; j < Ndvr; j++){
        Matrix[i*Ndvr + j] = t_kin(i+1, j+1, dx, mu);
//       cout << t_kin(i+1, j+1, dx, mu) << "\n";
        if(i == j) Matrix[i*Ndvr + j] += V[i];
      }
    }
    diagonalisation(); 
    VEIGEN_ql[L]          = VEIGEN[0];
    VEIGEN_ql[LLeg   + L] = VEIGEN[1];
    VEIGEN_ql[2*LLeg + L] = VEIGEN[2];
    VEIGEN_ql[3*LLeg + L] = VEIGEN[3];
  }

  //WRITE OUT THE EIGENSTATES
  ofstream out; out.open("ro-vib_eigenvals.dat");
  out << "#QUANTUM MECHANICAL RO-VIBRATIONAL EIGENSTATES in eV\n";
  out << "#nu,   L \n";
  for(int nu = 0; nu < 4; nu++){
    out << "nu = " << nu << " ;  " << 0 << " <= L < " << LLeg << "\n";
    for(int L = 0; L < LLeg; L++){
      out << "\t\t " << setprecision(10) << fixed << setw(15) << VEIGEN_ql[nu*LLeg + L]*27.211383 << "\n";
    }
  }
  out.close();
  out.open("LLEG_POTENTIALS.dat");
  //WRITE OUT THE POTENTIALS
  for(int i = 0; i < Ndvr; i++){
    double x = i*dx + Xval[0];
    out << x << " ";
    for(int L = 0; L < LLeg; L++){
      spline_interpolation(Xval, Yval, MOMS, Npoints, x, &vpot, &dvdx );
      GET_EROT(x, L, &Erot, &dErot);
      vpot += Erot;
      out << vpot << " ";
    }
    out << "\n";
  }
  out.close();

  for(int nu = 0; nu < 4; nu++){
    for(int L = 0; L < LLeg; L++){
      DETERMINE_RMIN_RMAX(Xval, Yval, MOMS, nu, L);
//      cout << nu << " " << L << " " << RMIN[nu*LLeg + L]*0.52917721 << " " << RMAX[nu*L + L]*0.52917721 << "\n"; 
    }
  }
  //DETERMINE THE PARTITION FUNCTION
  double kB = 3.1668152e-06; //boltzmann factor Hartree/K
  PARTITION_FUNCTION = 0.; KbT = kB * TSTREAM ; //hartree/Kelvin * STREAM TEMPERATURE
  KbTrot = (-181.102 + 0.648034*TSTREAM)*kB; //ROTATIONAL THERMAL ENERGY
  if(KbTrot < kB) KbTrot = kB; //to make sure   

  for(int nu = 0; nu < 4; nu++){
    for(int L = 0; L < LLeg; L++){
      double fac11 = (VEIGEN_ql[nu*LLeg] - VEIGEN_ql[0])/KbT;
      double fac22 = (VEIGEN_ql[nu*LLeg + L] - VEIGEN_ql[nu*LLeg])/KbTrot;
      double valli   = exp(-fac11); //vibrational contribution
             valli  *= (2.0*L+1.0) * exp(-fac22); //rotational contribution assuming that Trot neq Tvib
      PARTITION_FUNCTION += valli;
//      PARTITION_FUNCTION += (2.0*L+1.0) * exp(-(VEIGEN_ql[nu*LLeg + L] - VEIGEN_ql[0])/KbT);
//    cout << nu << " " << L << " " << exp(-fac11) <<  " " << exp(-fac22) << " " << (VEIGEN_ql[nu*LLeg + L] - VEIGEN_ql[nu*LLeg]) << " " << KbTrot << "\n";
    }
 }
// cout << PARTITION_FUNCTION << "\n";
};
void QC_MOM::init_pos_mom(){
  double randnr; double dvdx; double vpot;
  int qvib; int lrot; //double veigen; 
  int qvib1, lrot1;
  int yes = 1; double val; 

  //RANDOM SELECTION OF THE QUNATUM NUMBER v,J,mJ 
  if(STATE_SWITCH == 2 || STATE_SWITCH == 3){
    do{
      if( STATE_SWITCH == 2 ){
        qvib1 = QVIBLM[0];
      }
      else{
        qvib1 = random()%4;
      }
      lrot1 = random()%LLeg;
  //    veigen = (VEIGEN_ql[qvib1*LLeg + lrot1] - VEIGEN_ql[0])/KbT;
  //    val = (2.*lrot1+1.)* exp(-veigen)/PARTITION_FUNCTION;
      double valli   = exp(-(VEIGEN_ql[qvib1*LLeg] - VEIGEN_ql[0])/KbT); //vibrational contribution
      valli  *= (2.0*lrot1+1.0) * exp(-(VEIGEN_ql[qvib1*LLeg + lrot1] - VEIGEN_ql[qvib1*LLeg])/KbTrot);
      val = valli/PARTITION_FUNCTION;

      randnr = (double)random()/RAND_MAX;
      if(randnr < val && yes == 1){
        yes  = 0;
        qvib = qvib1;
        lrot = lrot1;
      }
    }while(yes == 1);
  }
  else{
    qvib = QVIBLM[0];
    lrot = QVIBLM[1];
  }

//  lrot = 0; qvib = 0;

  add_qvib(qvib, lrot);

  int mJ = 0;
  if(STATE_SWITCH == 0){
    mJ = QVIBLM[2];
  }
  else{
    if(lrot > 0){
      mJ = random()%(2*lrot + 1) - lrot;
//      mJ = random()%( lrot + 1 );
//      mJ *= pow(-1, random()%2);
    }
  }


//  mJ = 0;
  

  //save qunatum number
  QM   = mJ;
  QVIB = qvib;
  QL   = lrot;

  //DETERMINE RANDOMLY THE ORIENTATION OF A DIATOMIC-MOLECULE  
  //AND SET THETA AND PHI VALUES
  DETERMINE_THETA_PHI(lrot, mJ);
 
  //DETERMINE R
  randnr = (double)random()/RAND_MAX;
  DETERMINE_R(qvib, lrot, &R, randnr);


  //int valL = Npoints*lrot; 
  double Erot, dErot;
  spline_interpolation(Xval, Yval, MOMS, Npoints, R, &vpot, &dvdx );
  GET_EROT(R, lrot, &Erot, &dErot);
  vpot += Erot;
  dvdx += dErot;
  double fac = VEIGEN_ql[qvib*LLeg + lrot] - vpot;// - (VMIN[lrot] - VMIN[0]);
//  cout << "INTERNALS SELECTED     " << lrot << " " << mJ << " " << Erot*27.211383 << " " << fac*27.211383 << " " << (vpot-Erot)*27.211383 << "\n";
   
  if(fac < 0) cout << "FAC IS : " << fac*27211.383 << " meV\n";
  pr = pow(-1, random()%2)*sqrt(fac * 2.*mu);

  //DETERMINE ROTATIONAL MOMEMTA
  pt = 0.; pp = 0.;
  if(lrot > 0){
     pp = (double)mJ;
       
     if(mJ != 0){
//       double cos_eta = cos(BETA)/sin(THETA);
//       cout << "COS_ETA " << cos_eta << " " << cos(BETA) << "\n";
//       double sin_eta = sin(acos(cos_eta));
//       pt = pow(-1, random()%2) * sqrt( lrot*(lrot+1.))*sin_eta ;
        pt =  pow(-1, random()%2) * sqrt( lrot*(lrot+1.)  - (double)(mJ*mJ)/(sin(THETA)*sin(THETA))) ;
     }
     else pt = pow(-1, random()%2) * sqrt(lrot*(lrot+1.));
  }
  
//  cout << "MJ " << " " << mJ << "  L " << lrot << " qvib " << qvib << "\n"; 

/*
  cout << "ROT ENERGY SELECTED 2: " 
       << ((double)(pt*pt)/(2.*mu*R*R) 
          + (double)(mJ*mJ)/(2.*mu*R*R*sin(THETA)*sin(THETA)))*27.211383 
       << " " << (double)(pt*pt)/(2.*mu*R*R) * 27.211383 << " " 
       << (double)(mJ*mJ)/(2.*mu*R*R*sin(THETA)*sin(THETA)) * 27.211383
       << "   "   << THETA*180./M_PI << " " << mu << " " << R << "\n";
*/

  double u = 2.*L_surf * ((double)random()/RAND_MAX);
  double v = 2.*L_surf * ((double)random()/RAND_MAX);
  X = (u + 0.5*v)*1.8897261;
  Y = sqrt(3.)*0.5*v*1.8897261;

  //AS DONE ACCORDING TO GEERT JAN
  //X = u*1.8897261;
  //Y = sqrt(3.)*0.5*v*1.8897261;

  //CALCULATE THE VELOCITY ALONG Z ACCORING TO A BEAM 

  yes = 1;
  if(BEAM_SWITCH == 0) yes = 0;
  while(yes == 1){
   randnr = (double)random()/RAND_MAX;
   double vz = (INT_B-INT_A)*randnr + INT_A;
   double fac = BEAM_FUNC(vz, ALPHA, VBEAM)/NORM_BEAM;
   randnr = (double)random()/RAND_MAX;
   if(randnr < fac && yes == 1){
     RAND_VBEAM = -vz/2187691.3;
     yes = 0;  
   }
  }
 
};
void QC_MOM::getposr(double* x, double* px, int i){
  if(i == 0){
    *x = X;	  
  }	
  else if(i == 1){
    *x = Y;
  }	  
  else if(i == 3){
    *x = R; *px = pr;
  }	  
  else if(i== 4){
    *x = THETA;  *px = pt;
  }	  
  else if(i== 5){
    *x = PHI; *px = pp;
  }
 
};
void QC_MOM::DETERMINE_RMIN_RMAX(double* xvals, double* yvals, double* moms, int qvib, int lrot){
  double vpot, Erot; double dvdx, dErot; 
  double rmin; double rmax; //double DVPOT = fabs(VMIN[lrot] - VMIN[0]);
  double veigen = VEIGEN_ql[qvib*LLeg + lrot];
  double req = REQ[lrot]; double UNW = 0.;

  //MAKE A PROPAGATION TO OBTAIN A TIME FOR ONE PERIOD;
  spline_interpolation(xvals, yvals, moms, Npoints, req, &vpot, &dvdx );
  GET_EROT(req, lrot, &Erot, &dErot);
  vpot += Erot;
  dvdx += dErot;
  UNW = vpot;
  rmin = req; double pr = -sqrt((veigen - vpot) * 2.*mu);
  double x = req; double dt = 0.001; double conv = 0.;
  do{
    x += (pr - 0.5*dvdx*dt)*dt/mu;
    double dvdx_old = dvdx;
    spline_interpolation(xvals, yvals, moms, Npoints, x, &vpot, &dvdx );
    GET_EROT(x, lrot, &Erot, &dErot);
    vpot += Erot;
    dvdx += dErot;
    pr += -0.5*(dvdx_old + dvdx)*dt;
    if(x < rmin) rmin = x;
    conv = fabs(veigen - vpot);
  }while(conv > 1e-07);

  x = req; rmax = 0.;
  spline_interpolation(xvals, yvals, moms, Npoints, req, &vpot, &dvdx );
  GET_EROT(req, lrot, &Erot, &dErot);
  vpot += Erot;
  dvdx += dErot;
  pr = sqrt((veigen - vpot) * 2.*mu);
  do{
    x += (pr - 0.5*dvdx*dt)*dt/mu;
    double dvdx_old = dvdx;
    spline_interpolation(xvals, yvals, moms, Npoints, x, &vpot, &dvdx );
    GET_EROT(x, lrot, &Erot, &dErot);
    vpot += Erot;
    dvdx += dErot;
    pr += -0.5*(dvdx_old + dvdx)*dt;
    if(rmax < x) rmax = x;
    conv = fabs(veigen - vpot);
  }while(conv > 1e-07);
  RMIN[qvib*LLeg + lrot] = rmin;
  RMAX[qvib*LLeg + lrot] = rmax;
};

void QC_MOM::DETERMINE_R(int qvib, int lrot, double* R_out, double rand_nr){
  //MAKE A SHORT MD FROM RMIN TO RMAX 
  //AND DETERMINE FINAL TIME AND NUMBER OF STEPS  
  double rmin = RMIN[qvib*LLeg + lrot];
  double rmax = RMAX[qvib*LLeg + lrot];
  //double DVPOT = fabs(VMIN[lrot] - VMIN[0]);
  //double veigen = VEIGEN_ql[qvib*LLeg + lrot];
  double vpot, Erot; double dvdx, dErot; double x;



  spline_interpolation(Xval, Yval, MOMS, Npoints, rmin, &vpot, &dvdx );
  GET_EROT(rmin, lrot, &Erot, &dErot);
  vpot += Erot;
  dvdx += dErot;

  pr = 0.; 

  x = rmin; int count = 0;
  double dt = 0.5; //FROM RMIN TO RMAX
  double conv = 1.0;
  do{
    x += (pr - 0.5*dvdx*dt)*dt/mu;
    double dvdx_old = dvdx;
    spline_interpolation(Xval, Yval, MOMS, Npoints, x, &vpot, &dvdx );
    GET_EROT(x, lrot ,&Erot, &dErot);
    vpot += Erot;
    dvdx += dErot;

    pr += -0.5*(dvdx_old + dvdx)*dt;
    if(fabs(x-rmax) < 1e-05) conv = 0.;
    count++;
  }while(conv > 1e-07);
  spline_interpolation(Xval, Yval, MOMS, Npoints, rmin, &vpot, &dvdx );
  GET_EROT(rmin, lrot, &Erot, &dErot);
  pr = 0.;
  x = rmin;
  vpot += Erot;
  dvdx += dErot;

  int N_MD = count;

  double* MOMS_MD = new double[N_MD];
  double* MD_Y    = new double[N_MD];
  double* MD_X    = new double[N_MD];
  dt = 1.0; //THIS MAKES SURE THAT WE HAVE ONE PERIOD FROM RMIN -> RMAX -> RMIN
  for(int i = 0; i < N_MD; i++){
    MD_Y[i] = x ;
    MD_X[i] = dt * i;
    x += (pr - 0.5*dvdx*dt)*dt/mu;
    double dvdx_old = dvdx;
    spline_interpolation(Xval, Yval, MOMS, Npoints, x, &vpot, &dvdx );
    GET_EROT(x, lrot, &Erot, &dErot);
    vpot += Erot;
    dvdx += dErot;
    pr += -0.5*(dvdx_old + dvdx)*dt;
  }
  double MD_TIME = dt*(N_MD-1);
  spline_init(MD_X, MD_Y, N_MD, MOMS_MD);
  double randnr = ((double)random()/RAND_MAX);
  double t = MD_TIME*randnr;
  spline_interpolation(MD_X, MD_Y, MOMS_MD, N_MD, t, &x, &dvdx );
  *R_out = x;
  delete MOMS_MD; delete MD_Y; delete MD_X;
};

void QC_MOM::DETERMINE_THETA_PHI(int lrot, int mJ){

  double alpha, beta, gamma;
  double cos_theta; double sin_theta; double cos_phi; double sin_phi;
  double cos_alpha; double sin_alpha; double cos_gamma; double sin_gamma;
  double cos_beta; double sin_beta; double val;

   if(lrot > 0){
    cos_beta = 0.;
    if(fabs(mJ) > 0){
      cos_beta = ((double)mJ) /(sqrt(lrot*(lrot+1.0)));
      if(fabs(cos_beta) > 1.0) cos_beta = (int)cos_beta;
    }
//    SIN_BETA = sqrt(1.0 - cos_beta*cos_beta)
    beta = acos(cos_beta); sin_beta = sin(beta);
    BETA = beta;
    alpha = 2.*M_PI*((double)random()/RAND_MAX);
    gamma = 2.*M_PI*((double)random()/RAND_MAX);

    cos_gamma = cos(gamma); sin_gamma = sin(gamma);
    cos_alpha = cos(alpha); sin_alpha = sin(alpha);
    cos_theta = -sin_beta * cos_gamma;
    if(fabs(cos_theta) > 1.0) cos_theta = (int)cos_theta;
    THETA = acos(cos_theta); sin_theta = sin(THETA);

    sin_phi = 0.0; cos_phi = 1.0;
    if(fabs(sin_theta) > 1e-09){
      sin_phi = cos_beta*sin_alpha*cos_gamma + cos_alpha*sin_gamma;
      sin_phi /= sin_theta;
      cos_phi = cos_beta*cos_alpha*cos_gamma - sin_alpha*sin_gamma;
      cos_phi /= sin_theta;
    }
    if(fabs(cos_phi) > 1.0) cos_phi = (int)cos_phi;
//    if(fabs(sin_phi) > 1.0) sin_phi = (int)sin_phi;
    val = acos(cos_phi);
    if(fabs(sin_phi) != sin_phi) val = 2.*M_PI - val;
    PHI = val;
  }
  else{
    BETA = M_PI/2.;
    cos_theta = -1. + 2.*((double)random()/RAND_MAX);
    if(fabs(cos_theta ) > 1.0) cos_theta = (int)cos_theta;
    THETA = acos(cos_theta);
    PHI   =  2.*M_PI*((double)random()/RAND_MAX);
  }

};