#include "constants.txt"


//******************************************//
//*******Analytical solution Function*******//
//******************************************//
double OscP_JC(double *x, double *par)
{// begin

  double P = 0.0;
  
  double Enu = x[0];  // GeV
  
  double za  = par[0];
  double t1  = par[1];
  double t2  = par[2]; 
  int   idx  = par[3]; // 0 mu->e, 1 mu->mu, 2 mu->tau
  int   sel  = par[4]; // 0 Magnus 1st; 1 Magnus 2nd order
  
  double dCP = deltaCP*c_pi/180.;
  TComplex gamma33   = TComplex(cos(dCP), sin(dCP));
  TComplex gamma33cc = TComplex(cos(dCP),-sin(dCP));
  //**Gamma matrix
  TComplex GMatrix[3][3];
  TComplex GMatrixcc[3][3];
  for(int i = 0 ; i < 3 ; i++)
    for(int j = 0 ; j < 3 ; j++){
      GMatrix[i][j]   = TComplex(0.0,0.0);
      GMatrixcc[i][j] = TComplex(0.0,0.0);
    }
  
  GMatrix[0][0]   = TComplex(1.0,0.0);
  GMatrix[1][1]   = TComplex(1.0,0.0);
  GMatrix[2][2]   = gamma33;
  GMatrixcc[0][0] = TComplex(1.0,0.0);
  GMatrixcc[1][1] = TComplex(1.0,0.0);
  GMatrixcc[2][2] = gamma33cc;

  double re, im;
  
  //-the phase factors: alpha1, alpha2, alpha3 (Integrated approximated eigenvalues)
  
  TF1 *fffI3t = new TF1("fffI3t",alphaL_JC,-tf_max,tf_max,4); 
  fffI3t->SetParameter(0,za); 
  fffI3t->SetParameter(1,Enu); 
  fffI3t->SetParameter(2,t1); 
  fffI3t->SetParameter(3,1);   // select I1(=alphaL_1 = alpha1)
  double alpha1 = fffI3t->Eval(t2);
  delete fffI3t;
  
  TF1 *fffI12t = new TF1("fffI12t",alphaL_JC,-tf_max,tf_max,4); 
  fffI12t->SetParameter(0,za); 
  fffI12t->SetParameter(1,Enu); 
  fffI12t->SetParameter(2,t1); 
  fffI12t->SetParameter(3,3);   // select I_3(= alpha2)
  double alpha2 = fffI12t->Eval(t2);
  delete fffI12t;
  
  TF1 *fffI1t = new TF1("fffI1t",alphaH_JC,-tf_max,tf_max,4); 
  fffI1t->SetParameter(0,za); 
  fffI1t->SetParameter(1,Enu); 
  fffI1t->SetParameter(2,t1); 
  fffI1t->SetParameter(3,3);   // select I3(=alphaH_3 = alpha3)
  double alpha3 = fffI1t->Eval(t2);
  delete fffI1t;
  
  //**P(tf,t0) matrix. Added by M.A.AceroO. (9.Apr.2011)
  TComplex PMatrix[3][3];
  for(int i = 0 ; i < 3 ; i++)
    for(int j = 0 ; j < 3 ; j ++)
      PMatrix[i][j] = TComplex(0.0,0.0);
  PMatrix[0][0] = TComplex(cos(alpha1),-sin(alpha1));
  PMatrix[1][1] = TComplex(cos(alpha2),-sin(alpha2));
  PMatrix[2][2] = TComplex(cos(alpha3),-sin(alpha3));
  
  //- the mixing matrix in matter
  TComplex Um[3][3];
  TComplex Umcc[3][3];
  
  TF1 *ffUm = new TF1("ffUm",Umatt_JC,-tf_max,tf_max,5);
  ffUm->SetParameter(0,za);
  ffUm->SetParameter(1,Enu);

  
  for (int i=0;i<3;i++){
    for (int j=0;j<3;j++){
      
      ffUm->SetParameter(2,i);  // flav idx
      ffUm->SetParameter(3,j);  // mass idx
      ffUm->SetParameter(4,0);  // Re part
      re = ffUm->Eval(t1);
      ffUm->SetParameter(4,1);  // Im part 
      im = ffUm->Eval(t1);
      Um[i][j] = TComplex(re,im);
      Umcc[i][j] = TComplex(re,-im);
    }
  }

  delete ffUm;
  
  TComplex ULop0[3];
  TComplex ULop0cc[3];  
  
  TF1 *ffuLOp0 = new TF1("ffuLOp0",u_33_11_12_JC,0,Enu_max,6);
  ffuLOp0->SetParameter(0,za);
  ffuLOp0->SetParameter(1,t1); //*Lower limit set to the initial time
  ffuLOp0->SetParameter(2,0.); //*Upper limit set to 'tbar'
  ffuLOp0->SetParameter(5,sel);  

  
  for (int i=0;i<3;i++){
    ffuLOp0->SetParameter(3,i);  
    ffuLOp0->SetParameter(4,0);   // Re part
    re = ffuLOp0->Eval(Enu);
    ffuLOp0->SetParameter(4,1);   // Im part
    im = ffuLOp0->Eval(Enu);
    
    ULop0[i]   = TComplex(re,im);
    ULop0cc[i] = TComplex(re,-im);    
  }

    delete ffuLOp0;
  
  //**Added by M.A.AceroO. (8.Apr.2011)*
  TComplex UL0_Mat[3][3];
  for(int i = 0 ; i < 3 ; i++)
    for(int j = 0 ; j < 3 ; j++){
      UL0_Mat[i][j]   = TComplex(0.0,0.0);
    }      
  UL0_Mat[0][0] = ULop0[1];
  UL0_Mat[1][1] = ULop0cc[1];
  UL0_Mat[0][1] = ULop0[2];
  UL0_Mat[1][0] =-ULop0cc[2];
  UL0_Mat[2][2] = TComplex(1.0,0.0);
  
  
  //************************************
  //**Added by M.A.AceroO. (27.Jul.2011)*
  double real = 0.0;
  double imag = 0.0;
  TF1 *fun_Ph21t = new TF1("fun_Ph21t",Phi21ofT_JC,-tf_max,tf_max,3);
  fun_Ph21t->SetParameter(0,za); 
  fun_Ph21t->SetParameter(1,Enu); 
  //fun_Ph21t->SetParameter(2,t1);
  fun_Ph21t->SetParameter(2,tlohi[0]);
  double Phi21 = fun_Ph21t->Eval(0.);

  delete fun_Ph21t;

  real = cos(2.0*Phi21);
  imag = sin(2.0*Phi21);
  TComplex Phi_phase   = TComplex(real,-imag);
  TComplex Phi_phase_CC= TComplex(real, imag);

  TComplex ULf_Mat[3][3];
  for(int i = 0 ; i < 3 ; i++)
    for(int j = 0 ; j < 3 ; j++){
      ULf_Mat[i][j]   = TComplex(0.0,0.0);
    }      
  ULf_Mat[0][0] = UL0_Mat[0][0];
  ULf_Mat[1][1] = UL0_Mat[1][1];
  ULf_Mat[0][1] = UL0_Mat[1][0]*Phi_phase;
  ULf_Mat[1][0] = UL0_Mat[0][1]*Phi_phase_CC;
  ULf_Mat[2][2] = TComplex(1.0,0.0);
  
  //************************************
  
  TComplex UHop[3];
  TComplex UHopcc[3];  
  
  TF1 *ffuHOp = new TF1("ffuHOp",u_11_22_23_JC,0,Enu_max,6);
  ffuHOp->SetParameter(0,za);
  ffuHOp->SetParameter(1,t1);
  ffuHOp->SetParameter(2,t2);
  ffuHOp->SetParameter(5,sel);  
  
  for (int i=0;i<3;i++){
    ffuHOp->SetParameter(3,i);  
    ffuHOp->SetParameter(4,0);   // Re part
    re = ffuHOp->Eval(Enu);
    ffuHOp->SetParameter(4,1);   // Im part
    im = ffuHOp->Eval(Enu);
    
    UHop[i] = TComplex(re,im);
    UHopcc[i] = TComplex(re,-im);
  }
  
  delete ffuHOp;

  //**Added by M.A.AceroO. (8.Apr.2011)*
  TComplex UH_Mat[3][3];
  for(int i = 0 ; i < 3 ; i++)
    for(int j = 0 ; j < 3 ; j++){
      UH_Mat[i][j]   = TComplex(0.0,0.0);
    }      
  UH_Mat[0][0] = UHop[0];
  UH_Mat[1][1] = UHop[1];
  UH_Mat[2][2] = UHopcc[1];
  UH_Mat[1][2] = UHop[2];
  UH_Mat[2][1] =-UHopcc[2];
  
  //************************************
  
  //**Added by M.A.AceroO. (8.Apr.2011)*
  //**Matrix to define the Up operator
  //**Up = UL(tf,tbar)*UH(tf,t0)*UL(tbar,t0)
  TComplex UPoperator[3][3],temp[3][3];
  for(int i = 0 ; i < 3 ; i++){
    for(int j = 0 ; j < 3 ; j++){
      UPoperator[i][j] = TComplex(0.0,0.0);
    }
  }
  
  for(int i = 0 ; i < 3 ; i++)
    for(int j = 0 ; j < 3 ; j++)
      for(int k = 0 ; k < 3 ; k++)
	for(int l = 0 ; l < 3 ; l++)
	  UPoperator[i][j] += ULf_Mat[i][k]*UH_Mat[k][l]*UL0_Mat[l][j];
  
  //**Total Evolution operator UT(tf,t0). Added by M.A.AceroO. (9.Apr.2011)
  //**UT(tf,t0) = Gamma*PMatrix(tf,t0)*Up(tf,t0)*Gamma^+
  TComplex EvolOperator[3][3],A1[3][3],B1[3][3];
  for(int i = 0 ; i < 3 ; i++){
    for(int j = 0 ; j < 3 ; j++){
      EvolOperator[i][j] = TComplex(0.0,0.0);  
      A1[i][j]           = TComplex(0.0,0.0);  
      B1[i][j]           = TComplex(0.0,0.0);  
    }
  }
  
  for(int i = 0 ; i < 3 ; i++)
    for(int j = 0 ; j < 3 ; j++)
      for(int k = 0 ; k < 3 ; k++)
	A1[i][j] += GMatrix[i][k]*PMatrix[k][j];
  
  for(int k = 0 ; k < 3 ; k++)
    for(int l = 0 ; l < 3 ; l++)
      for(int m = 0 ; m < 3 ; m++)
	B1[k][l] += UPoperator[k][m]*GMatrixcc[m][l];
  
  for(int i = 0 ; i < 3 ; i++)
    for(int j = 0 ; j < 3 ; j++)
      for(int k = 0 ; k < 3 ; k++)
	EvolOperator[i][j] += A1[i][k]*B1[k][j];
  
  //calculate the oscillation amplitudes
  
  TComplex Amp[3]; //**Modified by M.A.AceroO. (9.Abr.2011)
  
  for (int i = 0 ; i < 3 ; i++){
    
    Amp[i] = 
      Um[i][0]*
      ( Umcc[1][0]*EvolOperator[0][0] + Umcc[1][1]*EvolOperator[0][1] + Umcc[1][2]*EvolOperator[0][2] )
      +
      Um[i][1]*
      ( Umcc[1][0]*EvolOperator[1][0] + Umcc[1][1]*EvolOperator[1][1] + Umcc[1][2]*EvolOperator[1][2] )
      +
      Um[i][2]*
      ( Umcc[1][0]*EvolOperator[2][0] + Umcc[1][1]*EvolOperator[2][1] + Umcc[1][2]*EvolOperator[2][2] );
  }
  
  if      (idx == 0) P = Amp[0].Rho2();
  else if (idx == 1) P = Amp[1].Rho2();
  else if (idx == 2) P = Amp[2].Rho2();
  
  printf("The oscilation probability is P = %5.2f \n",P);

  return P;
  
}// end

//-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
//-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-

double u_33_11_12_JC(double *x, double *par)
{// begin
  
  double Enu = x[0];  // GeV
  double U = 0.0;

  double za  = par[0];
  double t1  = par[1];
  double t2  = par[2];
  int  idx   = par[3]; // 0 u33, 1 u11, 2 u12
  int  reim  = par[4]; // 0 Re part, 1 Im part, 2 Im part of cc
  int  sel   = par[5]; // 0 1st order Magnus; 1 2nd order Magnus

  double dcp = deltaCP*c_pi/180.0;

  double s = 1.;
  if (sel == 0) s = 0.;

  TF1 *ggggPh21t = new TF1("ggggPh21t",Phi21ofT_JC,-tf_max,tf_max,3);
  ggggPh21t->SetParameter(0,za); 
  ggggPh21t->SetParameter(1,Enu); 
  ggggPh21t->SetParameter(2,t1);
  double Phi21 = ggggPh21t->Eval(t2); //** t2 = tbar = 0.0
  delete ggggPh21t;
  //**************************************************************************//

  TF1 *ggXi12Enu = new TF1("ggXi12Enu",Xi12LowofE_JC,0,Enu_max,4);
  ggXi12Enu->SetParameter(0,za);
  ggXi12Enu->SetParameter(1,t1);
  ggXi12Enu->SetParameter(2,t2); //upper limit set to tbar (= t2)

  ggXi12Enu->SetParameter(3,0);  // re part (xi_y)
  double Xiy = ggXi12Enu->Eval(Enu);
  ggXi12Enu->SetParameter(3,1);  // im part (xi_x)
  double Xix = -1.0*ggXi12Enu->Eval(Enu);

  delete ggXi12Enu;

  TF1 *ggXi3Enu = new TF1("ggXi3Enu",Xi3ofE_JC,0,Enu_max,4); 
  ggXi3Enu->SetParameter(0,za);
  ggXi3Enu->SetParameter(1,t1);
  ggXi3Enu->SetParameter(2,t2); //upper limit set to tbar (= t2)
  ggXi3Enu->SetParameter(3,0);
  double Xi3 = s*ggXi3Enu->Eval(Enu);
  delete ggXi3Enu;
  
  double Xi12 = Xiy**2 + Xix**2;

  //double Xi = sqrt( Xi12**2 + Xi3**2 );
  double Xi = sqrt( Xi12 + Xi3**2 );
  
  double sinXiovXi = 1.0;
  if ( Xi != 0 ) sinXiovXi = sin(Xi)/Xi;
  
  if (idx == 0) {
    if      (reim == 0) U =  cos(0.0);
    else if (reim == 1) U = -sin(0.0);
    else if (reim == 2) U =  sin(0.0);
  }
  else if (idx == 1) {
    if      (reim == 0) U =  cos(Xi);
    else if (reim == 1) U = -sinXiovXi*Xi3;
    else if (reim == 2) U =  sinXiovXi*Xi3;
  }
  else if (idx == 2) {
    if      (reim == 0) U = -sinXiovXi*Xiy;
    else if (reim == 1) U = -sinXiovXi*Xix;
    else if (reim == 2) U =  sinXiovXi*Xix;
  }
  //*****************************************************************************//

  return U;
} // end

//-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
//-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-

double u_11_22_23_JC(double *x, double *par)
{// begin

  double Enu = x[0];  // GeV
  double U = 0.0;

  double za  = par[0];
  double t1  = par[1];
  double t2  = par[2];
  int  idx   = par[3]; // 0 u11, 1 u22, 2 u23
  int  reim  = par[4]; // 0 Re part, 1 Im part, 2 Im part of cc
  int  sel   = par[5]; // 0 1st order Magnus; 1 2nd order Magnus

  double dcp = deltaCP*c_pi/180.0;

  double s = 1.;
  if (sel == 0) s = 0.;

  TF1 *ffI1t = new TF1("ffI1t",alphaH_JC,-tf_max,tf_max,4);
  ffI1t->SetParameter(0,za); 
  ffI1t->SetParameter(1,Enu); 
  ffI1t->SetParameter(2,t1); 
  ffI1t->SetParameter(3,1);    // selects alpha1_H = I1
  double I1 = ffI1t->Eval(t2);
  delete ffI1t;

  TF1 *ffI23t = new TF1("ffI23t",alphaH_JC,-tf_max,tf_max,4); 
  ffI23t->SetParameter(0,za); 
  ffI23t->SetParameter(1,Enu); 
  ffI23t->SetParameter(2,t1); 
  ffI23t->SetParameter(3,4); // selects I23
  double I23 = ffI23t->Eval(t2);
  delete ffI23t;
  
  TF1 *ffffPh32t = new TF1("ffffPh32t",Phi32ofT_JC,-tf_max,tf_max,3);
  ffffPh32t->SetParameter(0,za); 
  ffffPh32t->SetParameter(1,Enu); 
  ffffPh32t->SetParameter(2,0); // low lim of integ for phi set to "tbar"
  //ffffPh32t->SetParameter(2,t1);
  double Phi32 = ffffPh32t->Eval(t2);
  delete ffffPh32t;

  TF1 *ffXi12Enu = new TF1("ffXi12Enu",Xi12HighofE_JC,0,Enu_max,3);
  ffXi12Enu->SetParameter(0,za);
  ffXi12Enu->SetParameter(1,0);  // low lim of integ set to "tbar"
  //ffXi12Enu->SetParameter(1,t1);
  ffXi12Enu->SetParameter(2,t2);
  double Xi12 = ffXi12Enu->Eval(Enu);
  delete ffXi12Enu;

  TF1 *ffXi3Enu = new TF1("ffXi3Enu",Xi3ofE_JC,0,Enu_max,4); 
  ffXi3Enu->SetParameter(0,za);
  ffXi3Enu->SetParameter(1,t1);
  ffXi3Enu->SetParameter(2,t2);
  ffXi3Enu->SetParameter(3,1);
  double Xi3 = s*ffXi3Enu->Eval(Enu);
  delete ffXi3Enu;

  double Xi = sqrt( Xi12**2 + Xi3**2 );

  double sinXiovXi = 1.0;
  if ( Xi != 0 ) sinXiovXi = sin(Xi)/Xi;

  if (idx == 0) {
    if      (reim == 0) U =  cos(0.0);
    else if (reim == 1) U = -sin(0.0);
    else if (reim == 2) U =  sin(0.0);
  }
  else if (idx == 1) {
    if      (reim == 0) U = cos(Xi);
    else if (reim == 1) U =-sinXiovXi*Xi3;
    else if (reim == 2) U = sinXiovXi*Xi3;
  }
  else if (idx == 2) {
    if      (reim == 0) U = sin(Phi32)*sinXiovXi*Xi12;
    else if (reim == 1) U = cos(Phi32)*sinXiovXi*Xi12;
    else if (reim == 2) U =-cos(Phi32)*sinXiovXi*Xi12;
  }

  return U;
} // end

//-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
//-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-

double Xi3ofE_JC(double *x, double *par)
{// begin

  double Enu = x[0]; //GeV

  double za   = par[0];
  double t1   = par[1];
  double t2   = par[2];
  double sel  = par[3];  // 0 low, 1 high
  //int sel2    = par[4];

  double X = 0.0;
TF1 *fffThmijt = 0;
  if (sel == 0){
    fffThmijt = new TF1("fffThmijt",Thetam12_JC,-tf_max,tf_max,2);
    fffThmijt->SetParameter(0,za); 
    fffThmijt->SetParameter(1,Enu); 

  }else if (sel == 1){
    fffThmijt = new TF1("fffThmijt",Thetam13_JC,-tf_max,tf_max,2);
    fffThmijt->SetParameter(0,za); 
    fffThmijt->SetParameter(1,Enu); 

  }
  
  TF1 *ffXi3At = new TF1("ffXi3At",Xi3_AofT_JC,-tf_max,tf_max,4);
  ffXi3At->SetParameter(0,za); 
  ffXi3At->SetParameter(1,Enu); 
  ffXi3At->SetParameter(2,t1);
  ffXi3At->SetParameter(3,sel);

  //- Calculate integral using delta functions

  double II = 0.0;
  double z  = 1.0;
  for (int i=1; i<2*Nshells;i++){
    double tc = times[2*Nshells-1 -  i  ];
    double tn = times[2*Nshells-1 -(i-1)];
    if (i == 1) z = 0.0; else z = 1.0;
    //Checking:
    //printf("\n sel2 = %d, t1 = %f , t2 = %f \n",sel2,t1,t2);
    if (t1<tc)
      II += (
	     fffThmijt->Eval(0.5*(tc+tn))*(c_pi/180.)*
	     ( ffXi3At->Eval(tc)-z*ffXi3At->Eval(tn) )
	     );
    else if (t1<tn)
      II += (
	     fffThmijt->Eval(0.5*(tc+tn))*(c_pi/180.)*
	     ( -z*ffXi3At->Eval(tn) )
	     );

    if (t2<tc)
      II -= (
	     fffThmijt->Eval(0.5*(tc+tn))*(c_pi/180.)*
	     ( ffXi3At->Eval(tc)-z*ffXi3At->Eval(tn) )
	     );
    else if (t2<tn)
      II -= (
	     fffThmijt->Eval(0.5*(tc+tn))*(c_pi/180.)*
	     ( -z*ffXi3At->Eval(tn) )
	     );
  }
  X = II;  
delete fffThmijt;
  delete ffXi3At;

  return X;
}// end

//-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
//-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-

double Xi3_AofT_JC(double *x, double *par)
{// begin

  double t2 = x[0];

  double za  = par[0];
  double Enu = par[1];
  double t1  = par[2];
  int    sel = par[3]; // 0 low, 1 high
  
  double X = 0.0;
  TF1 *ffThmijt = 0;
  TF1 *fffhPhjkt = 0;
  if (sel == 0 ){
    ffThmijt = new TF1("ffThmijt",Thetam12_JC,-tf_max,tf_max,2);
    ffThmijt->SetParameter(0,za); 
    ffThmijt->SetParameter(1,Enu); 

    fffhPhjkt = new TF1("fffhPhjkt",Phi21ofT_JC,-tf_max,tf_max,3);
    fffhPhjkt->SetParameter(0,za);
    fffhPhjkt->SetParameter(1,Enu); 
    fffhPhjkt->SetParameter(2,t2);     
  }
  else if (sel == 1){
    ffThmijt = new TF1("ffThmijt",Thetam13_JC,-tf_max,tf_max,2);
    ffThmijt->SetParameter(0,za); 
    ffThmijt->SetParameter(1,Enu); 

    
    fffhPhjkt = new TF1("fffhPhjkt",Phi32ofT_JC,-tf_max,tf_max,3);
    fffhPhjkt->SetParameter(0,za);
    fffhPhjkt->SetParameter(1,Enu); 
    fffhPhjkt->SetParameter(2,t2); 
  }

  //- Calculate integral using delta functions

  double II = 0.0;
  double z  = 1.0;
  for (int i=1; i<2*Nshells;i++){
    double tc = times[2*Nshells-1 -  i  ];
    double tn = times[2*Nshells-1 -(i-1)];
    if (i == 1) z = 0.0; else z = 1.0;
    if (t1<tc)
      II += (
	     ffThmijt->Eval(0.5*(tc+tn))*(c_pi/180.)*
	     ( sin(fffhPhjkt->Eval(tc))-z*sin(fffhPhjkt->Eval(tn)) )
	     );
    else if (t1<tn)
      II += (
	     ffThmijt->Eval(0.5*(tc+tn))*(c_pi/180.)*
	     ( -z*sin(fffhPhjkt->Eval(tn)) )
	     );

    if (t2<tc)
      II -= (
	     ffThmijt->Eval(0.5*(tc+tn))*(c_pi/180.)*
	     ( sin(fffhPhjkt->Eval(tc))-z*sin(fffhPhjkt->Eval(tn)) )
	     );
    else if (t2<tn)
      II -= (
	     ffThmijt->Eval(0.5*(tc+tn))*(c_pi/180.)*
	     ( -z*sin(fffhPhjkt->Eval(tn)) )
	     );
  }
  X = II;

  delete ffThmijt;
  delete fffhPhjkt;

  return X;

}// end

//-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
//-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-

double Xi12LowofE_JC(double *x, double *par)//**Modified version by M.A.AceroO. (27.Jul.2011)
{// begin

  double Enu = x[0]; //GeV

  double za = par[0];
  double t1 = par[1];
  double t2 = par[2];
  int   sel = par[3]; // 0 re; 1 im

  double Xl = 0.0;
  double II = 0.0;
  double z  = 1.0;

  TF1 *ffThmijt = new TF1("ffThmijt",Thetam12_JC,-tf_max,tf_max,2);
  ffThmijt->SetParameter(0,za);
  ffThmijt->SetParameter(1,Enu);
  
  TF1 *fffPhjkt = new TF1("fffPhjkt",Phi21ofT_JC,-tf_max,tf_max,3);
  fffPhjkt->SetParameter(0,za);
  fffPhjkt->SetParameter(1,Enu);
  fffPhjkt->SetParameter(2,t1);  //**29.Jul.2011
  
  //- Calculate integral using delta functions
  
  for (int i=1; i<2*Nshells;i++){
    double tc = times[2*Nshells-1 -  i  ];
    double tn = times[2*Nshells-1 -(i-1)];
    
    if (i == 1) z = 0.0; else z = 1.0;
    if (t1<tc){
      if (sel == 1)
	II += (
	       ffThmijt->Eval(0.5*(tc+tn))*(c_pi/180.)*
	       ( sin(fffPhjkt->Eval(tc))-z*sin(fffPhjkt->Eval(tn)) )
	       );
      else if (sel == 0)
	II += (
	       ffThmijt->Eval(0.5*(tc+tn))*(c_pi/180.)*
	       ( cos(fffPhjkt->Eval(tc))-z*cos(fffPhjkt->Eval(tn)) )
	       );
    }
    else if (t1<tn){
      if (sel == 1)
	II += (
	       ffThmijt->Eval(0.5*(tc+tn))*(c_pi/180.)*
	       ( -z*sin(fffPhjkt->Eval(tn)) )
	       );
      else if (sel == 0)
	II += (
	       ffThmijt->Eval(0.5*(tc+tn))*(c_pi/180.)*
	       ( -z*cos(fffPhjkt->Eval(tn)) )
	       );
    }
    if (t2<tc){
      if (sel == 1)
	II -= (
	       ffThmijt->Eval(0.5*(tc+tn))*(c_pi/180.)*
	       ( sin(fffPhjkt->Eval(tc))-z*sin(fffPhjkt->Eval(tn)) )
	       );
      else if (sel == 0)
	II -= (
	       ffThmijt->Eval(0.5*(tc+tn))*(c_pi/180.)*
	       ( cos(fffPhjkt->Eval(tc))-z*cos(fffPhjkt->Eval(tn)) )
	       );
    }
    else if (t2<tn){
      if (sel == 1)
	II -= (
	       ffThmijt->Eval(0.5*(tc+tn))*(c_pi/180.)*
	       ( -z*sin(fffPhjkt->Eval(tn)) )
	       );
      else if (sel == 0)
	II -= (
	       ffThmijt->Eval(0.5*(tc+tn))*(c_pi/180.)*
	       ( -z*cos(fffPhjkt->Eval(tn)) )
	       );
    }
  }
  Xl = II;
  delete ffThmijt;
  delete fffPhjkt;
  return Xl;

}// end

//-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
//-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-

double Xi12HighofE_JC(double *x, double *par)
{// begin

  double Enu = x[0]; //GeV

  double za = par[0];
  double t1 = par[1];
  double t2 = par[2];

  double Xh = 0.0;
  double II = 0.0;
  double z  = 1.0;

  TF1 *ffThmijt = new TF1("ffThmijt",Thetam13_JC,-tf_max,tf_max,2);
  ffThmijt->SetParameter(0,za); 
  ffThmijt->SetParameter(1,Enu); 
  
  TF1 *fffPhjkt = new TF1("fffPhjkt",Phi32ofT_JC,-tf_max,tf_max,3); 
  fffPhjkt->SetParameter(0,za); 
  fffPhjkt->SetParameter(1,Enu); 
  fffPhjkt->SetParameter(2,t1); 
  
  //- Calculate integral using delta functions
  
  for (int i=1; i<2*Nshells;i++){
    double tc = times[2*Nshells-1 -  i  ];
    double tn = times[2*Nshells-1 -(i-1)];
    
    if (i == 1) z = 0.0; else z = 1.0;
    if (t1<tc)
      II += (
	     ffThmijt->Eval(0.5*(tc+tn))*(c_pi/180.)*
	     ( sin(fffPhjkt->Eval(tc))-z*sin(fffPhjkt->Eval(tn)) )
	     );
    else if (t1<tn)
      II += (
	     ffThmijt->Eval(0.5*(tc+tn))*(c_pi/180.)*
	     ( -z*sin(fffPhjkt->Eval(tn)) )
	     );
    
    if (t2<tc)
      II -= (
	     ffThmijt->Eval(0.5*(tc+tn))*(c_pi/180.)*
	     ( sin(fffPhjkt->Eval(tc))-z*sin(fffPhjkt->Eval(tn)) )
	     );
    else if (t2<tn)
      II -= (
	     ffThmijt->Eval(0.5*(tc+tn))*(c_pi/180.)*
	     ( -z*sin(fffPhjkt->Eval(tn)) )
	     );
  }
  Xh = 2*II;
  delete ffThmijt;
  delete fffPhjkt;
  return Xh;

}// end

//-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
//-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-

double Phi21ofT_JC(double *x, double *par)
{// begin

  double P = 0.0;
  double t = x[0];

  double za  = par[0];
  double Enu = par[1];
  double t1  = par[2];

  TF1 *llDm21_t = new TF1("llDm21_t",Deltam21_JC,-tf_max,tf_max,2); 
  llDm21_t->SetParameter(0,za);
  llDm21_t->SetParameter(1,Enu); 

  //- Calculate area under step function Deltam32

  for (int i=1; i<2*Nshells;i++){
    double tc = times[2*Nshells-1 -  i  ];
    double tn = times[2*Nshells-1 -(i-1)];
    if (t1<tc)
      P += llDm21_t->Eval(0.5*(tc+tn))*(tn-tc)*c_light/hbarc;
    else if (t1<tn)
      P += llDm21_t->Eval(0.5*(t1+tn))*(tn-t1)*c_light/hbarc;

    if (t<tc)
      P -= llDm21_t->Eval(0.5*(tc+tn))*(tn-tc)*c_light/hbarc;
    else if (t<tn)
      P -= llDm21_t->Eval(0.5*(t+tn))*(tn-t)*c_light/hbarc;
  }
delete llDm21_t;
  return P;

}// end

//-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
//-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-

double Phi32ofT_JC(double *x, double *par)
{// begin

  double P = 0.0;
  double t = x[0];

  double za  = par[0];
  double Enu = par[1];
  double t1  = par[2];

  TF1 *ffDm32_t = new TF1("ffDm32_t",Deltam32_JC,-tf_max,tf_max,2); 
  ffDm32_t->SetParameter(0,za); 
  ffDm32_t->SetParameter(1,Enu); 

  //- Calculate area under step function Deltam32

  for (int i=1; i<2*Nshells;i++){
    double tc = times[2*Nshells-1 -  i  ];
    double tn = times[2*Nshells-1 -(i-1)];
    if (t1<tc)
      P += ffDm32_t->Eval(0.5*(tc+tn))*(tn-tc)*c_light/hbarc;
    else if (t1<tn)
      P += ffDm32_t->Eval(0.5*(t1+tn))*(tn-t1)*c_light/hbarc;

    if (t<tc)
      P -= ffDm32_t->Eval(0.5*(tc+tn))*(tn-tc)*c_light/hbarc;
    else if (t<tn)
      P -= ffDm32_t->Eval(0.5*(t+tn))*(tn-t)*c_light/hbarc;
  }
  delete ffDm32_t;
  return P;


}// end

//-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
//-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-

double alphaL_JC(double *x, double *par)
{// begin

  // Idx = 1 is alpha_1
  // Idx = 2 is alpha_2
  // Idx = 3 is alpha_3 = I_3
  // Idx = 4 is I_12

  double A = 0.0;
  double t = x[0];

  double za  = par[0];
  double Enu = par[1];
  double t1  = par[2];
  int    idx = par[3];

  double th12 = theta_12*c_pi/180.;
  double th13 = theta_13*c_pi/180.;
  double s12  = sin(th12);
  double c12  = cos(th12);
  double s13  = sin(th13);
  double c13  = cos(th13);

  double D21 = dmsq21/(2*Enu);
  double D31 = dmsq31/(2*Enu);

  TF1 *llVt = new TF1("llVt", Vpotential_JC,-tf_max,tf_max,1);
  llVt->SetParameter(0,za);

  TF1 *llDm21t = new TF1("llDm21t",Deltam21_JC,-tf_max,tf_max,2); 
  llDm21t->SetParameter(0,za);
  llDm21t->SetParameter(1,Enu);

  TF1 *llDm32t = new TF1("llDm32t",Deltam32_JC,-tf_max,tf_max,2); 
  llDm32t->SetParameter(0,za);
  llDm32t->SetParameter(1,Enu);

  //- Calculate area under step function
  //** Modified by M.A.AceroO. (9.Apr.2011)
  //** idx = 3 selects alpha2 for the approximated total eigenvalue.

  double ZZ = 0.;
  for (int i=1; i<2*Nshells;i++){
    double tc = times[2*Nshells-1 -  i  ];
    double tn = times[2*Nshells-1 -(i-1)];
    if (t1<tc){
      if (idx == 1)
	ZZ = 0.5*(D21+llVt->Eval(0.5*(tc+tn))*c13**2-llDm21t->Eval(0.5*(tc+tn)));
      else if (idx == 2) 
	ZZ = 0.5*(D21+llVt->Eval(0.5*(tc+tn))*c13**2+llDm21t->Eval(0.5*(tc+tn)));
      else if (idx == 3) 
	ZZ = 0.5*(D31+llVt->Eval(0.5*(tc+tn))*s13**2+D21*c12**2 + llDm21t->Eval(0.5*(tc+tn))-llDm32t->Eval(0.5*(tc+tn)));
      else if (idx == 4)
	ZZ = 0.5*(D21+llVt->Eval(0.5*(tc+tn))*c13**2);
      A += ZZ*(tn-tc)*c_light/hbarc;
    }
    else if (t1<tn){
      if (idx == 1)
	ZZ = 0.5*(D21+llVt->Eval(0.5*(t1+tn))*c13**2-llDm21t->Eval(0.5*(t1+tn)));
      else if (idx == 2)
	ZZ = 0.5*(D21+llVt->Eval(0.5*(t1+tn))*c13**2+llDm21t->Eval(0.5*(t1+tn)));
      else if (idx == 3)
	ZZ = 0.5*(D31+llVt->Eval(0.5*(tc+tn))*s13**2+D21*c12**2 + llDm21t->Eval(0.5*(tc+tn))-llDm32t->Eval(0.5*(tc+tn)));
      else if (idx == 4)
	ZZ = 0.5*(D21+llVt->Eval(0.5*(t1+tn))*c13**2);
      A += ZZ*(tn-t1)*c_light/hbarc;
    }
    if (t<tc){
      if (idx == 1)
	ZZ = 0.5*(D21+llVt->Eval(0.5*(tc+tn))*c13**2-llDm21t->Eval(0.5*(tc+tn)));
      else if (idx == 2) 
	ZZ = 0.5*(D21+llVt->Eval(0.5*(tc+tn))*c13**2+llDm21t->Eval(0.5*(tc+tn)));
      else if (idx == 3) 
	ZZ = 0.5*(D31+llVt->Eval(0.5*(tc+tn))*s13**2+D21*c12**2 + llDm21t->Eval(0.5*(tc+tn))-llDm32t->Eval(0.5*(tc+tn)));
      else if (idx == 4)
	ZZ = 0.5*(D21+llVt->Eval(0.5*(tc+tn))*c13**2);
      A -= ZZ*(tn-tc)*c_light/hbarc;
    }
    else if (t<tn){
      if (idx == 1)
	ZZ = 0.5*(D21+llVt->Eval(0.5*(t+tn))*c13**2-llDm21t->Eval(0.5*(t+tn)));
      else if (idx == 2) 
	ZZ = 0.5*(D21+llVt->Eval(0.5*(t+tn))*c13**2+llDm21t->Eval(0.5*(t+tn)));
      else if (idx == 3) 
	ZZ = 0.5*(D31+llVt->Eval(0.5*(tc+tn))*s13**2+D21*c12**2 + llDm21t->Eval(0.5*(tc+tn))-llDm32t->Eval(0.5*(tc+tn)));
      else if (idx == 4)
	ZZ = 0.5*(D21+llVt->Eval(0.5*(t+tn))*c13**2);
      A -= ZZ*(tn-t)*c_light/hbarc;
    }
  }
  delete llVt;
  delete llDm21t;
  delete llDm32t;
  return A;


}// end

//-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
//-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-

double alphaH_JC(double *x, double *par)
{// begin

  // Idx = 1 is alpha_1 = I_1
  // Idx = 2 is alpha_2
  // Idx = 3 is alpha_3
  // Idx = 4 is I_23

  double A = 0.0;
  double t = x[0];

  double za  = par[0];
  double Enu = par[1];
  double t1  = par[2];
  int    idx = par[3];

  double th12 = theta_12*c_pi/180.;
  double s12  = sin(th12);
  double c12  = cos(th12);

  double D21 = dmsq21/(2*Enu);
  double D31 = dmsq31/(2*Enu);

  TF1 *ffVt = new TF1("ffVt", Vpotential_JC,-tf_max,tf_max,1);
  ffVt->SetParameter(0,za);

  TF1 *ffDm32t = new TF1("ffDm32t",Deltam32_JC,-tf_max,tf_max,2); 
  ffDm32t->SetParameter(0,za);
  ffDm32t->SetParameter(1,Enu);

  //- Calculate area under step function

  double ZZ = 0.;
  for (int i=1; i<2*Nshells;i++){
    double tc = times[2*Nshells-1 -  i  ];
    double tn = times[2*Nshells-1 -(i-1)];
    if (t1<tc){
      if (idx == 1)
	ZZ = D21*c12**2;
      else if (idx == 2) 
	ZZ = 0.5*(ffVt->Eval(0.5*(tc+tn))+D31+D21*s12**2-ffDm32t->Eval(0.5*(tc+tn)));
      else if (idx == 3) 
	ZZ = 0.5*(ffVt->Eval(0.5*(tc+tn))+D31+D21*s12**2+ffDm32t->Eval(0.5*(tc+tn)));
      else if (idx == 4)
	ZZ = 0.5*(ffVt->Eval(0.5*(tc+tn))+D31+D21*s12**2);
      A += ZZ*(tn-tc)*c_light/hbarc;
    }    
    else if (t1<tn){
      if (idx == 1)
	ZZ = D21*c12**2;
      else if (idx == 2)
	ZZ = 0.5*(ffVt->Eval(0.5*(t1+tn))+D31+D21*s12**2-ffDm32t->Eval(0.5*(t1+tn)));
      else if (idx == 3)
	ZZ = 0.5*(ffVt->Eval(0.5*(t1+tn))+D31+D21*s12**2+ffDm32t->Eval(0.5*(t1+tn)));
      else if (idx == 4)
	ZZ = 0.5*(ffVt->Eval(0.5*(t1+tn))+D31+D21*s12**2);
      A += ZZ*(tn-t1)*c_light/hbarc;
    }    
    if (t<tc){
      if (idx == 1)
	ZZ = D21*c12**2;
      else if (idx == 2)
	ZZ = 0.5*(ffVt->Eval(0.5*(tc+tn))+D31+D21*s12**2-ffDm32t->Eval(0.5*(tc+tn)));
      else if (idx == 3)
	ZZ = 0.5*(ffVt->Eval(0.5*(tc+tn))+D31+D21*s12**2+ffDm32t->Eval(0.5*(tc+tn)));    
      else if (idx == 4)
	ZZ = 0.5*(ffVt->Eval(0.5*(tc+tn))+D31+D21*s12**2);
      A -= ZZ*(tn-tc)*c_light/hbarc;
    }
    else if (t<tn){
      if (idx == 1)
	ZZ = D21*c12**2;
      else if (idx == 2)
	ZZ = 0.5*(ffVt->Eval(0.5*(t+tn))+D31+D21*s12**2-ffDm32t->Eval(0.5*(t+tn)));
      else if (idx == 3)
	ZZ = 0.5*(ffVt->Eval(0.5*(t+tn))+D31+D21*s12**2+ffDm32t->Eval(0.5*(t+tn)));
      else if (idx == 4)
	ZZ = 0.5*(ffVt->Eval(0.5*(t+tn))+D31+D21*s12**2);
      A -= ZZ*(tn-t)*c_light/hbarc;
    }
  }
  delete ffVt;
  delete ffDm32t;
  return A;


}// end

//-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
//-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-

double Deltam21_JC(double *x, double *par)
{// begin

  double D = 0.0;
  double t = x[0];

  double za  = par[0];
  double Enu = par[1];

  double th13 = theta_13*c_pi/180.;
  double c13  = cos(th13);

  TF1 *lVt = new TF1("lVt", Vpotential_JC,-tf_max,tf_max,1);
  lVt->SetParameter(0,za);
  double Vt = lVt->Eval(t);

  TF1 *llVl = new TF1("llVl", VlowRes_JC,0.0,Enu_max,0);
  double Vl = llVl->Eval(Enu);

  double Bl = Vl*tan(2* theta_12*c_pi/180.);

  D = c13**2*sqrt( (Vt-Vl)**2 + Bl**2 );
delete lVt;
delete llVl;
  return D;
  
}//end

double Deltam32_JC(double *x, double *par)
{// begin

  double D = 0.0;
  double t = x[0];

  double za  = par[0];
  double Enu = par[1];
  

  TF1 *fVt = new TF1("fVt", Vpotential_JC,-tf_max,tf_max,1);
  fVt->SetParameter(0,za);
  double Vt = fVt->Eval(t);

  TF1 *ffVh = new TF1("ffVh", VhighRes_JC,0.0,Enu_max,0);
  double Vh = ffVh->Eval(Enu);

  double Bh = Vh*tan(2* theta_13*c_pi/180.);

  D = sqrt( (Vt-Vh)**2 + Bh**2 );
delete fVt;
delete ffVh;
  return D;
  
}//end

double Umatt_JC(double *x, double *par)
{// begin

  bool Fix_c12m = false;

  double U = 0.0;
  double t = x[0];

  double za   = par[0];
  double Enu  = par[1];
  int  idx_f  = par[2]; // 0  e, 1 mu, 2 tau
  int  idx_m  = par[3]; // 0 m1, 1 m2, 2  m3
  int  reim   = par[4]; // 0 Re part, 1 Im part, 2 Im part of cc

  double th23 = theta_23*c_pi/180.;
  double s23  = sin(th23);
  double c23  = cos(th23);

  double thm12 = c_pi/2;
  if (!Fix_c12m){
    TF1 *ffThm12t = new TF1("ffThm12t",Thetam12_JC,-tf_max,tf_max,2);
    ffThm12t->SetParameter(0,za); 
    ffThm12t->SetParameter(1,Enu); 
    thm12 = ffThm12t->Eval(t)*c_pi/180.;
    delete ffThm12t;
  }
  
  TF1 *ffThm13t = new TF1("ffThm13t",Thetam13_JC,-tf_max,tf_max,2);
  ffThm13t->SetParameter(0,za); 
  ffThm13t->SetParameter(1,Enu); 
  double thm13 = ffThm13t->Eval(t)*c_pi/180.;
  delete ffThm13t;

  double c12m = cos(thm12);
  double s12m = sin(thm12);
  double c13m = cos(thm13);
  double s13m = sin(thm13);
  double d    = deltaCP*c_pi/180.;

  // The Mixing matrix in matter:

  if      (idx_f == 0){// electron row
    
    if      (idx_m == 0){
      if      (reim == 0) U = c12m*c13m;
      else if (reim == 1) U = 0.0;
      else if (reim == 2) U = 0.0;
    }
    else if (idx_m == 1){
      if      (reim == 0) U = s12m*c13m;
      else if (reim == 1) U = 0.0;
      else if (reim == 2) U = 0.0;
    }
    else if (idx_m == 2){
      if      (reim == 0) U = s13m*cos(d);
      else if (reim == 1) U =-s13m*sin(d);
      else if (reim == 2) U = s13m*sin(d);
    }
    
  }
  else if (idx_f == 1){// muon row
    
    if      (idx_m == 0){
      if      (reim == 0) U = -s12m*c23-c12m*s23*s13m*cos(d);
      else if (reim == 1) U = -c12m*s23*s13m*sin(d);
      else if (reim == 2) U =  c12m*s23*s13m*sin(d);
    }
    else if (idx_m == 1){
      if      (reim == 0) U = c12m*c23-s12m*s23*s13m*cos(d);
      else if (reim == 1) U = -s12m*s23*s13m*sin(d);
      else if (reim == 2) U =  s12m*s23*s13m*sin(d);
    }
    else if (idx_m == 2){
      if      (reim == 0) U = s23*c13m;
      else if (reim == 1) U = 0.0;
      else if (reim == 2) U = 0.0;
    }

  }
  else if (idx_f == 2){// tau row
        
    if      (idx_m == 0){
      if      (reim == 0) U = s12m*s23-c12m*c23*s13m*cos(d);
      else if (reim == 1) U = -c12m*c23*s13m*sin(d);
      else if (reim == 2) U =  c12m*c23*s13m*sin(d);
    }
    else if (idx_m == 1){
      if      (reim == 0) U = -c12m*s23-s12m*c23*s13m*cos(d);
      else if (reim == 1) U = -s12m*c23*s13m*sin(d);
      else if (reim == 2) U =  s12m*c23*s13m*sin(d);
    }
    else if (idx_m == 2){
      if      (reim == 0) U = c23*c13m;
      else if (reim == 1) U = 0.0;
      else if (reim == 2) U = 0.0;
    }
  }
    
  return U;

}// end

//-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
//-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-

double Thetam12_JC(double *x, double *par)
{// begin

  double T = 0.0;
  double t = x[0];

  double za  = par[0];
  double Enu = par[1];

  TF1 *gfVt = new TF1("gfVt", Vpotential_JC,-tf_max,tf_max,1);
  gfVt->SetParameter(0,za);
  double Vt = gfVt->Eval(t);
  
  TF1 *fVl = new TF1("fVl", VlowRes_JC,0.0,Enu_max,0);
  double Vl = fVl->Eval(Enu);
  
  double Bl = Vl*tan(2* theta_12*c_pi/180.);

  T = 0.5*acos( (Vl-Vt) / sqrt( (Vl-Vt)**2 + Bl**2 ) );
  T = T*180.0/c_pi;
  delete gfVt;
  delete fVl;

  return T;

}//end

//-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
//-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-

double Thetam13_JC(double *x, double *par)
{// begin

  double T = 0.0;
  double t = x[0];

  double za  = par[0];
  double Enu = par[1];

  TF1 *fVt = new TF1("fVt", Vpotential_JC,-tf_max,tf_max,1);
  fVt->SetParameter(0,za);
  double Vt = fVt->Eval(t);
  delete fVt;
  
  TF1 *fVh = new TF1("fVh", VhighRes_JC,0.0,Enu_max,0);
  double Vh = fVh->Eval(Enu);
  delete fVh;
  
  double Bh = Vh*tan(2* theta_13*c_pi/180.);

  T = 0.5*acos( (Vh-Vt) / sqrt( (Vh-Vt)**2 + Bh**2 ) );
  T = T*180.0/c_pi;

  return T;

}//end

//-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
//-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-

double VlowRes_JC(double *x, double *par)
{//begin

  double E = x[0];

  //printf("Energy = %f \n",E);

  double th12 = theta_12*c_pi/180.;
  double th13 = theta_13*c_pi/180.;
  double c13  = cos(th13);
  //printf("cos(th13) = %f \n",c13);

  double V = (dmsq21/(2.0*E))*cos(2.0*th12)/c13**2;
  //printf("VlowR = %f \n",V);

  return V;

}//end

//-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
//-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-

double VhighRes_JC(double *x, double *par)
{//begin

  double E = x[0];

  double th12 = theta_12*c_pi/180.;
  double th13 = theta_13*c_pi/180.;
  double s12  = sin(th12);

  double V = (dmsq31-dmsq21*s12**2)*cos(2*th13)/(2*E);

  return V;

}//end

//-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
//-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-

double Vpotential_JC(double *x, double *par)
{// begin

  double V;
  double t = x[0];

  double za= par[0];

  //- Cet array of crossing times for the zenith angle
  GetCrossingTimes(za);
  

  V = 0.;
  for (int i=Nshells-1; i>=0; i--)
    if      (fabs(t)<= tcross[i]) V = V_shell[i];

  
  return V;
  
}//end

//-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
//-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-

void GetCrossingTimes(double za)
{// begin

  if (first)
    printf("Making array of crossing times for za=%4.1f deg \n",za);

   //- calculate path lengths from the zenith angle
  za = za*c_pi/180.;
  for (int i = 0; i<Nshells; i++){
    l_ss[i] = 0.0;
    if (sin(za)<R_shell[i]/R_Earth) 
      l_ss[i] = 2*R_shell[i]*
	sqrt(1-(R_Earth*sin(za)/R_shell[i])**2);
  }// for i
  
  //- calculate crossing times for each shell from the middle, out
  for (int i = 0; i<Nshells;i++) 
    tcross[i] = 0.5*l_ss[i]/c_light;
  
  //- construct array of crossing times from entry to exit 
  //- (zero is at midpoint)
  for (int i = 0; i<Nshells; i++){
    times[i] = -tcross[Nshells-1-i];
    times[Nshells+i] = tcross[i];
  }

  for (int i = 0; i<2*Nshells; i++)
    if (first) 
      printf("times[%d] = %9.6f\n",i,times[i]);  

  tlohi[0]= -tcross[Nshells-1];
  tlohi[1]=  tcross[Nshells-1]; 
  
  if (first) 
    printf("------------------------------------------ Done! \n");
  

  first = false;
  return;
  
}// end

//-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
//-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-