#define NDEBUG
#include <assert.h>

//#include <time.h>

using namespace std;
#include "re_tab.hxx"

ClassImp(TABLE)
ClassImp(RE_TAB)

int RE_TAB::ZA[92]={
1,4,7,9,11,12,14,16,19,20,23,24,27,28,31,32,35,40,39,40,45,48,51,52,55,56,
59,58,63,64,69,74,75,80,79,84,85,88,90,90,93,96,99,102,103,106,107,114,115,
120,121,130,127,130,134,138,139,140,141,142,147,152,153,158,160,164,165,166,
169,174,175,180,181,184,187,192,193,195,197,202,205,208,209,210,210,222,223,
226,227,232,231,238};

//______________________________________________________________________________
RE_TAB::RE_TAB(const char *mat_prefix, 
               const char *tab_source,
               const char *tab_dir)
              {

  /// Default location for energy calibration tables:
  TString  sReTabDirLocalDefPath  = "calfiles/KRATTA/Energy_calibration/eloss";
  
  g_mat_prefix         = mat_prefix; 
  g_tab_source         = tab_source;
  g_tab_dir            = tab_dir;

  TString  sAsyeosrootMainDir     = gSystem->ExpandPathName("$VMCWORKDIR");
  g_tab_dir_default_path  = sAsyeosrootMainDir + "/" + sReTabDirLocalDefPath;
  g_az_spec            = Form("%s/tabs_atima/az.dat", g_tab_dir_default_path.Data() );
  
  MPART[0]             = MEl/AMU; // Electron mass (amu), use Z=1, A=-1
  MPART[1]             = MEl/AMU; // Electron mass (amu), use Z=1, A=-2
  MPART[2]             = MPi/AMU; // Pion     mass (amu), use Z=1, A=-3
  MPART[3]             = MMu/AMU; // Muon     mass (amu), use Z=1, A=-4
  MPART[4]             = MK/AMU;  // Kaon     mass (amu), use Z=1, A=-5

  INIT_APZP();
}

//______________________________________________________________________________
RE_TAB::~RE_TAB(){}

//______________________________________________________________________________
void RE_TAB::INIT_APZP(){

  int i;

  //printf("TABS for: %s\n",mat_prefix);

  for(i=0;i<NAZ;i++){
    A[i]=0;
    Z[i]=0;
    M[i]=0;
    sprintf(SYM[i],"  ");
    TAB[i].AZ=0;
    TAB[i].Z=0;
    TAB[i].A=0;
    TAB[i].M=0;
    sprintf(TAB[i].SYM,"  ");
  }
  NAZMAXI=0;
  NAZMAXI_AZ=0;

  string dummy;

  ifstream az_file(g_az_spec,ios::in);
  if(az_file.fail()){
    printf("[RE_TAB:] >>> %s not found\n",g_az_spec);
    exit(1);
  }
  getline(az_file,dummy);

  char sym[3];
  int a,z,naz=0;
  float m,amab,rho;
  while(az_file >> z >> sym >> a >> m >> amab >> rho){
    A[naz]    = a;
    Z[naz]    = z;
    M[naz]    = m;
    sprintf(SYM[naz],"%s",&sym[0]);
    naz++;
  }
  az_file.close();
  NAZMAXI_AZ=naz;
  if(NAZMAXI_AZ>=NAZ){
    printf("[RE_TAB:] increase NAZ dimension\n");
    exit(1);
  }
}
//______________________________________________________________________________

float RE_TAB::RANGE(int Zp, int Ap, float energy){

  float range=-1;
                  
  int naz=NAPZP(Zp,Ap);             // Zp,Ap requested species

  if(!TAB[naz].AZ) INIT_TABLE(naz); // initialize table for the nearest A and Z

  if((Ap==TAB[naz].A && Zp==TAB[naz].Z) || Ap==0){ //Zp,Ap match the ones in the table
    range = TAB[naz].AZ->RANGE(energy);
  }
  else{                                 // have to scale
    float M1 = TAB[naz].M;
    float Z1 = TAB[naz].Z;
    float M2 = float(Ap);
    float Z2 = float(Zp);
    if(Ap<0 && Zp==1 && TAB[naz].A==1 && Z1==1) M2 = MPART[-Ap-1];//Pi, Mu, K etc.
    range  = (M2*Z1*Z1)/(M1*Z2*Z2)*TAB[naz].AZ->RANGE(energy*M1/M2);
  }

  return range;
}

//______________________________________________________________________________

float RE_TAB::ENERGY(int Zp, int Ap, float range){
  
  float energy=-1;

  int naz=NAPZP(Zp,Ap);             // Zp,Ap requested species

  if(!TAB[naz].AZ) INIT_TABLE(naz); // initialize table for the nearest A and Z

  if((Ap==TAB[naz].A && Zp==TAB[naz].Z) || Ap==0){ //Zp,Ap match the ones in the table
    energy = TAB[naz].AZ->ENERGY(range);
  }
  else{                                 // have to scale
    float M1 = TAB[naz].M;
    float Z1 = TAB[naz].Z;
    float M2 = float(Ap);
    float Z2 = float(Zp);
    if(Ap<0 && Zp==1 && TAB[naz].A==1 && Z1==1) M2 = MPART[-Ap-1];//Pi, Mu, K etc.
    energy = M2/M1*TAB[naz].AZ->ENERGY(range*(M1*Z2*Z2)/(M2*Z1*Z1));
  }

  return energy;
}

// NOTE: THIS function: RE_TAB::INIT_TABLE is localy dependeed! 
// Correct that!
//______________________________________________________________________________
// LUKASIK VERSION

//void RE_TAB::INIT_TABLE(int naz){

//char in_fname[100];
//char in_dname[100];
//if(strstr(g_tab_dir,"default")){
  //if(strstr(g_tab_source,"ATIMA_OLD")) 
    //sprintf(in_dname,"/u/lukasik/eloss/tabs_atima_old");
  //else if(strstr(g_tab_source,"ATIMA")) 
    //sprintf(in_dname,"/u/lukasik/eloss/tabs_atima");
  //else if(strstr(g_tab_source,"PRAL-ZIEGLER")) 
    //sprintf(in_dname,"/u/lukasik/eloss/tabs_pral_ziegler");
  //else if(strstr(g_tab_source,"PRAL-SRIM")) 
    //sprintf(in_dname,"/u/lukasik/eloss/tabs_pral_srim");
  //else if(strstr(g_tab_source,"PRAL")) 
    //sprintf(in_dname,"/u/lukasik/eloss/tabs_pral");
  //else if(strstr(g_tab_source,"VEDA")) 
    //sprintf(in_dname,"/u/lukasik/eloss/tabs_veda");
  //else if(strstr(g_tab_source,"ZIEGLER")) 
    //sprintf(in_dname,"/u/lukasik/eloss/tabs_ziegler");
  //else if(strstr(g_tab_source,"SRIM")) 
    //sprintf(in_dname,"/u/lukasik/eloss/tabs_srim");
  //else{
    //sprintf(in_dname,"unknown");
    //printf("available tables: ATIMA | PRAL | PRAL-PROJ | VEDA | ZIEGLER | SRIM\n");
  //}
  //sprintf(in_fname,"%s/%s_%s_%d",in_dname,g_mat_prefix,&SYM[naz][0],A[naz]);
//}
//else
//sprintf(in_fname,"%s/%s_%s_%d",g_tab_dir,g_mat_prefix,&SYM[naz][0],A[naz]);


//TAB[naz].AZ = new TABLE(&in_fname[0],g_tab_source,
                        //&SYM[naz][0],Z[naz],A[naz],M[naz]);
//TAB[naz].Z  = Z[naz];
//TAB[naz].A  = A[naz];
//TAB[naz].M  = M[naz];
//sprintf(TAB[naz].SYM,"%s",&SYM[naz][0]);
//NAZMAXI++;
////printf(">>> %3d out of %3d (%7.3f%%) tables in use : %s\tloaded\n",
////       NAZMAXI,NAZMAXI_AZ,NAZMAXI*100./float(NAZMAXI_AZ),in_fname);
//if(NAZMAXI>=NAZ){printf("NAZMAXI>=NAZ, RE_TAB::INIT_TABLE\n");exit(1);}
//}

////______________________________________________________________________________
void RE_TAB::INIT_TABLE(int naz){

  char in_fname[200];

  TString re_tab_dir_source ="";

  if(strstr( g_tab_dir, "default" )){
    
    if     (strstr(g_tab_source,"ATIMA_OLD"))    re_tab_dir_source = "tabs_atima_old";
    else if(strstr(g_tab_source,"ATIMA"))        re_tab_dir_source = "tabs_atima";
    else if(strstr(g_tab_source,"PRAL-ZIEGLER")) re_tab_dir_source = "tabs_pral_ziegler";
    else if(strstr(g_tab_source,"PRAL-SRIM"))    re_tab_dir_source = "tabs_pral_srim";
    else if(strstr(g_tab_source,"PRAL"))         re_tab_dir_source = "tabs_pral";
    else if(strstr(g_tab_source,"VEDA"))         re_tab_dir_source = "tabs_veda";
    else if(strstr(g_tab_source,"ZIEGLER"))      re_tab_dir_source = "tabs_ziegler";
    else if(strstr(g_tab_source,"SRIM"))         re_tab_dir_source = "tabs_srim";
    else{                                        re_tab_dir_source = "unknown";
      printf("[RE_TAB:] Available tables: ATIMA | PRAL | PRAL-PROJ | VEDA | ZIEGLER | SRIM\n");
    }
    TString in_dname = g_tab_dir_default_path + "/" + re_tab_dir_source;
    sprintf(in_fname,"%s/%s_%s_%d",in_dname.Data() ,g_mat_prefix,&SYM[naz][0],A[naz]);
  }else{
     sprintf(in_fname,"%s/%s_%s_%d",g_tab_dir,g_mat_prefix,&SYM[naz][0],A[naz]);
  }
  
  TAB[naz].AZ = new TABLE(&in_fname[0],g_tab_source,
                          &SYM[naz][0],Z[naz],A[naz],M[naz]);
  TAB[naz].Z  = Z[naz];
  TAB[naz].A  = A[naz];
  TAB[naz].M  = M[naz];
  sprintf(TAB[naz].SYM,"%s",&SYM[naz][0]);
  NAZMAXI++;
  //printf(">>> %3d out of %3d (%7.3f%%) tables in use : %s\tloaded\n",
  //       NAZMAXI,NAZMAXI_AZ,NAZMAXI*100./float(NAZMAXI_AZ),in_fname);
  if(NAZMAXI>=NAZ){printf("[RE_TAB:]  NAZMAXI>=NAZ, RE_TAB::INIT_TABLE\n");exit(1);}
}  

int RE_TAB::GetNSPMAXI(int Zp, int Ap){
int naz=NAPZP(Zp,Ap);
if(!TAB[naz].AZ) INIT_TABLE(naz);
return TAB[naz].AZ->GetNSPMAXI();
}
//______________________________________________________________________________

float RE_TAB::GetR(int Zp, int Ap, int i){
int naz=NAPZP(Zp,Ap);
if(!TAB[naz].AZ) INIT_TABLE(naz);
return TAB[naz].AZ->GetR(i);
}
//______________________________________________________________________________

float RE_TAB::GetE(int Zp, int Ap, int i){
int naz=NAPZP(Zp,Ap);
if(!TAB[naz].AZ) INIT_TABLE(naz);
return TAB[naz].AZ->GetE(i);
}
//______________________________________________________________________________

float* RE_TAB::GetR(int Zp, int Ap){
int naz=NAPZP(Zp,Ap);
if(!TAB[naz].AZ) INIT_TABLE(naz);
return TAB[naz].AZ->GetR();
}
//______________________________________________________________________________

float* RE_TAB::GetE(int Zp, int Ap){
int naz=NAPZP(Zp,Ap);
if(!TAB[naz].AZ) INIT_TABLE(naz);
return TAB[naz].AZ->GetE();
}
//______________________________________________________________________________

int RE_TAB::GetA(int Zp, int Ap){
int naz=NAPZP(Zp,Ap);
if(!TAB[naz].AZ) INIT_TABLE(naz);
return TAB[naz].A;
}
//______________________________________________________________________________

int    RE_TAB::GetZ(int Zp, int Ap){
int naz=NAPZP(Zp,Ap);
if(!TAB[naz].AZ) INIT_TABLE(naz);
return TAB[naz].Z;
}
//______________________________________________________________________________

float  RE_TAB::GetM(int Zp, int Ap){
int naz=NAPZP(Zp,Ap);
float Mp = 0;
if(!TAB[naz].AZ) INIT_TABLE(naz);
if((Ap==TAB[naz].A && Zp==TAB[naz].Z) || Ap==0){ //Zp,Ap match the ones in the table
  Mp = TAB[naz].M;
}
else{                                 // have to scale
  Mp = float(Ap);
}
if(Ap<0 && Zp==1 && TAB[naz].A==1 && TAB[naz].Z==1) Mp = MPART[-Ap-1];//Pi, Mu, K etc.
  return Mp;
}
//______________________________________________________________________________

char*  RE_TAB::GetSym(int Zp, int Ap){
int naz=NAPZP(Zp,Ap);
if(!TAB[naz].AZ) INIT_TABLE(naz);
return &TAB[naz].SYM[0];
}
//______________________________________________________________________________

int RE_TAB::NAPZP(int Zp, int Ap){
  // search for nearest A and Z present in the az.dat file
  // returns consecutive number of
  // apropriate table[naz] for requested Ap, Zp
  // and actual Atab and Ztab for this table[naz]

  // Zp,Ap requested species
  // Ztab  nearest species present in the tables

  int Ztab=-1;
  int i, naz=0;
  int dz=1000;
  int da=1000;

  if(Ap==0) Ap=ZA[(Zp<=92) ? Zp-1 : 91];
  for(i=0;i<NAZMAXI_AZ;i++){
    if(abs(Zp-Z[i])<=dz){
      dz   = abs(Zp-Z[i]);
      Ztab = Z[i];
    }
  }
  for(i=0;i<NAZMAXI_AZ;i++){
    if(Z[i]==Ztab && abs(Ap-A[i])<=da){
      da   = abs(Ap-A[i]);
      naz  = i;
      //Atab = A[i];
    }
  }
  if(Ztab==-1){printf("[RE_TAB:] Error in RE_TAB::NAPZP %d %d\n",Zp,Ap); exit(1);}
  return naz;
}

//______________________________________________________________________________





// RANGE-ENERGY TABLE 



//______________________________________________________________________________

TABLE::TABLE( char *tab_fname,
        const char *tab_source,
              char *SPEC, 
	       int  Z, 
	       int  A, 
	     float  M){

READ_TABLE(tab_fname,tab_source,SPEC,Z,A,M);
DERIVXY();
}
//______________________________________________________________________________

TABLE::~TABLE(){}
//______________________________________________________________________________

void TABLE::READ_TABLE(      char *tab_fname,
                       const char *tab_source,
                             char *SPEC, 
		              int  Z, 
		              int  A, 
		            float  M){

int i;

//printf("TABS for: %s\n",mat_prefix);

for(i=0;i<DIM;i++) E[i]=R[i]=S[i]=0;


ifstream in_file(tab_fname,ios::in);
if(in_file.fail()){
  printf("[RE_TAB:] >>> %s not found\n",tab_fname);
  exit(1);
}

if(strstr(tab_source,"ATIMA_OLD")){
  string dummy;

  int a;
  char dum1[20], dum2[20],spec[3];
  float e, dee, r, strag;


  for(i=0;i<10;i++){
    getline(in_file,dummy);
  }
  in_file >> dum1 >> spec >> a >> dum1 >> dum2;
  if(!strstr(spec,SPEC) || a!=A){
    printf("[RE_TAB:] conflict !!!\n");
    exit(1);
  }

  assert(a!=A); // assert test

  for(i=0;i<2;i++){
    getline(in_file,dummy);
  }
  E[0] = 0.;
  R[0] = 0.;
  i=1;
  while(in_file >> e >> r >> strag >> dee){
    E[i] = e*M;
    R[i] = r;
    i++;
    if(i==DIM){
      printf("[RE_TAB:] increase table DIMension\n");
      exit(1);
    }
  }
  NSPMAXI = i;
}
else if(strstr(tab_source,"ATIMA")){
  string dummy;

  int a, z;
  float m;
  char spec[3],ch;
  float e, dee, r, strag;


  for(i=0;i<2;i++){
    getline(in_file,dummy);
  }
//  in_file >>  spec >> z >> a >> m >> dum1 >> dum2 >> dum1 >> dum2;
  in_file >>  spec >> z >> a >> m;
  getline(in_file,dummy); // get the rest
 // printf("%s %d %d %f %s %s\n",spec,z,a,m,dum1,dum2);
  if(!strstr(spec,SPEC) || a!=A){
    printf("[RE_TAB:] Conflict !!!\n");
    exit(1);
  }

  assert(a!=A); // assert test

  for(i=0;i<5;i++){
    getline(in_file,dummy);
  //printf("%d %s\n",i,(char*)dummy.c_str());
  }
  E[0] = 0.;
  R[0] = 0.;
  i=1;
  while(in_file >> e >> ch >> dee >> ch >> r >> ch >> strag >> ch){
    E[i] = e*M;
    R[i] = r;
  //printf("%f %f %f %f\n",e,dee,r,strag);
    i++;
    if(i==DIM){
      printf("[RE_TAB:] Increase table DIMension\n");
      exit(1);
    }
  }
  NSPMAXI = i;
}
else if(strstr(tab_source,"PRAL-ZIEGLER")|| strstr(tab_source,"PRAL-SRIM")){
  char CH;
  string dummy;
  float e, r, Ap, Zp, At, Zt;
  float r_proj,SIGMAX,SIGMAZ,SE,SN,QN;
  for(i=0;i<2;i++){
    getline(in_file,dummy);
    //in_file.get(CH);
  }
  in_file >> Ap >> Zp >> At >> Zt;
  in_file.get(CH);
  getline(in_file,dummy);
  //in_file.get(CH);
  getline(in_file,dummy);
  //in_file.get(CH);

  if(Zp!=Z || Ap!=A){
    printf("[RE_TAB:] Conflict !!!\n");
    exit(1);
  }

  E[0] = 0.;
  R[0] = 0.;
  i=1;
  while(in_file >> e >> r >> r_proj >> SIGMAX >> SIGMAZ >> SE >> SN >> QN){
    E[i] = e*Ap;
    R[i] = r;
    //printf("%d %f MeV/n %f um\n",i,e,r);
    i++;
    if(i==DIM){
      printf("[RE_TAB:] Increase table DIMension\n");
      exit(1);
    }
  }
  NSPMAXI = i;
}
else if(strstr(tab_source,"PRAL")){
  char CH;
  string dummy;
  float e, r, Ap, Zp, At, Zt;
  float r_proj,SIGMAX,SIGMAZ,SE,SN,QN;
  bool pral_proj = strstr(tab_source,"PRAL-PROJ");
  //printf("%s %d\n",tab_source,pral_proj);
  for(i=0;i<2;i++){
    getline(in_file,dummy);
    //in_file.get(CH);
  }
  in_file >> Ap >> Zp >> At >> Zt;
  in_file.get(CH);
  getline(in_file,dummy);
  //in_file.get(CH);
  getline(in_file,dummy);
  //in_file.get(CH);

  if(Zp!=Z || Ap!=A){
    printf("[RE_TAB:] Conflict !!!\n");
    exit(1);
  }

  E[0] = 0.;
  R[0] = 0.;
  i=1;
  while(in_file >> e >> r >> r_proj >> SIGMAX >> SIGMAZ >> SE >> SN >> QN){
    E[i] = e*Ap;
    R[i] = (pral_proj) ? r_proj : r;
    //printf("%d %f MeV/n %f um\n",i,e,r);
    i++;
    if(i==DIM){
      printf("[RE_TAB:] Increase table DIMension\n");
      exit(1);
    }
  }
  NSPMAXI = i;
}
else if(strstr(tab_source,"VEDA") || strstr(tab_source,"ZIEGLER")){
  char CH;
  string dummy;
  float e, r, Ap, Zp, At, Zt;

  for(i=0;i<2;i++){
    getline(in_file,dummy);
    //in_file.get(CH);
  cout <<  dummy << "\n";
  }
  in_file >> Ap >> Zp >> At >> Zt;
  in_file.get(CH);
  getline(in_file,dummy);
  //in_file.get(CH);
  cout <<  dummy << "\n";
  cout <<  Ap << " "<< Zp << " " << At << " " <<  Zt<< "\n";


  E[0] = 0.;
  R[0] = 0.;
  i=1;
  while(in_file >> e >> r){
    E[i] = e*Ap;
    R[i] = r;
    //printf("%d %f MeV/n %f um\n",i,e,r);
    i++;
    if(i==DIM){
      printf("[RE_TAB:] Increase table DIMension\n");
      exit(1);
    }
  }
  NSPMAXI = i;
}
else if(strstr(tab_source,"SRIM")){
  char CH;
  string dummy;
  char kev[20], um[20], dum1[20], dum2[20];
  float e, dee, den, r, los, las;

  //in_file.setf(ios::skipws);
  while(in_file.get(CH)){
    if(CH!='-'){
      getline(in_file,dummy);
      //in_file.get(CH);
      //printf(">>%s\n",dummy);
    }
    else{
      getline(in_file,dummy);
      //in_file.get(CH);
      //printf(">>%s\n",dummy);
      break;
    }
  }
  E[0] = 0.;
  R[0] = 0.;
  i=1;
  while(in_file.get(CH)){
    //printf("%c\n",CH);
    if(CH!='-'){
      in_file.putback(CH);
      in_file >> e >> kev >> dee >> den >> r >> 
      um >> los >> dum1 >> las >> dum2;
      getline(in_file,dummy);
      //in_file.get(CH);
      //printf("%f %s %e %e %f %s %f %s %f %s\n",
      //e,kev,dee,den,r,um,los,dum1,las,dum2);
      if     (strstr(kev,"MeV")) e  = e;        // MeV -> MeV
      else if(strstr(kev,"keV")) e /= 1000.;    // keV -> MeV
      else if(strstr(kev,"GeV")) e *= 1000.;    // GeV -> MeV
      else if(strstr(kev,"TeV")) e *= 1000000.; // TeV -> MeV
      else if(strstr(kev,"eV"))  e /= 1000000.; // eV  -> MeV
      else{
        printf("[RE_TAB:] Unknown energy units, exiting...\n");
        exit(1);
      }
      if     (strstr(um,"um"))  r  = r;	       // um -> um
      else if(strstr(um,"mm"))  r *= 1000.;    // mm -> um
      else if(strstr(um,"A"))   r /= 10000.;   // A  -> um
      else if(strstr(um,"cm"))  r *= 10000.;   // cm -> um
      else if(strstr(um,"km"))  r *= 1.e9; // km  -> um
      else if(strstr(um,"m"))   r *= 1000000.; // m  -> um
      else{
        printf("[RE_TAB:] Unknown range units, exiting...\n");
        exit(1);
      }
      //printf("%f MeV %f um\n",e,r);
      E[i] = e;
      R[i] = r;
      i++;
      if(i==DIM){
        printf("[RE_TAB:] Increase table DIMension\n");
        exit(1);
      }
    }
    else break;
  }
  NSPMAXI = i;
}
else{
  printf("[RE_TAB:] R-E tables not defined (ATIMA | PRAL | PRAL-PROJ | VEDA | "
  "ZIEGLER | SRIM | PRAL-ZIEGLER | PRAL-SRIM)\n");
  exit(1);
}

in_file.close();
}
//______________________________________________________________________________

void TABLE::DERIVXY(){

//C******************
//C     CALCUL DES DERIVEES SECONDES
//C*****************

int N,N1,N2,N3,NM;
float H1,H2,H3,G1,G2,G3,HH,HHP,W,U,WP,UP;

N=NSPMAXI-1;
N1=N-1;
N2=N-2;
N3=N-3;
H1=E[N]-E[N1];
H2=E[N1]-E[N2];
H3=E[N2]-E[N3];
G1=R[N]-R[N1];
G2=R[N1]-R[N2];
G3=R[N2]-R[N3];
HH=2*H1+3*H2+H2*H2/H1;
HHP=HH;
W=6*(G1/H1-G2/H2)/HH;
U=(H2*H2/H1-H1)/HH;
HHP=2*H2+3*H3+H3*H3/H2;
WP=6*(G2/H2-G3/H3)/HHP;
UP=(H3*H3/H2-H2)/HHP;
S[N]=(W*(H1+H2)-H1*(W*UP+WP))/(H2-U*(H1+H2)+U*UP*H1);
NM=N1;
for(int I=NM;I>0;I--){
  S[I]=U*S[I+1]+W;
  if(I==1) break;
  N1=N2;
  N2=I-2;
  H1=H2;
  H2=E[N1]-E[N2];
  G1=G2;
  G2=R[N1]-R[N2];
  HH=2*H1+3*H2+H2*H2/H1;
  W=6*(G1/H1-G2/H2)/HH;
  U=(H2*H2/H1-H1)/HH;
}
S[0]=(1+H2/H1)*S[1]-H2*S[2]/H1;
return;
      
}
//______________________________________________________________________________

float TABLE::RANGE(float XX){

//C********************
//C     CALCULATION OF THE RANGE FOR A GIVEN ENERGY
//C********************

int N,I,J;
float DG,DD,H,rRANGE;

N=NSPMAXI-1;
for(I=1;I<=N;I++){
  if(XX <= E[I]) break;
}
J=I;
I--;
DG=XX-E[I];
DD=E[J]-XX;
H=E[J]-E[I];
rRANGE=S[I]*DD*DD*DD/(6*H)+S[J]*DG*DG*DG/(6*H)+(R[J]/H-S[J]*H/6)*DG
+(R[I]/H-S[I]*H/6.)*DD;

return rRANGE;
      
}

//______________________________________________________________________________

float TABLE::ENERGY(float YY){

//C********************
//C     CALCULATION OF THE ENERGY FOR A GIVEN RANGE
//C********************

//      IMPLICIT REAL*8 (A-H,O-Z)
const double DEG2RAD = 0.0174533;
const double RAD2DEG = 57.2958;

float rENERGY;
int N,I,J,K;

double RAC[3];
//double RMIN=-.001, RMAX=1.001;
double RMIN=-.005, RMAX=1.005;
double H,H2,B1,B2,B3,B4,A,B,C,A0,A1,A2,A3,DEL,DELSQ,R0,Q,D,RR,DSQ;
double RHO,S1R,S2R,PHY,CPHY,SPHY,S1I;


if(YY < 0.){
  //printf("PROBLEM 1 DANS rENERGY\n");
  //exit(1);
}
N=NSPMAXI-1;
if(YY>R[N]){
  ////printf("PROBLEM 2 DANS rENERGY\n");
  ////printf("%f > %f (R[N])\n",YY,R[N]);
  //rENERGY = 0.;
  //for(int i=0;i<4;i++) rENERGY += COF_ER[i]*pow(YY,i); // pol3
  ////rENERGY = (E[N]-E[N-1])/(R[N]-R[N-1])*(YY-R[N-1])+E[N-1]; // pol1
  //printf("extrapolating %f > R[N]=%f, E=%f %f\n",YY,R[N],rENERGY,
  //(E[N]-E[N-1])/(R[N]-R[N-1])*(YY-R[N-1])+E[N-1]);
  //return rENERGY;
  ////exit(1);
//printf("extrapolating %f > R[N]=%f\n",YY,R[N]);
}
//C     DO 1 I=2,N
//C     if(YY-R[I])2,2,1
for(I=N;I>=0;I--){
  if(YY>R[I]) break;
}
//STOP 'PROBLEM 3 DANS rENERGY'
J=I+1;
//I=I-1;
//printf("SPDEDX %d %f %f %f\n",I,YY,R[J],R[I]);
H=E[J]-E[I];
H2=H*H/6.;
B1=S[I]*H2;
B2=S[J]*H2;
B3=R[J]-S[J]*H2;
B4=R[I]-S[I]*H2;
A=B1-B2;
A3=A;
A2=3.*B2;
A1=(B4-B3-3.*B2);
A0=(B2+B3-YY);
if(fabs(A3)+fabs(A2)+fabs(A1)+fabs(A0) == 0.){
  rENERGY=E[J];
  return rENERGY;
}
if(fabs(A3)<=fabs(0.0001*A2)){
//C------------------
//C         EQUATION DU 2EME DEGRE A*X**2+B*X+C=0
//C------------------
  A=A2;
  B=A1;
  C=A0;
  if(fabs(A)<=fabs(0.0001*B)){
    if(fabs(B)<=fabs(0.0001*C)){
      rENERGY=E[J];
      return rENERGY;
    }
    RAC[0]=-C/B;
    if(RAC[0]>=RMIN && RAC[0]<=RMAX){
      rENERGY=E[J]-H*RAC[0];
      return rENERGY;
    }
    else{
      printf("[RE_TAB:] PROBLEME 1 DANS ENERGY(Y) >> %f\n",RAC[0]);
      //return 0;
      exit(1);
    }
  }
  else{
    DEL=B*B-4.*A*C;
    if(DEL>=0.){
      DELSQ=sqrt(DEL);
      RAC[0]=(-B+DELSQ)/(2.*A);
      RAC[1]=(-B-DELSQ)/(2.*A);
      for(K=0;K<2;K++){
        if(RAC[K] >= RMIN && RAC[K] <= RMAX){
          rENERGY=E[J]-H*RAC[K];
          return rENERGY;
        }
      }
      printf("[RE_TAB:] PROBLEME 2 DANS ENERGY(Y) %f %f\n",RAC[0],RAC[1]);
      exit(1);
    }
    else{
      printf("[RE_TAB:] PROBLEME 3 DANS ENERGY(Y) %f\n",DEL);
      exit(1);
    }
  }
}
//C------------------
//C         EQUATION DU 3EME DEGRE X**3+A2*X**2+A1*X+A0=0
//C------------------
A2=3.*B2/A;
A1=(B4-B3-3*B2)/A;
A0=(B2+B3-YY)/A;
R0=(A1*A2-3.*A0)/6.-A2*A2*A2/27.;
Q=A1/3-A2*A2/9.;
D=Q*Q*Q+R0*R0;
if(D > 0.){
//C---------------
//C           1 RACINE REELLE
//C--------------
  DSQ=sqrt(D);
  RR=R0+DSQ;
  RHO=pow(fabs(RR),double(1.)/double(3.));
  if(RR < 0.){
    S1R=-RHO;
  }
  else{
    S1R=RHO;
  }
  RR=R0-DSQ;
  RHO=pow(fabs(RR),double(1.)/double(3.));
  if(RR < 0.){
    S2R=-RHO;
  }
  else{
    S2R=RHO;
  }
  RAC[0]=(S1R+S2R)-A2/3.;
  //if(RAC[0] >= RMIN && RAC[0] <= RMAX){
    rENERGY=E[J]-H*RAC[0];
    return rENERGY;
  //}
  //else{
  //  printf("PROBLEME 4 DANS ENERGY(Y) %f\n",RAC[0]);
  //  exit(1);
  //}
}
else{
//C-----------------
//C           3 RACINES REELLES
//C-----------------
  DSQ=sqrt(fabs(D));
  PHY=atan2(DSQ,R0)*RAD2DEG;
  if(PHY < 0.)PHY=PHY+360.;
  PHY=PHY/3.;
  CPHY=cos(PHY*DEG2RAD);
  SPHY=sin(PHY*DEG2RAD);
  RHO=pow((R0*R0+fabs(D)),double(1.)/double(6.));
  S1R=RHO*CPHY;
  S1I=RHO*SPHY;
  S2R=S1R;
  RAC[0]=2*S1R-A2/3.;
  RAC[1]=-S1R-A2/3.-sqrt(3.)*S1I;
  RAC[2]=-S1R-A2/3.+sqrt(3.)*S1I;
  for(K=0;K<3;K++){
    if(RAC[K] >= RMIN && RAC[K] <= RMAX){
      rENERGY=E[J]-H*RAC[K];
      return rENERGY;
    }
  }
  //printf("PROBLEME 5 DANS ENERGY(Y) %f %f %f\n",RAC[0],RAC[1],RAC[2]);
  //printf("PROBLEME 5 DANS ENERGY(Y) AI %f %f %f\n",A2,A1,A0);
  return E[J]-H*RAC[0];
  //exit(1);
}
}

//______________________________________________________________________________

int RE_TAB::A_EPAX(int Zp, int ZT, int AT){
// [2013-11-26 SK:] Parameter Zp was previously named Z, but to avoid
// shadowing a class member [-Wshadow] was changed to Zp

//--------------------------------------------------
//     CF K.SUMMERER PHYS. REV. 42(1990)2546
//--------------------------------------------------

  const int ZMIN=101;
  const int AMIN=301;
  const int AT0 = 0;

  int IAMIN,AA[ZMIN],ZZ[AMIN],IA,IZ;
  float Ap,Z_BETA,DEL,F_DEL;
  int a_epax;

  for(int i=0;i<ZMIN;i++) AA[i]=0;
  for(int i=0;i<AMIN;i++) ZZ[i]=0;

  if(AT != AT0){
    for(IA=5;IA<=AT-1;IA++){
      Ap=IA;
      F_DEL=1.;
      if(Ap/AT > 0.86) F_DEL=-51.*(Ap/AT-0.86)*(Ap/AT-0.86)+1;
      if(Ap < 66) DEL=2.041E-4*Ap*Ap;
      else       DEL=2.703E-2*Ap-0.895;
      Z_BETA=Ap/(1.98+0.0155*pow(double(Ap),double(0.6667)));
      ZZ[IA]=int(Z_BETA+DEL*F_DEL+0.5);
    }
    ZZ[AT]=ZT;
    IAMIN=5;
    for(IZ=3;IZ<=ZT-1;IZ++){
      for(IA=IAMIN;IA<=AT-1;IA++){
        if(IZ>=ZZ[IA] && IZ<=ZZ[IA+1]){
      AA[IZ]=int((IZ-ZZ[IA])/float(ZZ[IA+1]-ZZ[IA])+IA+0.5);
          IAMIN=IA;
          break;
        }
      }
    }
    AA[ZT]=AT;
  }
  a_epax=AA[Zp];
  if(a_epax == 0) a_epax = int(2.*Zp);
  if(Zp > ZT) a_epax = int((Zp-100.+sqrt(10000.+40.*Zp+Zp*Zp))/0.6); // GREEN
  if(Zp==1) a_epax = 1; 
  if(Zp==4) a_epax = 9;

  return a_epax;
}
//______________________________________________________________________________

//      IF(A1.EQ.0.)  A1 = (ZZ1-100.+SQRT(10000.+40.*ZZ1+ZZ1**2))/0.6
//______________________________________________________________________________