docs: def _get_ln_y_ref(self, rup, dists, C):
          """
          Get an intensity on a reference soil.
  
          Implements eq. 13a.
          """
          # reverse faulting flag
          Frv = 1. if 30 <= rup.rake <= 150 else 0.
          # normal faulting flag
          Fnm = 1. if -120 <= rup.rake <= -60 else 0.
          # hanging wall flag
  
          Fhw = np.zeros_like(dists.rx)
          idx = np.nonzero(dists.rx >= 0.)
          Fhw[idx] = 1.
  
          # a part in eq. 11
          mag_test1 = np.cosh(2. * max(rup.mag - 4.5, 0))
  
          # centered DPP
          centered_dpp = self._get_centered_cdpp(dists)
          # centered_ztor
          centered_ztor = self._get_centered_ztor(rup, Frv)
          #
          dist_taper = np.fmax(1 - (np.fmax(dists.rrup - 40,
                                    np.zeros_like(dists)) / 30.),
                               np.zeros_like(dists))
          dist_taper = dist_taper.astype(np.float64)
          ln_y_ref = (
              # first part of eq. 11
              C['c1']
              + (C['c1a'] + C['c1c'] / mag_test1) * Frv
              + (C['c1b'] + C['c1d'] / mag_test1) * Fnm
              + (C['c7'] + C['c7b'] / mag_test1) * centered_ztor
              + (C['c11'] + C['c11b'] / mag_test1) *
              np.cos(math.radians(rup.dip)) ** 2
              # second part
              + C['c2'] * (rup.mag - 6)
              + ((C['c2'] - C['c3']) / C['cn'])
              * np.log(1 + np.exp(C['cn'] * (C['cm'] - rup.mag)))
              # third part
              + C['c4']
              * np.log(dists.rrup + C['c5']
                       * np.cosh(C['c6'] * max(rup.mag - C['chm'], 0)))
              + (C['c4a'] - C['c4'])
              * np.log(np.sqrt(dists.rrup ** 2 + C['crb'] ** 2))
              # forth part
              + (C['cg1'] + C['cg2'] / (np.cosh(max(rup.mag - C['cg3'], 0))))
              * dists.rrup
              # fifth part
              + C['c8'] * dist_taper
              * min(max(rup.mag - 5.5, 0) / 0.8, 1.0)
              * np.exp(-1 * C['c8a'] * (rup.mag - C['c8b']) ** 2) * centered_dpp
              # sixth part
              + C['c9'] * Fhw * np.cos(math.radians(rup.dip)) *
              (C['c9a'] + (1 - C['c9a']) * np.tanh(dists.rx / C['c9b']))
              * (1 - np.sqrt(dists.rjb ** 2 + rup.ztor ** 2)
                 / (dists.rrup + 1.0))
          )
  
          return ln_y_ref