code: def burst_range_spectrum(psd, snr=8, energy=1e-2):
      
      # calculate frequency dependent range in parsecs
      a = (constants.G * energy * constants.M_sun * 0.4 /
           (pi**2 * constants.c))**(1/2.)
      dspec = psd ** (-1/2.) * a / (snr * psd.frequencies)
  
      # convert to output unit
      rspec = dspec.to('Mpc')
  
      # rescale 0 Hertz (which has 0 range always)
      if rspec.f0.value == 0.0:
          rspec[0] = 0.0
  
      return rspec
docs: Calculate the frequency-dependent GW burst range from a strain PSD
  
      Parameters
      ----------
      psd : `~gwpy.frequencyseries.FrequencySeries`
          the instrumental power-spectral-density data
  
      snr : `float`, optional
          the signal-to-noise ratio for which to calculate range,
          default: `8`
  
      energy : `float`, optional
          the relative energy output of the GW burst,
          default: `0.01` (GRB-like burst)
  
      Returns
      -------
      rangespec : `~gwpy.frequencyseries.FrequencySeries`
          the burst range `FrequencySeries` [Mpc (default)]