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 = (dists.rx >= 0)
          # aftershock flag. always zero since we only consider main shock
          AS = 0
  
          ln_y_ref = (
              # first line of eq. 13a
              C['c1']
              + (C['c1a'] * Frv
                 + C['c1b'] * Fnm
                 + C['c7'] * (rup.ztor - 4))
              * (1 - AS)
              + (C['c10'] + C['c7a'] * (rup.ztor - 4)) * AS
              # second line
              + C['c2'] * (rup.mag - 6)
              + ((C['c2'] - C['c3']) / C['cn'])
              * np.log(1 + np.exp(C['cn'] * (C['cm'] - rup.mag)))
              # third line
              + C['c4']
              * np.log(dists.rrup
                       + C['c5']
                       * np.cosh(C['c6'] * max(rup.mag - C['chm'], 0)))
              # fourth line
              + (C['c4a'] - C['c4'])
              * np.log(np.sqrt(dists.rrup ** 2 + C['crb'] ** 2))
              # fifth line
              + (C['cg1'] + C['cg2'] / (np.cosh(max(rup.mag - C['cg3'], 0))))
              * dists.rrup
              # sixth line
              + C['c9'] * Fhw
              * np.tanh(dists.rx
                        * (np.cos(np.radians(rup.dip)) ** 2)
                        / C['c9a'])
              * (1 - np.sqrt(dists.rjb ** 2 + rup.ztor ** 2)
                 / (dists.rrup + 0.001))
          )
          return ln_y_ref