def METHOD_NAME(mean_stds, probs, rho, imti, imls, cs_poes,
            phi_b, invtime, c, _c=None):
    # `mean_stds` is an array with the values of mean of log(gm) and the
    # corresponding standard deviation. The size of this array is 4 x M x U
    # where M is the number if IMLs and U is the number of ruptures considered.
    # `probs` is an array with the probabilities of at least one occurrence of
    # a given rupture. `rho` is a vector of size M with the correlation
    # coefficients. `imti` is the index of the conditioning IMT. `imls` is a
    # list with the IMLs (corresponding to different probabilities of
    # exceedance) for the conditioning IMT. `cs_poes` are the probabilities of
    # exceedance characterising the values in `imls`. `phi_b` is float.
    # `invtime` is the investigation time [yr]. `c` has shape M x ? x ?
    M, O, _ = c.shape
    mu = mean_stds[0]  # shape (M, U)
    sig = mean_stds[1]  # shape (M, U)
    for p, iml in enumerate(imls):
        # Calculate the contribution of each rupture to the total
        # probability of occurrence for the reference IMT. Both
        # `eps` and `poes` have shape U
        eps = (numpy.log(iml) - mu[imti]) / sig[imti]
        poes = truncnorm_sf(phi_b, eps)
        # Converting to rates and dividing by the rate of exceedance of the
        # reference IMT and level.  This is eq. 13 of the OQ Engine Underlying
        # Hazard Science Book.
        ws = -numpy.log((1. - probs) ** poes) / invtime
        # Normalizing by the AfE for the investigated IMT and level
        ws /= -numpy.log(1. - cs_poes[p]) / invtime
        # weights not summing up to 1
        c[:, 0, p] = ws.sum()
        # For each intensity measure type
        for m in range(len(mu)):
            # Equation 14 in Lin et al. (2013)
            term1 = mu[m] + rho[m] * eps * sig[m]
            c[m, 1, p] = ws @ term1
            # This is executed only if we already have the final CS
            if _c is not None:
                # Equation 15 in Lin et al. (2013)
                term2 = sig[m] * (1. - rho[m]**2)**0.5
                term3 = term1 - _c[m, 1, p]
                c[m, 2, p] = ws @ (term2**2 + term3**2)