#include <math.h>

double userfunc_(double *x)
{
  return 0.0;
}


double userfunc1_(double *x)
{
  /* coulomb force */
  return pow(*x,-1.5);
}


double userfunc2_(double *x)
{
  /* coulomb potential */
  return pow(*x,-0.5);
}


double userfunc3_(double *x)
{
  /* LJ force */
  return 2.0*pow(*x,-7.0) - pow(*x,-4.0);
}


double userfunc4_(double *x)
{
  /* LJ potential */
  return pow(*x,-6.0) - pow(*x,-3.0);
}


double userfunc5_(double *x)
{
  return sin(M_PI*(*x));
}


double userfunc6_(double *x)
{
  return cos(M_PI*(*x));
}


double userfunc7_(double *x)
{
  /* Ewald real space force */
  return M_2_SQRTPI*exp(-*x)*pow(*x,-1.0)
    + erfc(sqrt(*x))*pow(*x,-1.5);
}


double userfunc8_(double *x)
{
  /* Ewald real space potential */
  return erfc(sqrt(*x))*pow(*x,-0.5);
}


double userfunc9_(double *x)
{
  return sin(M_PI*(*x));
}


double userfunc10_(double *x)
{
  return 0.0;
}


double userfunc11_(double *x)
{
  /* Ewald real space force with cutoff */
  double rcut2,ret;

  rcut2=32.0;
  if(*x<rcut2){
      ret=M_2_SQRTPI*exp(-*x)*pow(*x,-1.0)
	  + erfc(sqrt(*x))*pow(*x,-1.5);
  }
  else{
      ret=0.0;
  }
  return ret;
}


double userfunc12_(double *x)
{
  double ret,rcut,ccut;

  rcut=64.0;
  ccut=3.0*rcut/5.0;
  if(*x<ccut)      ret=pow(*x,-1.5);
  else if(*x<rcut) ret=-1.5*pow(ccut,-2.5)*(*x-ccut)+pow(ccut,-1.5);
  else             ret=0.0;
  
  return ret;
}