#include <math.h>
#include <algorithm>
#include <cassert>
#include "catima/calculations.h"
#include "catima/constants.h"
#include "catima/data_ionisation_potential.h"
#include "catima/data_atima.h"
#include "catima/generated_LS_coeff.h"
#include "catima/nucdata.h"
#include "catima/storage.h"
#ifdef GLOBAL
extern "C"
{
    #include "global/globallib.h"
}
#endif


namespace catima{

double reduced_energy_loss_unit(const Projectile &p, const Target &t){
    double zpowers = pow(p.Z,0.23)+pow(t.Z,0.23);
    double asum = p.A + t.A;
    return 32.53*t.A*1000*p.T*p.A/(p.Z*t.Z*asum*zpowers); //projectile energy is converted from MeV/u to keV
}


double dedx_n(const Projectile &p, const Material &mat){
    double w;
    double sum=0.0;
    for(int i=0;i<mat.ncomponents();i++){
        auto t = mat.get_element(i);
        w = mat.weight_fraction(i);
        sum += w*dedx_n(p,t);
    }
    return sum;   
}

double dedx_n(const Projectile &p, const Target &t){
    
    double zpowers = pow(p.Z,0.23)+pow(t.Z,0.23);
    double asum = p.A + t.A;
    double epsilon = 32.53*t.A*1000*p.T*p.A/(p.Z*t.Z*asum*zpowers); //projectile energy is converted from MeV/u to keV
    double sn=0;
    
    if(epsilon<=30){
        sn = log(1+(1.1383*epsilon))/ (2*(epsilon + 0.01321*pow(epsilon,0.21226) + 0.19593*pow(epsilon,0.5)));
    }
    else{
        sn = log(epsilon)/(2*epsilon);
    }
    sn = 100*8.4621*p.Z*t.Z*p.A*sn*Avogadro/(asum*zpowers*t.A);
    return sn;
}

double bethek_dedx_e(Projectile &p,const Material &mat, const Config &c){
    double w;
    double sum=0.0;
    for(int i=0;i<mat.ncomponents();i++){
        auto t = mat.get_element(i);
        w = mat.weight_fraction(i);
        sum += w*bethek_dedx_e(p,t,c,mat.I());
    }
    return sum;   
}

double bethek_dedx_e(Projectile &p, const Target &t, const Config &c, double I){
    assert(t.Z>0 && p.Z>0);
    assert(t.A>0 && p.A>0);
    assert(p.T>0.0);
    if(p.T==0)return 0.0;
    double gamma=1.0 + p.T/atomic_mass_unit;
    double beta2=1.0-1.0/(gamma*gamma);
    double beta = sqrt(beta2);
    double zp_eff = z_effective(p,t,c);    
    assert(zp_eff>=0);
    double Ipot = (I>0.0)?I:ipot(t.Z);
    assert(Ipot>0);
    double f1 = dedx_constant*pow(zp_eff,2.0)*t.Z/(beta2*t.A);
    double f2 = log(2.0*electron_mass*1000000*beta2/Ipot);

    double eta = beta*gamma;
    if(!(c.dedx&corrections::no_shell_correction) &&  eta>=0.13){ //shell corrections
        double cor = (+0.422377*pow(eta,-2)
                    +0.0304043*pow(eta,-4)
                    -0.00038106*pow(eta,-6))*1e-6*pow(Ipot,2) 
                  +(+3.858019*pow(eta,-2) 
                    -0.1667989*(pow(eta,-4))
                    +0.00157955*(pow(eta,-6)))*1.0e-9*pow(Ipot,3); 
        f2 = f2 -cor/t.Z;
    }
    f2+=2*log(gamma) -beta2;
    
    double barkas=1.0;
    if(!(c.dedx&corrections::no_barkas)){
        barkas = bethek_barkas(zp_eff,eta,t.Z);
        }
    
    double delta = bethek_density_effect(beta, t.Z);
    
    double LS = 0.0;
    if(!(c.dedx&corrections::no_lindhard)){
        //double LS = bethek_lindhard(p);
        LS = precalculated_lindhard(p);
        }
    double result  = (f2)*barkas + LS - delta/2.;
    result *=f1;
    
    return result;
}


double bethek_barkas(double zp_eff,double eta, double zt){
    double V2FVA[4]={0.33,0.30,0.26,0.23};
    double    VA[4]={1.,2.,3.,4.};
    double v1 = eta/(fine_structure*sqrt(zt));
    double v2fv;
    if(v1 >= 4){
        v2fv = 0.45 / sqrt(v1);
        }
    else if((v1 >= 1) && v1 < 4){//VALUES FROM THE JACKSON MC CARTHY FUNCTION //PHYS. REV. B 6 4131 P4136
        int i;
        for(i=1; i<4; i++){
	          if( VA[i] >= v1) break;
          }
	     v2fv = V2FVA[i-1]+(v1-VA[i-1])*(V2FVA[i]-V2FVA[i-1])/(VA[i]-VA[i-1]);
    }
    else{
        v2fv=0;
    }
    return 1.0+2.0 * zp_eff * v2fv /(v1*v1*sqrt(zt));
}

double bethek_density_effect(double beta, int zt){
    double gamma = 1/sqrt(1-(beta*beta));
    double x = log(beta * gamma) / 2.3025851;
    int i = zt-1;
    double del = 0;
    
    if(zt>97){  // check if data are available, if not take highest z data
        i=96;
        }
    if (x < density_effect::x0[i] ){
        if(density_effect::del_0[i] > 0.)del = density_effect::del_0[i] * pow(10.0,(2.*(x-density_effect::x0[i])));
        }
    else {
        del = 4.6052 * x - density_effect::c[i];
        if ( density_effect::x0[i]<= x &&  x <= density_effect::x1[i] ) del +=  density_effect::a[i] * pow((density_effect::x1[i] - x),density_effect::m[i]);
        }
    return del;
}


double bethek_lindhard(const Projectile &p){
    const double compton=3.05573356675e-3; // 1.18 fm / Compton wavelength
    double rho = exp(log(p.A)/3.0)*compton;
    double gamma=1.0 + p.T/atomic_mass_unit;
    double beta2=1.0-1.0/(gamma*gamma);
    double beta = sqrt(beta2);
    double eta = p.Z*fine_structure/beta;
    double beta_gamma_R = beta*gamma*rho;
    double sum = 0;
    int n=1;
    
    if(gamma < 10.0/rho){
        double dk[3];
        double dmk = 0;
        double dkm1 = 0;
        while(n<5000){
            double k0 = n;
            int max = (n==1)?3:2;
            for(int i=0;i<max;i++){
                double k;
                if(i==0)k=k0;
                if(i==1)k=-k0 - 1.0;
                if(i==2)k=-k0;
                double l = (k>0)?k:-k-1.0;
                double signk = (k>0)?1:((k<0)?-1:0);
                double sk = sqrt(k*k-fine_structure*fine_structure*p.Z*p.Z);
                std::complex<double> cexir_n (k,-eta/gamma);
                std::complex<double> cexir_den (sk,-eta);
                std::complex<double> cexir = std::sqrt(cexir_n/cexir_den);
                std::complex<double> csketa (sk + 1.0, eta);
                std::complex<double> cpiske(0.0,(M_PI*(l-sk)/2.0) - lngamma(csketa).imag());
                std::complex<double> cedr = cexir*std::exp(cpiske);
                double H=0;
                
                //std::complex<double> ceds(0.0,0.0);
                // finite struct part
                std::complex<double> cmsketa (-sk + 1.0, eta);
                std::complex<double> cexis_den (-sk,-eta);
                std::complex<double> cexis = std::sqrt(cexir_n/cexis_den);
                std::complex<double> cpimske(0.0,(M_PI*(l+sk)/2.0) - lngamma(cmsketa).imag());
                std::complex<double> ceds = cexis*std::exp(cpimske);
                std::complex<double> cmbeta_gamma_R(0,-beta_gamma_R);
                std::complex<double> c2beta_gamma_R(0,2.0*beta_gamma_R);
                std::complex<double> c2sk_1 (2.0*sk+1,0);
                std::complex<double> cm2sk_1 (-2.0*sk+1,0);
                std::complex<double> clambda_r = cexir*std::exp(cmbeta_gamma_R)*hyperg(csketa,c2sk_1,c2beta_gamma_R);
                std::complex<double> clambda_s = cexis*std::exp(cmbeta_gamma_R)*hyperg(cmsketa,cm2sk_1,c2beta_gamma_R);
                std::complex<double> cGrGs = lngamma(cm2sk_1);
                double GrGs = clambda_r.imag()/clambda_s.imag();
                GrGs *= exp( lngamma(csketa).real() 
                             - lngamma(cmsketa).real() 
                             - lngamma(c2sk_1).real()
                             + cGrGs.real() 
                             + 2.0*sk*log(2.0*beta_gamma_R));
                if(cos(cGrGs.imag()) < 1.0)GrGs*=-1;
                if(fabs(GrGs)>1.0e-9){
                    double FrGr = sqrt((gamma-1)/(gamma+1)) * clambda_r.real()/clambda_r.imag();
                    double FsGs = sqrt((gamma-1)/(gamma+1)) * clambda_s.real()/clambda_s.imag();
                    double gz = -1.0*signk*(rho*gamma + 1.5*p.Z*fine_structure);
                    double z1 = -1.0*signk*p.Z;
                    double b0 = 1.0;
                    double a0 = (1.0 + 2.0*fabs(k))*b0/(rho-gz);
                    double a1 = 0.5*(gz+rho)*b0;
                    double an = a1;
                    double anm1 = a0;
                    double bnm1 = b0;
                    double asum = a0;
                    double bsum = b0;
                    double nn = 1.0;
                    while(fabs(anm1/asum)>1e-6 && fabs(anm1/asum)>1e-6){
                        double bn = ((rho-gz)*an + fine_structure*z1*anm1/2.0)/(2.0*nn+2.0*fabs(k)+1.0);
                        double anp1 = ((gz+rho)*bn - fine_structure*z1*bnm1/2.0)/(2.0*nn + 2.0);
                        asum += an;
                        bsum += bn;
                        nn += 1.0;
                        anm1 = an;
                        an = anp1;
                        bnm1 = bn;
                    }
                    double figi= (k>0) ? asum/bsum : bsum/asum;
                    H = (FrGr - figi)/(figi-FsGs)* GrGs;
                    
                }
                else 
                    H = 0;
                
                
                dk[i] = std::arg(cedr + H*ceds);
            }
            if(n>1)dk[2] = dmk;
        
            double sdm2 = sin(dk[2]-dk[1]);
            double term1 = k0*(k0+1.0)*sdm2*sdm2/(eta*eta*(2.0*k0 + 1.0));
            if(n>1){
                double sd2 = sin(dk[0]-dkm1);
                term1 += k0*(k0-1.0)*sd2*sd2/(eta*eta*(2.0*k0 - 1.0));
            }
            double sdd = sin(dk[0]-dk[2]);
            double term2 = k0*sdd*sdd/(eta*eta*(4.0*k0*k0 - 1.0));
            double term3 = term1 - 1.0/k0;
            sum += term2 + term3;
            n += 1;
            dmk = dk[1];
            dkm1 = dk[0];
            
        }// end of while n<100
        
    }
    else{ // ultrarelativistic limit
        sum = -log(beta_gamma_R) - 0.2;
    }   
    return sum + (0.5*beta2);
}

double bethek_lindhard_X(const Projectile &p){
    const double compton=3.05573356675e-3; // 1.18 fm / Compton wavelength
    double rho = exp(log(p.A)/3.0)*compton;
    double gamma=1.0 + p.T/atomic_mass_unit;
    double beta2=1.0-1.0/(gamma*gamma);
    double beta = sqrt(beta2);
    double eta = p.Z*fine_structure/beta;
    double beta_gamma_R = beta*gamma*rho;
    double sum = 0;
    int n=1;
    
    if(1){
        double dk[4];
        double dmk = 0;
        double dmkp1 = 0;
        double dkm1 = 0;
        double dkm2 = 0;
        
        while(n<1000){
            double k0 = n;
            //int max = (n==1)?4:2;
            int max = 4;
            for(int i=0;i<max;i++){
                double k;
                if(i==0)k=k0;
                if(i==1)k=-k0 - 1.0;
                if(i==2 && n==1)k=-k0;
                if(i==3)k=-k0 - 2.0;
                double l = (k>0)?k:-k-1.0;
                double signk = (k>0)?1:((k<0)?-1:0);
                double sk = sqrt(k*k-fine_structure*fine_structure*p.Z*p.Z);
                std::complex<double> cexir_n (k,-eta/gamma);
                std::complex<double> cexir_den (sk,-eta);
                std::complex<double> cexir = std::sqrt(cexir_n/cexir_den);
                std::complex<double> csketa (sk + 1.0, eta);
                std::complex<double> cpiske(0.0,(M_PI*(l-sk)/2.0) - lngamma(csketa).imag());
                std::complex<double> cedr = cexir*std::exp(cpiske);
                double H=0;
                
                //std::complex<double> ceds(0.0,0.0);
                // finite struct part
                std::complex<double> cmsketa (-sk + 1.0, eta);
                std::complex<double> cexis_den (-sk,-eta);
                std::complex<double> cexis = std::sqrt(cexir_n/cexis_den);
                std::complex<double> cpimske(0.0,(M_PI*(l+sk)/2.0) - lngamma(cmsketa).imag());
                std::complex<double> ceds = cexis*std::exp(cpimske);
                std::complex<double> cmbeta_gamma_R(0,-beta_gamma_R);
                std::complex<double> c2beta_gamma_R(0,2.0*beta_gamma_R);
                std::complex<double> c2sk_1 (2.0*sk+1,0);
                std::complex<double> cm2sk_1 (-2.0*sk+1,0);
                std::complex<double> clambda_r = cexir*std::exp(cmbeta_gamma_R)*hyperg(csketa,c2sk_1,c2beta_gamma_R);
                std::complex<double> clambda_s = cexis*std::exp(cmbeta_gamma_R)*hyperg(cmsketa,cm2sk_1,c2beta_gamma_R);
                std::complex<double> cGrGs = lngamma(cm2sk_1);
                double GrGs = clambda_r.imag()/clambda_s.imag();
                GrGs *= exp( lngamma(csketa).real() 
                             - lngamma(cmsketa).real() 
                             - lngamma(c2sk_1).real()
                             + cGrGs.real() 
                             + 2.0*sk*log(2.0*beta_gamma_R));
                if(cos(cGrGs.imag()) < 1.0)GrGs*=-1;
                if(fabs(GrGs)>1.0e-9){
                    double FrGr = sqrt((gamma-1)/(gamma+1)) * clambda_r.real()/clambda_r.imag();
                    double FsGs = sqrt((gamma-1)/(gamma+1)) * clambda_s.real()/clambda_s.imag();
                    double gz = -1.0*signk*(rho*gamma + 1.5*p.Z*fine_structure);
                    double z1 = -1.0*signk*p.Z;
                    double b0 = 1.0;
                    double a0 = (1.0 + 2.0*fabs(k))*b0/(rho-gz);
                    double a1 = 0.5*(gz+rho)*b0;
                    double an = a1;
                    double anm1 = a0;
                    double bnm1 = b0;
                    double asum = a0;
                    double bsum = b0;
                    double nn = 1.0;
                    while(fabs(anm1/asum)>1e-6 && fabs(anm1/asum)>1e-6){
                        double bn = ((rho-gz)*an + fine_structure*z1*anm1/2.0)/(2.0*nn+2.0*fabs(k)+1.0);
                        double anp1 = ((gz+rho)*bn - fine_structure*z1*bnm1/2.0)/(2.0*nn + 2.0);
                        asum += an;
                        bsum += bn;
                        nn += 1.0;
                        anm1 = an;
                        an = anp1;
                        bnm1 = bn;
                    }
                    double figi= (k>0) ? asum/bsum : bsum/asum;
                    H = (FrGr - figi)/(figi-FsGs)* GrGs;
                    
                }
                else 
                    H = 0;
                
                
                dk[i] = std::arg(cedr + H*ceds);
            }
            if(n>1)dk[2] = dmk;
        
            double strterm1p = 0;
            double strterm1n = 0;
            double strterm2 = 0;
            double strterm3 = 0;
            double eta2 = eta*eta;
            
            double sdm2 = sin(dk[0]-dkm2);
            if(n>2){
                strterm1p = sdm2*sdm2*(k0-1)*(k0-2)/((2.0*k0 - 1.0)*(2.0*k0-3.0));
            }
            
            sdm2 = sin(dk[2]-dk[3]);
            strterm1n = sdm2*sdm2*(-k0-1)*(-k0-2)/((-2.0*k0 - 1.0)*(-2.0*k0-3.0));
            
            if(n>1){
                double sd2 = sin(dk[0]-dmkp1);
                strterm2 += (k0-1.0)*sd2*sd2/((2.0*k0 - 3.0)*(4.0*k0*k0 - 1.0));
            }
            
            double sdd = sin(dk[0]-dk[1]);
            strterm3 = sdd*sdd*(k0+1.0)*((1/(4.0*k0*k0 -1.0))+(1/(4*(k0+1.0)*(k0+1.0) - 1.0)))/(2.0*k0 + 1.0);
            
            //sum += k0*(strterm1p + strterm1n + (strterm2*2) + strterm3)/eta2;
            sum += k0*(strterm1p + strterm1n + (strterm2*2) + strterm3)/eta2;
            sum += - (2.0/k0);
            //std::cout<<n<<" "<<strterm1p<<" "<<strterm1n<<" "<<strterm2<<" "<<strterm3<<" "<<sum<<std::endl;
            n += 1;
            dmk = dk[1];
            dkm2 = dkm1;
            dkm1 = dk[0];
            dmkp1 = dk[2];
            
        }// end of while n<100
        
    }
    else{ // ultrarelativistic limit
        
    }
    return 2*bethek_lindhard(p) - sum - beta2;
    //return sum;
}



double sezi_p_se(double energy,const Target &t){
    double sp = -1;
    double e = 1000*energy; //e in keV/u
    int i = t.Z - 1;
    
    if(e<=25)e=25;
    //double sl = (proton_stopping_coef[i][0]*pow(e,proton_stopping_coef[i][1])) + (proton_stopping_coef[i][2]*pow(e,proton_stopping_coef[i][3]));
    //double sh = proton_stopping_coef[i][4]/pow(e,proton_stopping_coef[i][5]) * log( (proton_stopping_coef[i][6]/e) + (proton_stopping_coef[i][7]*e));
    double sl = (proton_stopping_coef[i][0]*catima::power(e,proton_stopping_coef[i][1])) + (proton_stopping_coef[i][2]*catima::power(e,proton_stopping_coef[i][3]));
    double sh = proton_stopping_coef[i][4]/catima::power(e,proton_stopping_coef[i][5]) * log( (proton_stopping_coef[i][6]/e) + (proton_stopping_coef[i][7]*e));
    sp = sl*sh/(sl+sh);
    e=1000*energy;
    if(e<=25){
        //sp *=(t.Z>6)?pow(e/25,0.45):pow(e/25,0.25);
        sp *=(t.Z>6)?catima::power(e/25,0.45):catima::power(e/25,0.25);
    }
    
    return 100*sp*Avogadro/t.A;
}

double sezi_dedx_e(const Projectile &p, const Target &t){
    double e=p.T*1000; // e in keV/u
    double se = 0;
    
    if(p.Z==1){
        return sezi_p_se(p.T,t);
    }
    else if(p.Z == 2){
        double a=0;
        double b=0;
        //double zeta = 0;
        
        if(e<=1)e=1;
        // He Zeff
        b = log(e);
        a = 0.2865 + b*(0.1266+ b*(-0.001429+ b*(0.02402 + b*(-0.01135 + b*0.001475))));
        double heh = 1.0 - exp(-std::min(30.,a));
        b = 7.6 - std::max(0., b);
        a = (1.0 + (0.007 + 0.00005*t.Z)*exp(- b*b ));
        heh *= a*a;
        //zeta = sqrt(heh);
        se = sezi_p_se(p.T,t)*heh*4.0; //scale proton stopping
        if(e==1)se*= sqrt(p.T*1000.0);  //vel proportional
        return se;
    }
    else{ // heavy ion
        double h1,h2,h3,h4;
        double a,q,b;
        double l1,l0,l;
        double YRmin = 0.130; //  YRmin = VR / ZP**0.67 <= 0.13 OR VR <= 1.0
        double VRmin = 1.0;
        double v=0;
        double vfermi = atima_vfermi[(int)t.Z-1];
        double yr=0;
        double zeta = 0;
        double se;
        
        v = sqrt(e/25.0)/vfermi;
        double v2=v*v;
        
        double vr = (v >= 1)? v*vfermi*(1.+ 1./(5.*v2)) : 3.0*vfermi/4.0*(1.0+v2*(2.0/3.0-v2/15.0));
        
        h1= 1./catima::power(p.Z,0.6667);
        yr = std::max(YRmin,vr*h1);
        yr = std::max(yr, VRmin*h1);
        
        //--  CALCULATE ZEFF
        a = -0.803*catima::power(yr,0.3) + 1.3167*catima::power(yr,0.6) + 0.38157*yr + 0.008983*yr*yr;
        q = std::min(1.0, std::max(0.0 , (1.0 - exp(-std::min(a, 50.0))))); //-- Q = IONIZATION LEVEL OF THE ION AT RELATIVE VELOCITY YR
        
        //-- IONIZATION LEVEL TO EFFECTIVE CHARGE
        h1 = 1./ catima::power(p.Z,0.3333);
        
        b = (std::min(0.43, std::max(0.32,0.12 + 0.025*p.Z)))*h1;
        l0 = (.8 - q * std::min(1.2,0.6 +p.Z/30.0))*h1;
        if(q < 0.2){
            l1 = 0;
            }
        else{
            if (q < std::max(0.0,0.9-0.025*p.Z)){
                l1 = b*(q-0.2)/fabs(std::max(0.0,0.9-0.025*p.Z)-0.2000001);
                }
	     	 else{
             	if(q < std::max(0.0,1.0 - 0.025*std::min(16.,p.Z))) l1 = b;
                else l1 = b*(1.0 - q)/(0.025*std::min(16.,p.Z));
                }
        }
        // calculate screening
        l = std::max(l1,l0*atima_lambda_screening[(int)p.Z-1]);
        h1 =4.0*l*vfermi/1.919;
        zeta = q + (1./(2.*(vfermi*vfermi)))*(1. - q)* log(1. + h1*h1);
         
         // ZP**3 EFFECT AS IN REF. 779?
        a = 7.6 - std::max(0.0, log(e));
        zeta = zeta*(1. + (1./(p.Z*p.Z))*(0.18 + .0015*t.Z)*exp(-a*a));
        
        h1= 1./catima::power(p.Z,0.6667);
        if (yr <= ( std::max(YRmin, VRmin*h1))){
            VRmin=std::max(VRmin, YRmin/h1);
            //--C        ..CALCULATE VELOCITY STOPPING FOR YR < YRmin
            double vmin =.5*(VRmin + sqrt(std::max(0.0,VRmin*VRmin - .8*vfermi*vfermi)));
            double eee = 25.0*vmin*vmin;
            double eval = 1;
            if((t.Z == 6) || (((t.Z == 14) || (t.Z == 32)) &&  (p.Z <= 19))) eval = 0.35;
            else eval = 0.5;
            
            h1 = zeta *p.Z;
            h4 = catima::power(e / eee,eval);
            se = sezi_p_se(eee*0.001,t) * h1*h1*h4;
            return se;
        }
        else {
            se = sezi_p_se(p.T,t)*catima::power(zeta*p.Z,2.0);
            return se;
        }
        return 0;
    }
};

double sezi_dedx_e(const Projectile &p, const Material &mat){
    double w;
    double sum=0.0;
    for(int i=0;i<mat.ncomponents();i++){
        auto t = mat.get_element(i);
        w = mat.weight_fraction(i);
        sum += w*sezi_dedx_e(p,t);
    }
    return sum;
}


double gamma_from_T(double T){
    return 1.0 + T/atomic_mass_unit;
};

double beta_from_T(double T){
    double gamma = gamma_from_T(T);
    return sqrt(1.0-1.0/(gamma*gamma));
}

double energy_straggling_firsov(double z1,double energy, double z2, double m2){
    double gamma = gamma_from_T(energy);
    double beta2=1.0-1.0/(gamma*gamma);
    double factor=4.8184E-3*pow(z1+z2,8.0/3.0)/m2;
    return factor*beta2/fine_structure/fine_structure;
    }

double angular_scattering_variance(Projectile &p, Target &t){
    if(p.T<=0)return 0.0;
    double e=p.T;
    double _p = sqrt(e*(e+2*atomic_mass_unit))*p.A;
    double beta = _p/((e+atomic_mass_unit)*p.A);
    double lr = radiation_length(t.Z,t.A);
    return 198.81 * pow(p.Z,2)/(lr*pow(_p*beta,2));
}

/// radioation lengths are taken frm Particle Data Group 2014
double radiation_length(int z, double m){
    double lr = 0;
    if(z==1){return 63.04;}
    if(z==2){return 94.32;}
    if(z==3){return 82.78;}
    if(z==4){return 65.19;}
    if(z==6){return 42.7;}
    if(z==7){return 37.99;}
    if(z==8){return 34.24;}
    if(z==9){return 32.93;}
    if(z==10){return 28.94;}
    if(z==13){return 24.01;}
    if(z==14){return 21.82;}
    if(z==17){return 19.28;}
    if(z==18){return 19.55;}
    if(z==22){return 16.16;}
    if(z==26){return 13.84;}
    if(z==29){return 12.86;}
    if(z==32){return 12.25;}
    if(z==50){return 8.82;}
    if(z==54){return 8.48;} 
    if(z==74){return 6.76;}
    if(z==78){return 6.54;}
    if(z==79){return 6.46;}
    if(z==82){return 6.37;}
    if(z==92){return 6.00;}
    
    double z2 = z*z;
    double z_13 = 1.0/pow(z,1./3.);
    double z_23 = z_13*z_13;
    double a2 = fine_structure*fine_structure*z2;
    double a4 = a2*a2;
    double a6 = a4*a2;
    lr= 716.405*m/(z2* (log(184.15*z_13) + log(1194.0*z_23)/z - -1.202*a2 + 1.0369*a4 - 1.008*a6/(1+a2) ) );
    return lr;
}


double precalculated_lindhard(const Projectile &p){
    double T = p.T;
    int z = (int)p.Z ;
    if(z>LS_MAX_Z)z=LS_MAX_Z;
    if(p.T<ls_coefficients::ls_energy_table(0))T=ls_coefficients::ls_energy_table(0);
    
    double da = (p.A - element_atomic_weight(z))/element_atomic_weight(z);
    z = z-1;
    //catima::Interpolator ls_a(ls_coefficients::ls_energy_points,ls_coefficients::ls_coefficients_a[z],LS_NUM_ENERGY_POINTS,interpolation_t::linear);
    //catima::Interpolator ls_ahi(ls_coefficients::ls_energy_points,ls_coefficients::ls_coefficients_ahi[z],LS_NUM_ENERGY_POINTS,interpolation_t::linear);
    //catima::Interpolator ls_a(ls_coefficients::ls_energy_table.values,ls_coefficients::ls_coefficients_a[z],LS_NUM_ENERGY_POINTS,interpolation_t::cspline);
    //catima::Interpolator ls_ahi(ls_coefficients::ls_energy_table.values,ls_coefficients::ls_coefficients_ahi[z],LS_NUM_ENERGY_POINTS,interpolation_t::cspline);
    double v1 = EnergyTable_interpolate(ls_coefficients::ls_energy_table,T,ls_coefficients::ls_coefficients_a[z]);
    double v2 = EnergyTable_interpolate(ls_coefficients::ls_energy_table,T,ls_coefficients::ls_coefficients_ahi[z]);

    //double dif = ls_ahi(T) - ls_a(T);
    //return ls_a(T)+(dif*da/ls_coefficients::a_rel_increase);
    double dif = v2 - v1;
    return v1+(dif*da/ls_coefficients::a_rel_increase);
    
}

double precalculated_lindhard_X(const Projectile &p){
    double T = p.T;
    int z = (int)p.Z ;
    if(z>LS_MAX_Z)z=LS_MAX_Z;
    //if(p.T<ls_coefficients::ls_energy_table(0))T=ls_coefficients::ls_energy_table(0);
    if(p.T<ls_coefficients::ls_energy_table(0))
		return 1.0;
    double da = (p.A - element_atomic_weight(z))/element_atomic_weight(z);
    z = z-1;
    
    //catima::Interpolator ls_X_a(ls_coefficients::ls_energy_table.values,ls_coefficients::ls_X_coefficients_a[z],LS_NUM_ENERGY_POINTS,interpolation_t::linear);
    //catima::Interpolator ls_X_ahi(ls_coefficients::ls_energy_table.values,ls_coefficients::ls_X_coefficients_ahi[z],LS_NUM_ENERGY_POINTS,interpolation_t::linear);
    double v1 = EnergyTable_interpolate(ls_coefficients::ls_energy_table,T,ls_coefficients::ls_X_coefficients_a[z]);
    double v2 = EnergyTable_interpolate(ls_coefficients::ls_energy_table,T,ls_coefficients::ls_X_coefficients_ahi[z]);
    
    //double dif = ls_X_ahi(T) - ls_X_a(T);
    //return ls_X_a(T)+(dif*da/ls_coefficients::a_rel_increase);
    double dif = v2 - v1;
    return v1+(dif*da/ls_coefficients::a_rel_increase);
}

double dedx_variance(Projectile &p, Target &t, const Config &c){
    double gamma = gamma_from_T(p.T);
    double cor=0;
    double beta = beta_from_T(p.T);
    double beta2 = beta*beta;
    //double zp_eff = z_effective(p,t,c);
    double zp_eff = z_eff_Pierce_Blann(p.Z, beta);
    double f = domega2dx_constant*pow(zp_eff,2)*t.Z/t.A;

    if(c.dedx_straggling == omega::atima){
        cor = 24.89 * pow(t.Z,1.2324)/(electron_mass*1e6 * beta2)*
			log( 2.0*electron_mass*1e6*beta2/(33.05*pow(t.Z,1.6364)));
	    cor = std::max(cor, 0.0 );
    }
	//double X = bethek_lindhard_X(p);
    double X = precalculated_lindhard_X(p);
    X *= gamma*gamma;
    if(p.T<30.0)
		return std::min(f*(X+cor), energy_straggling_firsov(p.Z, p.T, t.Z,t.A));
	else
		return f*(X+cor);
}

double z_effective(const Projectile &p,const Target &t, const Config &c){
    if(c.z_effective == z_eff_type::none){
        return p.Q;
    }

    double gamma=1.0 + p.T/atomic_mass_unit;
    double beta = sqrt(1.0-1.0/(gamma*gamma));
    if(c.z_effective == z_eff_type::pierce_blann){
        return z_eff_Pierce_Blann(p.Z, beta);
    }
    else if(c.z_effective == z_eff_type::anthony_landorf){
        return z_eff_Anthony_Landford(p.Z, beta, t.Z);
    }
    
    else if(c.z_effective == z_eff_type::hubert){
        return z_eff_Hubert(p.Z, p.T, t.Z);
    }
    else if(c.z_effective == z_eff_type::winger){
        return z_eff_Winger(p.Z, beta, t.Z);
    }
    else if(c.z_effective == z_eff_type::global){
        return z_eff_global(p.Z, p.T, t.Z);
    }
    else if(c.z_effective == z_eff_type::atima14){
        return z_eff_atima14(p.Z, p.T, t.Z);
    }
    else if(c.z_effective == z_eff_type::schiwietz){
        return z_eff_Schiwietz(p.Z, beta, t.Z);
    }
    else{
        return 0.0;
    }
}

double z_eff_Pierce_Blann(double z, double beta){
    return z*(1.0-exp(-0.95*fine_structure_inverted*beta/pow(z,2.0/3.0))); 
}

double z_eff_Anthony_Landford(double pz, double beta, double tz){
    double B = 1.18-tz*(7.5e-03 - 4.53e-05*tz);
    double A = 1.16-tz*(1.91e-03 - 1.26e-05*tz);
    return pz*(1.0-(A*exp(-fine_structure_inverted*B*beta/pow(pz,2.0/3.0)))); 
}

double z_eff_Hubert(double pz, double E, double tz){
    double lntz = log(tz);
    double x1,x2,x3,x4;

    if(E<2.5)
        return 0.0;
    if(tz == 4){
        x1 = 2.045 + 2.0*exp(-0.04369*pz);
        x2 = 7.0;
        x3 = 0.2643;
        x4 = 0.4171;
    }
    else if(tz==6){
        x1 = 2.584 + 1.91*exp(-0.03958*pz);
        x2 = 6.933;
        x3 = 0.2433;
        x4 = 0.3969;
    }
    else{
        x1 = (1.164 + 0.2319*exp(-0.004302*tz)) + 1.658*exp(-0.0517*pz);
        x2 = 8.144 + 0.09876*lntz;
        x3 = 0.314 + 0.01072*lntz;
        x4 = 0.5218 + 0.02521*lntz;
    }
     
    return pz*(1-x1*exp(-x2*catima::power(E,x3)*catima::power(pz,-x4)));
}

double z_eff_Winger(double pz, double beta, double tz){
	double c0,c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13;
	double x, lnz, lnzt, a0,a1,a2,a3,a4; 

	c0 = 0.4662;
	c1 = 0.5491;
	c2 = 0.7028;
	c3 = 0.1089;
	c4 = 0.001644;
	c5 = 0.5155;
	c6 = 0.05633;
	c7 = 0.005447;
	c8 = 0.8795;
	c9 = 1.091;
	c10= 0.0008261;
	c11= 2.848;
	c12= 0.2442;
	c13= 0.00009293;
	// the following is from Helmut to correct winger formula for H and He target
    if(tz==1){
        tz = 2.6;
    }
    if(tz==2){
        tz = 2.8;
    }
	x = beta /0.012 /pow(pz,0.45);
	lnz =log(pz);
	lnzt=log(tz);

	a0 = -c0;
	a1 = -c1  * exp( c2 *lnz - c3 *lnz*lnz +c4*lnz*lnz*lnz -c5*lnzt + c6 *lnzt*lnzt);
	a2 =  c7  * exp( c8 *lnz - c9 *lnzt);
	a3 = -c10 * exp( c11*lnz - c12*lnz*lnz*lnz);
	a4 = -c13;

    return pz * (1. - exp(a0 +a1*x +a2*x*x +a3*x*x*x +a4*x*x*x*x));	
	}

double z_eff_global(double pz, double E, double tz){
    if(E>2000)
        return pz;
    else
        #ifdef GLOBAL
        return global_qmean(pz, tz, E);
        #else
        return -1;
        #endif
}

double z_eff_Schiwietz(double pz, double beta, double tz){
    double scaled_velocity;
    double c1, c2;
    double x;

    scaled_velocity = catima::power(pz,-0.543)*beta/bohr_velocity;
    c1 = 1-0.26*exp(-tz/11.0)*exp(-(tz-pz)*(tz-pz)/9);
    c2 = 1+0.030*scaled_velocity*log(tz);
    x = c1*catima::power(scaled_velocity/c2/1.54,1+(1.83/pz));
    return pz*((8.29*x) + (x*x*x*x))/((0.06/x) + 4 + (7.4*x) + (x*x*x*x) );

}

double z_eff_atima14(double pz, double T, double tz){
    double qpb;
    double qhigh,qwinger,qglobal,qmean=0;
    double c1 = 1.4;
    double c2 = 0.28;
    double c3 = 0.04;
    double beta = beta_from_T(T);
    double gamma = gamma_from_T(T);
    double emax, emin;
    qpb = z_eff_Pierce_Blann(pz,beta);
    #ifdef GLOBAL

    if(T>30.0 && T<1500.0 && pz>28){
        qglobal = z_eff_global(pz,T,tz);
        qglobal = (qglobal - qpb)*c1/catima::power(tz+1,c2)*(1-exp(-c3*T)) + qpb;
    }
    
    emax = 1500.;
    emin = 1000.;
    if(T>emax){
        qhigh = pz * (1.0-exp(-180.0*beta*catima::power(gamma,0.18)*catima::power(pz,-0.82)*catima::power(tz,0.1)));
        qmean = qhigh;
    }
    else if(T>=emin && T<=emax){
        qhigh = pz * (1.0-exp(-180.0*beta*catima::power(gamma,0.18)*catima::power(pz,-0.82)*catima::power(tz,0.1)));
        if(pz<=28){
            qwinger = z_eff_Winger(pz,beta,tz);
            qmean = ((emax-T)*qwinger + (T-emin)*qhigh)/(emax-emin);
        }
        else{
            qmean = ((emax-T)*qglobal + (T-emin)*qhigh)/(emax-emin);
        }
    }
    else{
        emax = 70.0;
        emin = 30.0;
        if(pz<=28){
            qwinger = z_eff_Winger(pz,beta,tz);
            qmean = qwinger;
        }
        else{
            if(T>=emax){
                qmean = qglobal;
                }
            else if(T<emin){
                qwinger = z_eff_Winger(pz,beta,tz);
                qmean = qwinger;
                }
            else
                {
                qwinger = z_eff_Winger(pz,beta,tz);
                qmean = ((emax-T)*qwinger + (T-emin)*qglobal)/(emax-emin);
                }
        }
    }
    #endif
    return qmean;
    }

std::complex<double> hyperg(const std::complex<double> &a,
                                                const std::complex<double> &b,
                                                const std::complex<double> &z){
        double dm = 0.0;
        std::complex<double> term(1.0, 0.0);
        std::complex<double> sumterm(1.0, 0.0);
        std::complex<double> previousterm;
        do {
            previousterm = term;
            dm += 1.0;
            std::complex<double> Cm(dm-1.0, 0.0);
            term = previousterm * ((a + Cm)/(b + Cm)) * (z/dm);
            sumterm += term;
        } while( std::abs(term) > 1.0e-6 && std::abs(previousterm) > 1.0e-6 );
        return(sumterm);
}

std::complex<double> lngamma( const std::complex<double> &z )
    {
        const static double coeff[6] = {76.18009172947146,
            -86.50532032941677,
            24.01409824083091,
            -1.231739572450155,
            0.1208650973866179e-2,
            -0.5395239384953e-5};
        double x, y;
        if(z.real() > 0) {
            x=z.real()-1.0;
            y=z.imag();
        } else {
            x=-z.real();
            y=-z.imag();
        }
        double r = sqrt((x+5.5)*(x+5.5)+y*y);
        double aterm1=y*log(r);
        double aterm2=(x+0.5)*atan2(y,(x+5.5))-y;
        double lterm1=(x+0.5)*log(r);
        double lterm2=-y*atan2(y,(x+5.5)) - (x+5.5) + 0.5*log(2.0*M_PI);
        double num=0.0;
        double denom=1.000000000190015;
        for(int j=1;j<7;j++){
            double fj=(double)j;
            double cterm=coeff[j-1]/((x+fj)*(x+fj)+y*y);
            num+=cterm;
            denom+=(x+fj)*cterm;
        }
        num*=-y;
        double aterm3=atan2(num,denom);
        double lterm3 = 0.5*log(num*num + denom*denom);
        std::complex<double> result(lterm1+lterm2+lterm3,aterm1+aterm2+aterm3);
        if(z.real() < 0){
            std::complex<double> lpi(log(M_PI), 0.0);
            result = lpi - (result + std::log(std::sin(M_PI*z)));
        }
        return(result);
    }

}