docs: def _get_corr_stddevs(C, tau_ss, stddev_types, num_sites, phi_ss, NL=None,
                        tau_value=None):
      """
      Return standard deviations adjusted for single station sigma
      as the total standard deviation - as proposed to be used in
      the Swiss Hazard Model [2014].
      """
      stddevs = []
      temp_stddev = phi_ss * phi_ss
  
      if tau_value is not None and NL is not None:
          temp_stddev = temp_stddev + tau_value * tau_value * ((1 + NL) ** 2)
      else:
          temp_stddev = temp_stddev + C[tau_ss] * C[tau_ss]
  
      for stddev_type in stddev_types:
          if stddev_type == const.StdDev.TOTAL:
              stddevs.append(np.sqrt(temp_stddev) + np.zeros(num_sites))
      return stddevs