docs: def c_Prada(z, m, h=h, Om_M=Om_M, Om_L=Om_L):
      """Concentration from c(M) relation published in Prada et al. (2012).
  
      Parameters
      ----------
      z : float or array_like
          Redshift(s) of halos.
      m : float or array_like
          Mass(es) of halos (m200 definition), in units of solar masses.
      h : float, optional
          Hubble parameter. Default is from Planck13.
      Om_M : float, optional
          Matter density parameter. Default is from Planck13.
      Om_L : float, optional
          Cosmological constant density parameter. Default is from Planck13.
  
      Returns
      ----------
      ndarray
          Concentration values (c200) for halos.
  
      Notes
      ----------
      This c(M) relation is somewhat controversial, due to its upturn in
      concentration for high masses (normally we expect concentration to
      decrease with increasing mass). See the reference below for discussion.
  
      References
      ----------
      Calculation based on results of N-body simulations presented in:
  
      F. Prada, A.A. Klypin, A.J. Cuesta, J.E. Betancort-Rijo, and J.
      Primack, "Halo concentrations in the standard Lambda cold dark matter
      cosmology," Monthly Notices of the Royal Astronomical Society, Volume
      423, Issue 4, pp. 3018-3030, 2012.
      """
  
      z, m = _check_inputs(z, m)
  
      # EQ 13
      x = (1. / (1. + z)) * (Om_L / Om_M)**(1. / 3.)
  
      # EQ 12
      intEQ12 = np.zeros(len(x))  # integral
      for i in range(len(x)):
          # v is integration variable
          temp = integrate.quad(lambda v: (v / (1 + v**3.))**(1.5), 0, x[i])
          intEQ12[i] = temp[0]
  
      Da = 2.5 * ((Om_M / Om_L)**(1. / 3.)) * (np.sqrt(1. + x**3.) /
                                               (x**(1.5))) * intEQ12
  
      # EQ 23
      y = (1.e+12) / (h * m)
      sigma = Da * (16.9 * y**0.41) / (1. + (1.102 * y**0.2) + (6.22 * y**0.333))
  
      # EQ 21 & 22 (constants)
      c0 = 3.681
      c1 = 5.033
      alpha = 6.948
      x0 = 0.424
      s0 = 1.047  # sigma_0^-1
      s1 = 1.646  # sigma_1^-1
      beta = 7.386
      x1 = 0.526
  
      # EQ 19 & 20
      cmin = c0 + (c1 - c0) * ((1. / np.pi) * np.arctan(alpha * (x - x0)) + 0.5)
      smin = s0 + (s1 - s0) * ((1. / np.pi) * np.arctan(beta * (x - x1)) + 0.5)
  
      # EQ 18
      cmin1393 = c0 + (c1 - c0) * ((1. / np.pi) * np.arctan(alpha *
                                                            (1.393 - x0)) + 0.5)
      smin1393 = s0 + (s1 - s0) * ((1. / np.pi) * np.arctan(beta *
                                                            (1.393 - x1)) + 0.5)
      B0 = cmin / cmin1393
      B1 = smin / smin1393
  
      # EQ 15
      sigma_prime = B1 * sigma
  
      # EQ 17
      A = 2.881
      b = 1.257
      c = 1.022
      d = 0.06
  
      # EQ 16
      Cs = A * ((sigma_prime / b)**c + 1.) * np.exp(d / (sigma_prime**2.))
  
      # EQ 14
      concentration = B0 * Cs
  
      return concentration