Scalar Field Potential

/**
 * Compute the scalar field potential and all its derivatives up to
 * the requested order in a single function call. This avoids redundant
 * extraction of c1–c4, redundant switch dispatch, and allows sharing
 * common subexpressions (exp, pow, trig) across derivatives.
 *
 * @param pba        Input: pointer to background structure
 * @param phi        Input: scalar field value
 * @param V          Output: V(phi)
 * @param dV         Output: V'(phi) (pure, without coupling)
 * @param ddV        Output: V''(phi)
 * @param d3V        Output: V'''(phi) (only written if d3V != NULL)
 * @param d4V        Output: V''''(phi) (only written if d4V != NULL)
 */
void V_scf_derivs(
    struct background *pba,
    double phi,
    double *V,
    double *dV,
    double *ddV,
    double *d3V,
    double *d4V)
{
  double c1 = pba->scf_parameters[0];
  double c2 = pba->scf_parameters[1];
  double c3 = pba->scf_parameters[2];
  double c4 = pba->scf_parameters[3];

  switch (pba->scf_potential)
  {
  case 1: /* power-law: V = c1^(4-c2) * phi^c2 + c3 */
  {
    double pref = pow(c1, 4. - c2);
    double phic2 = pow(phi, c2);
    double phic2m1 = phic2 / phi;   /* phi^(c2-1) */
    double phic2m2 = phic2m1 / phi; /* phi^(c2-2) */
    *V = pref * phic2 + c3;
    *dV = c2 * pref * phic2m1;
    *ddV = c2 * (c2 - 1.) * pref * phic2m2;
    if (d3V)
      *d3V = (c2 - 2.) * (c2 - 1.) * c2 * pref * phic2m2 / phi;
    if (d4V)
      *d4V = (c2 - 3.) * (c2 - 2.) * (c2 - 1.) * c2 * pref * phic2m2 / (phi * phi);
    break;
  }
  case 2: /* cosine: V = c1 * cos(c2*phi) */
  {
    double s = sin(c2 * phi);
    double c = cos(c2 * phi);
    *V = c1 * c;
    *dV = -c1 * c2 * s;
    *ddV = -c1 * c2 * c2 * c;
    if (d3V)
      *d3V = c1 * c2 * c2 * c2 * s;
    if (d4V)
      *d4V = c1 * c2 * c2 * c2 * c2 * c;
    break;
  }
  case 3: /* hyperbolic: V = c1 * [1 - tanh(c2*phi)] */
  {
    double t = tanh(c2 * phi);
    double ch = cosh(c2 * phi);
    double sech2 = 1.0 / (ch * ch);
    *V = c1 * (1. - t);
    *dV = -c1 * c2 * sech2;
    *ddV = 2. * c1 * c2 * c2 * t * sech2;
    if (d3V)
      *d3V = -2. * c1 * c2 * c2 * c2 * (cosh(2. * c2 * phi) - 2.) * sech2 * sech2;
    if (d4V)
    {
      double ch5 = ch * ch * ch * ch * ch;
      *d4V = 2. * c1 * c2 * c2 * c2 * c2 * (sinh(3. * c2 * phi) - 11. * sinh(c2 * phi)) / ch5;
    }
    break;
  }
  case 4: /* pNG: V = c1^4 * [1 + cos(phi/c2)] */
  {
    double c1_4 = pow(c1, 4.);
    double invc2 = 1.0 / c2;
    double s = sin(phi * invc2);
    double c = cos(phi * invc2);
    *V = c1_4 * (1. + c);
    *dV = -c1_4 * invc2 * s;
    *ddV = -c1_4 * invc2 * invc2 * c;
    if (d3V)
      *d3V = c1_4 * invc2 * invc2 * invc2 * s;
    if (d4V)
      *d4V = c1_4 * invc2 * invc2 * invc2 * invc2 * c;
    break;
  }
  case 5: /* iPL: V = c1^(4+c2) * phi^(-c2) */
  {
    double pref = pow(c1, 4. + c2);
    double phimc2 = pow(phi, -c2);
    double inv_phi = 1.0 / phi;
    *V = pref * phimc2;
    *dV = -c2 * pref * phimc2 * inv_phi;
    *ddV = c2 * (c2 + 1.) * pref * phimc2 * inv_phi * inv_phi;
    if (d3V)
      *d3V = -c2 * (c2 + 1.) * (c2 + 2.) * pref * phimc2 * inv_phi * inv_phi * inv_phi;
    if (d4V)
      *d4V = c2 * (c2 + 1.) * (c2 + 2.) * (c2 + 3.) * pref * phimc2 * inv_phi * inv_phi * inv_phi * inv_phi;
    break;
  }
  case 6: /* exponential: V = c1 * exp(-c2*phi) */
  {
    double e = c1 * exp(-c2 * phi); /* = V */
    *V = e;
    *dV = -c2 * e;
    *ddV = c2 * c2 * e;
    if (d3V)
      *d3V = -c2 * c2 * c2 * e;
    if (d4V)
      *d4V = c2 * c2 * c2 * c2 * e;
    break;
  }
  case 7: /* SqE: V = c1^(c2+4) * phi^(-c2) * exp(c1*phi^2) */
  {
    double pref = pow(c1, c2 + 4.);
    double phimc2 = pow(phi, -c2);
    double ep = exp(c1 * phi * phi);
    double base = pref * phimc2 * ep;
    double phi2 = phi * phi;
    double inv_phi = 1.0 / phi;
    *V = base;
    *dV = base * inv_phi * (2. * c1 * phi2 - c2);
    double phi4 = phi2 * phi2;
    *ddV = base * inv_phi * inv_phi * (4. * c1 * c1 * phi4 + (-4. * c1 * c2 - 2. * c1) * phi2 + c2 * (c2 + 1.));
    if (d3V)
    {
      double phi6 = phi4 * phi2;
      *d3V = base * inv_phi * inv_phi * inv_phi * (8. * c1 * c1 * c1 * phi6 + (-12. * c1 * c1 * (c2 - 1) * phi4) + 6. * c1 * c2 * c2 * phi2 - c2 * (c2 + 1.) * (c2 + 2.));
    }
    if (d4V)
    {
      double phi4v = phi2 * phi2;
      double phi6 = phi4v * phi2;
      double phi8 = phi6 * phi2;
      *d4V = base * inv_phi * inv_phi * inv_phi * inv_phi * (16. * c1 * c1 * c1 * c1 * phi8 + (16. * c1 * c1 * c1 * phi6 * (3. - 2. * c2)) + (12. * c1 * c1 * phi4v * (2. * (c2 - 1.) * c2 + 1.)) - 4 * c1 * phi2 * c2 * (c2 + 1.) * (2. * c2 + 1.) + c2 * (c2 + 1.) * (c2 + 2.) * (c2 + 3.));
    }
    break;
  }
  case 8: /* Bean: V = c1 * [(c4-phi)^2 + c2] * exp(-c3*phi) */
  {
    double e = exp(-c3 * phi);
    double dp = c4 - phi;
    double dp2 = dp * dp;
    double bracket = dp2 + c2;
    *V = c1 * bracket * e;
    *dV = -c1 * e * (2. * dp + c3 * bracket);
    *ddV = c1 * e * (2. - 4. * dp * c3 + c3 * c3 * bracket);
    if (d3V)
      *d3V = -c1 * c3 * e * (c3 * (c3 * (c2 + phi * phi) - 6. * phi) + 6.);
    if (d4V)
      *d4V = c1 * c3 * c3 * e * (c3 * (c3 * (c2 + phi * phi) - 8. * phi) + 12.);
    break;
  }
  case 9: /* DoubleExp: V = c1 * (exp(-c2*phi) + c3 * exp(-c4*phi)) */
  {
    double e2 = exp(-c2 * phi);
    double e4 = exp(-c4 * phi);
    *V = c1 * (e2 + c3 * e4);
    *dV = -c1 * (c2 * e2 + c3 * c4 * e4);
    *ddV = c1 * (c2 * c2 * e2 + c3 * c4 * c4 * e4);
    if (d3V)
      *d3V = -c1 * (c2 * c2 * c2 * e2 + c3 * c4 * c4 * c4 * e4);
    if (d4V)
      *d4V = c1 * (c2 * c2 * c2 * c2 * e2 + c3 * c4 * c4 * c4 * c4 * e4);
    break;
  }
  default:
    *V = 0.;
    *dV = 0.;
    *ddV = 0.;
    if (d3V)
      *d3V = 0.;
    if (d4V)
      *d4V = 0.;
    break;
  }
}