{

if (l == 0 || l == 1) return 1.0;

else return l*factorial(l-1);

{

if (l == -1 || l == 0) return 1.0;

else return l * dualfac (l-2);

(

const double *rad,

const double *rad_f,

const int& mesh,

const double &yp1, // if yp1 > ypmax, consider the second derivative

const double &ypn,

double* y2

)

{

ModuleBase::timer::start("Mathzone_Add1","SplineD2");

double dx1 = 0.0, dx2 = 0.0, dy1 = 0.0, dy2 = 0.0, p = 0.0, qn = 0.0, sig = 0.0, un = 0.0;

double * u = nullptr;

u = new double[mesh-1];

const double ypmax = 99999.00;

if (yp1 > ypmax)

{

y2[0] = 0.0;

u[0] = 0.0;

}

else

{

y2[0] = -0.5;

dx1 = rad[1] - rad[0];

dy1 = rad_f[1] - rad_f[0];

u[0] = ( 3.0 / dx1 ) * (dy1 / dx1 - yp1);

}

for(int i = 1; i < mesh-1 ; i++)

{

dx1 = rad[i] - rad[i-1];

dx2 = rad[i+1] - rad[i-1];

dy1 = rad_f[i+1] - rad_f[i];

dy2 = rad_f[i] - rad_f[i-1];

sig = dx1 / dx2;

p = sig * y2[i-1] + 2.0;

y2[i] = (sig - 1.0) / p;

u[i] = dy1 / (dx2 - dx1) - dy2 / dx1;

u[i] = (6.0 * u[i] / dx2 - sig * u[i-1]) / p;

}

if (ypn > ypmax)

{

qn = un = 0.0;

}

else

{

dx1 = rad[mesh-1] - rad[mesh-2];

dy1 = rad_f[mesh-1] - rad_f[mesh-2];

qn = 0.5;

un = 3.0 / dx1 * (ypn - dy1 / dx1);

}

y2[mesh-1] = (un - qn * u[mesh-2]) / (qn * y2[mesh-2] + 1.0);

for(int i = mesh-2; i >= 0; i--)

{

y2[i] = y2[i] * y2[i+1] + u[i];

}

delete[] u;

ModuleBase::timer::end("Mathzone_Add1","SplineD2");

(

const double * const rad,

const double * const rad_f,

const double * const y2,

const int& mesh,

const double * const r,

const int& rsize,

double * const y,

double * const dy

)

{

ModuleBase::timer::start("Mathzone","cubic_spline");

#ifdef _OPENMP

#pragma omp parallel for schedule(static)

#endif

for(int m = 0; m < rsize ; m++)

{

int klo = 0;

int khi = mesh-1;

while (khi - klo > 1)

{

const int k = (khi + klo) / 2 ;

if(rad[k] > r[m]) khi = k;

else klo = k;

}

const double h = rad[khi] - rad[klo];

if(h == 0.0)

{

std::cout << "Cubic_Spline_Interpolation, h == 0.0 so that cannot be divided" << std::endl;

//ModuleBase::WARNING_QUIT("Cubic_Spline_Interpolation","h == 0.0 so that cannot be divided");

exit(0);

}

const double a = (rad[khi] - r[m]) / h;

const double b = (r[m] - rad[klo]) / h;

const double dy_tmp = (rad_f[khi] - rad_f[klo]) / h -

(3.0 * a * a - 1.0) / 6.0 * h * y2[klo] +

( 3.0 * b * b - 1.0) / 6.0 * h * y2[khi];

dy[m] = dy_tmp;

const double y_tmp = a * rad_f[klo] + b * rad_f[khi] + ((a*a*a - a) * y2[klo] + (b*b*b - b) * y2[khi]) * (h*h) / 6.0;

y[m] = y_tmp;

//const double ddy_tmp = a * y2[klo] + b * y2 [khi];

//ddy[m] = ddy_tmp;

}

ModuleBase::timer::end("Mathzone","cubic_spline");

(

const double* old_phi,

const int& msh,

const double& dr,

const double& newr

)

{

double h1 = 0.0, h2 = 0.0, h3 = 0.0, f1 = 0.0, f2 = 0.0, f3 = 0.0, f4 = 0.0;

double g1 = 0.0, g2 = 0.0, x1 = 0.0, x2 = 0.0, y1 = 0.0, y2 = 0.0, f = 0.0;

double result = 0.0;

double rmax = (msh-1) * dr;

if (newr < 0.0)

{

throw std::runtime_error("newr should >= 0. "+std::string(__FILE__)+" line "+std::to_string(__LINE__));

}

else if ( rmax <= newr )

{

f = 0.0;

}

else

{

// the old position.

const int m = static_cast<int> (newr / dr) + 1;

// Spline like interpolation

if (m == 1)

{

h2 = dr;

h3 = dr;

f2 = old_phi[m-1];

f3 = old_phi[m];

f4 = old_phi[m+1];

h1 = -(h2+h3);

f1 = f4;

}

else if (m == (msh-1))

{

h1 = dr;

h2 = dr;

f1 = old_phi[m-2];

f2 = old_phi[m-1];

f3 = old_phi[m];

h3 = -(h1+h2);

f4 = f1;

}

else

{

h1 = dr;

h2 = dr;

h3 = dr;

f1 = old_phi[m-2];

f2 = old_phi[m-1];

f3 = old_phi[m];

f4 = old_phi[m+1];

}

// Calculate the value at newr

g1 = ((f3-f2)*h1/h2 + (f2-f1)*h2/h1)/(h1+h2);

g2 = ((f4-f3)*h2/h3 + (f3-f2)*h3/h2)/(h2+h3);

x1 = newr - (m-1)*dr;

x2 = newr - m*dr;

y1 = x1/h2;

y2 = x2/h2;

f = y2*y2*(3.0*f2 + h2*g1 + (2.0*f2 + h2*g1)*y2)

+ y1*y1*(3.0*f3 - h2*g2 - (2.0*f3 - h2*g2)*y1);

}

result = f;

return result;

(

const double *radf,

const int &mesh,

const double &dr,

const int &nd,

double *phind

)

{

ModuleBase::timer::start("Mathzone_Add1", "Uni_Deriv_Phi");

int FFT_NR = 2*mesh-1; // FFT_NR = 16019

// std::cout << "\n mesh=" << mesh << ", radf[8010]=" << radf[8010] << ", radf[8009]=" << radf[8009] ;

// mesh=8010, radf[8010]=4.396478951532926e-01, radf[8009]=0.000000000000000e+00

fftw_complex *fft_phir = new fftw_complex[FFT_NR];

fftw_complex *fft_phik = new fftw_complex[FFT_NR];

fftw_complex *fft_ndphik = new fftw_complex[FFT_NR];

fftw_complex *fft_ndphir = new fftw_complex[FFT_NR];

fftw_plan p1;

fftw_plan p2;

////CAREFUL: POINT 0 is OF GOOD IMPORTANCE

//for (int ir = 0; ir < FFT_NR/2; ++ir)

//{

//fft_phir[ir].re = radf[ir];

//fft_phir[ir].im = 0.0;

//}

//for (int ir = FFT_NR/2; ir < FFT_NR; ++ir)

//{

//int jr = FFT_NR - ir;

//fft_phir[ir].re = radf[jr];

//fft_phir[ir].im = 0.0;

//}

//

// second call: different value at [8010]; FFT_NR = 16019, FFT_NR/2=8009 ::

// CAREFUL: POINT 0 is OF GOOD IMPORTANCE

for (int ir = 0; ir < FFT_NR/2; ++ir) // ik = 0 1 ... 8008

{

c_re(fft_phir[ir]) = radf[ir];

c_im(fft_phir[ir]) = 0.0;

}

for (int ir = FFT_NR/2; ir < FFT_NR; ++ir) // ir = 8009 8010 ... 16018

{

//int jr = FFT_NR - ir ; // jr = 8010 8009 ... 1

int jr = FFT_NR - ir -1 ; // -> 8009 8008 ... 0

c_re(fft_phir[ir]) = radf[jr];

c_im(fft_phir[ir]) = 0.0;

}

// FFTW

#if defined __FFTW3

//std::cout << "\n Call FFTW3 ";

p1 = fftw_plan_dft_1d(FFT_NR, fft_phir, fft_phik, FFTW_FORWARD, FFTW_ESTIMATE);

fftw_execute(p1);

//fftw_destroy_plan(p1);

#elif defined __FFTW2

//std::cout << "\n Call FFTW2 ";

p1 = fftw_create_plan(FFT_NR, FFTW_FORWARD, FFTW_ESTIMATE);

fftw_one(p1, fft_phir, fft_phik);

//fftw_destroy_plan(p1);

#endif

double dk_uniform = ModuleBase::TWO_PI / FFT_NR / dr;

//for (int ik = 0; ik < FFT_NR/2; ik++)

//{

//double kp = ik * dk_uniform;

//fft_ndphik[ik].re = pow(kp, nd) * fft_phik[ik].re;

//fft_ndphik[ik].im = 0.0;

//}

//for (int ik = FFT_NR/2; ik < FFT_NR; ik++)

//{

//double kp = -(FFT_NR - ik)* dk_uniform;

//fft_ndphik[ik].re = pow(kp, nd) * fft_phik[ik].re;

//fft_ndphik[ik].im = 0.0;

//}

for (int ik = 0; ik < FFT_NR/2; ik++) // ik = 0 1 ... 8008

{

double kp = ik * dk_uniform;

c_re(fft_ndphik[ik]) = pow(kp, nd) * c_re(fft_phik[ik]);

c_im(fft_ndphik[ik]) = 0.0;

}

for (int ik = FFT_NR/2; ik < FFT_NR; ik++) // ik = 8009 8010 ... 16018

{

//double kp = -(FFT_NR - ik )* dk_uniform; //(...) = 8010 8009 ... 1

double kp = -(FFT_NR - ik -1)* dk_uniform; //(...) -> 8009 8008 ... 0

c_re(fft_ndphik[ik]) = pow(kp, nd) * c_re(fft_phik[ik]);

c_im(fft_ndphik[ik]) = 0.0;

}

#if defined __FFTW3

p2 = fftw_plan_dft_1d(FFT_NR, fft_ndphik, fft_ndphir, FFTW_BACKWARD, FFTW_ESTIMATE);

fftw_execute(p2);

//fftw_destroy_plan(p2);

#elif defined __FFTW2

p2 = fftw_create_plan(FFT_NR, FFTW_BACKWARD, FFTW_ESTIMATE);

fftw_one(p2, fft_ndphik, fft_ndphir);

//fftw_destroy_plan(p2);

#endif

bool is_re = true;

double fac = 0.0;

if (nd % 4 == 0)

{

is_re = true;

fac = 1.0;

}

else if (nd % 4 == 1)

{

is_re = false;

fac = -1.0;

}

else if (nd % 4 == 2)

{

is_re = true;

fac = -1.0;

}

else

{

is_re = false;

fac = 1.0;

}

for (int ir = 0; ir < mesh; ir++)

{

if (is_re)

{

phind[ir] = fac * c_re(fft_ndphir[ir]) / FFT_NR;

}

else

{

phind[ir] = fac * c_im(fft_ndphir[ir]) / FFT_NR;

}

}

fftw_destroy_plan (p1);

fftw_destroy_plan (p2);

delete [] fft_phir;

delete [] fft_phik;

delete [] fft_ndphik;

delete [] fft_ndphir;

ModuleBase::timer::end("Mathzone_Add1", "Uni_Deriv_Phi");