docs: def intervals_between(lon1, lat1, depth1, lon2, lat2, depth2, length):
      """
      Find a list of points between two given ones that lie on the same
      great circle arc and are equally spaced by ``length`` km.
  
      :param float lon1, lat1, depth1:
          Coordinates of a point to start placing intervals from. The first
          point in the resulting list has these coordinates.
      :param float lon2, lat2, depth2:
          Coordinates of the other end of the great circle arc segment
          to put intervals on. The last resulting point might be closer
          to the first reference point than the second one or further,
          since the number of segments is taken as rounded division of
          length between two reference points and ``length``.
      :param length:
          Required distance between two subsequent resulting points, in km.
      :returns:
          Tuple of three 1d numpy arrays: longitudes, latitudes and depths
          of resulting points respectively.
  
      Rounds the distance between two reference points with respect
      to ``length`` and calls :func:`npoints_towards`.
      """
      assert length > 0
      hdist = geodetic_distance(lon1, lat1, lon2, lat2)
      vdist = depth2 - depth1
      # if this method is called multiple times with coordinates that are
      # separated by the same distance, because of floating point imprecisions
      # the total distance may have slightly different values (for instance if
      # the distance between two set of points is 65 km, total distance can be
      # 64.9999999999989910 and 65.0000000000020322). These two values bring to
      # two different values of num_intervals (32 in the first case and 33 in
      # the second), and this is a problem because for the same distance we
      # should have the same number of intervals. To reduce potential differences
      # due to floating point errors, we therefore round total_distance to a
      # fixed precision (7)
      total_distance = round(numpy.sqrt(hdist ** 2 + vdist ** 2), 7)
      num_intervals = int(round(total_distance / length))
      if num_intervals == 0:
          return numpy.array([lon1]), numpy.array([lat1]), numpy.array([depth1])
      dist_factor = (length * num_intervals) / total_distance
      return npoints_towards(
          lon1, lat1, depth1, azimuth(lon1, lat1, lon2, lat2),
          hdist * dist_factor, vdist * dist_factor, num_intervals + 1)