docs: def npoints_towards(lon, lat, depth, azimuth, hdist, vdist, npoints):
      """
      Find a list of specified number of points starting from a given one
      along a great circle arc with a given azimuth measured in a given point.
  
      :param float lon, lat, depth:
          Coordinates of a point to start from. The first point in a resulting
          list has these coordinates.
      :param azimuth:
          A direction representing a great circle arc together with a reference
          point.
      :param hdist:
          Horizontal (geodetic) distance from reference point to the last point
          of the resulting list, in km.
      :param vdist:
          Vertical (depth) distance between reference and the last point, in km.
      :param npoints:
          Integer number of points to return. First and last points count,
          so if there have to be two intervals, ``npoints`` should be 3.
      :returns:
          Tuple of three 1d numpy arrays: longitudes, latitudes and depths
          of resulting points respectively.
  
      Implements "completely general but more complicated algorithm" from
      http://williams.best.vwh.net/avform.htm#LL
      """
      assert npoints > 1
      rlon, rlat = numpy.radians(lon), numpy.radians(lat)
      tc = numpy.radians(360 - azimuth)
      hdists = numpy.arange(npoints, dtype=float)
      hdists *= (hdist / EARTH_RADIUS) / (npoints - 1)
      vdists = numpy.arange(npoints, dtype=float)
      vdists *= vdist / (npoints - 1)
  
      sin_dists = numpy.sin(hdists)
      cos_dists = numpy.cos(hdists)
      sin_lat = numpy.sin(rlat)
      cos_lat = numpy.cos(rlat)
  
      sin_lats = sin_lat * cos_dists + cos_lat * sin_dists * numpy.cos(tc)
      lats = numpy.degrees(numpy.arcsin(sin_lats))
  
      dlon = numpy.arctan2(numpy.sin(tc) * sin_dists * cos_lat,
                           cos_dists - sin_lat * sin_lats)
      lons = numpy.mod(rlon - dlon + numpy.pi, 2 * numpy.pi) - numpy.pi
      lons = numpy.degrees(lons)
  
      depths = vdists + depth
  
      # the first point should be left intact
      lons[0] = lon
      lats[0] = lat
      depths[0] = depth
  
      return lons, lats, depths