docs: def csnr(freqs, hc, hn, fmrg, fpeak, prefactor=1.0):
      """Calculate the SNR of a frequency domain waveform.
  
      SNRCalculation is a function that takes waveforms (frequencies and hcs)
      and a noise curve, and returns SNRs for all binary phases and the whole waveform.
  
      Arguments:
          freqs (1D or 2D array of floats): Frequencies corresponding to the waveforms.
              Shape is (num binaries, num_points) if 2D.
              Shape is (num_points,) if 1D for one binary.
          hc (1D or 2D array of floats): Characteristic strain of the waveforms.
              Shape is (num binaries, num_points) if 2D.
              Shape is (num_points,) if 1D for one binary.
          fmrg: (scalar float or 1D array of floats): Merger frequency of each binary separating
              inspiral from merger phase. (0.014/M) Shape is (num binaries,)
              if more than one binary.
          fpeak: (scalar float or 1D array of floats): Peak frequency of each binary separating
              merger from ringdown phase. (0.014/M) Shape is (num binaries,)
              if more than one binary.
          hn: (1D or 2D array of floats): Characteristic strain of the noise.
              Shape is (num binaries, num_points) if 2D.
              Shape is (num_points,) if 1D for one binary.
          prefactor (float, optional): Factor to multiply snr (not snr^2) integral values by.
              Default is 1.0.
  
      Returns:
          (dict): Dictionary with SNRs from each phase.
  
      """
      cfd = os.path.dirname(os.path.abspath(__file__))
      if 'phenomd.cpython-35m-darwin.so' in os.listdir(cfd):
          exec_call = cfd + '/phenomd.cpython-35m-darwin.so'
  
      else:
          exec_call = cfd + '/phenomd/phenomd.so'
  
      c_obj = ctypes.CDLL(exec_call)
  
      # check dimensionality
      remove_axis = False
      try:
          len(fmrg)
      except TypeError:
          remove_axis = True
          freqs, hc = np.array([freqs]), np.array([hc])
          hn, fmrg, fpeak = np.array([hn]), np.array([fmrg]), np.array([fpeak])
  
      # this implimentation in ctypes works with 1D arrays
      freqs_in = freqs.flatten()
      hc_in = hc.flatten()
      hn_in = hn.flatten()
  
      num_binaries, length_of_signal = hc.shape
  
      # prepare outout arrays
      snr_cast = ctypes.c_double*num_binaries
      snr_all = snr_cast()
      snr_ins = snr_cast()
      snr_mrg = snr_cast()
      snr_rd = snr_cast()
  
      # find SNR values
      c_obj.SNR_function(ctypes.byref(snr_all), ctypes.byref(snr_ins),
                         ctypes.byref(snr_mrg), ctypes.byref(snr_rd),
                         freqs_in.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
                         hc_in.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
                         hn_in.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
                         fmrg.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
                         fpeak.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
                         ctypes.c_int(length_of_signal), ctypes.c_int(num_binaries))
  
      # make into numpy arrays
      snr_all, snr_ins, = np.ctypeslib.as_array(snr_all), np.ctypeslib.as_array(snr_ins)
      snr_mrg, snr_rd = np.ctypeslib.as_array(snr_mrg), np.ctypeslib.as_array(snr_rd)
  
      # remove axis if one binary
      if remove_axis:
          snr_all, snr_ins, snr_mrg, snr_rd = snr_all[0], snr_ins[0], snr_mrg[0], snr_rd[0]
  
      # prepare output by multiplying by prefactor
      return ({'all': snr_all*prefactor, 'ins': snr_ins*prefactor,
              'mrg': snr_mrg*prefactor, 'rd': snr_rd*prefactor})