docs: def track2(lat1: float, lon1: float, lat2: float, lon2: float,
             ell: Ellipsoid = None, npts: int = 100, deg: bool = True):
      """
      computes great circle tracks starting at the point lat1, lon1 and ending at lat2, lon2
  
      Parameters
      ----------
  
      Lat1 : float or numpy.ndarray of float
          Geodetic latitude of first point (degrees)
      Lon1 : float or numpy.ndarray of float
          Geodetic longitude of first point (degrees)
      Lat2 : float or numpy.ndarray of float
          Geodetic latitude of second point (degrees)
      Lon2 : float or numpy.ndarray of float
          Geodetic longitude of second point (degrees)
      ell : Ellipsoid, optional
            reference ellipsoid
      npts : int, optional
          number of points (default is 100)
      deg : bool, optional
          degrees input/output  (False: radians in/out)
  
      Results
      -------
  
      lats : numpy.ndarray of float
          latitudes of points along track
      lons : numpy.ndarray of float
          longitudes of points along track
  
      Based on code posted to the GMT mailing list in Dec 1999 by Jim Levens and by Jeff Whitaker <jeffrey.s.whitaker@noaa.gov>
      """
  
      if ell is None:
          ell = Ellipsoid()
  
      if npts <= 1:
          raise ValueError('npts must be greater than 1')
  
      if npts == 2:
          return [lat1, lat2], [lon1, lon2]
  
      if deg is True:
          rlat1, rlon1, rlat2, rlon2 = np.radians([lat1, lon1, lat2, lon2])
      else:
          rlat1, rlon1, rlat2, rlon2 = lat1, lon1, lat2, lon2
  
      gcarclen = 2. * np.arcsin(np.sqrt((np.sin((rlat1 - rlat2) / 2))**2 +
                                        np.cos(rlat1) * np.cos(rlat2) * (np.sin((rlon1 - rlon2) / 2))**2))
      # check to see if points are antipodal (if so, route is undefined).
      if np.allclose(gcarclen, pi):
          raise ValueError('cannot compute intermediate points on a great circle whose endpoints are antipodal')
  
      distance, azimuth, _ = vdist(lat1, lon1, lat2, lon2)
      incdist = distance / (npts - 1)
  
      latpt = lat1
      lonpt = lon1
      lons = [lonpt]
      lats = [latpt]
      for n in range(npts - 2):
          latptnew, lonptnew, _ = vreckon(latpt, lonpt, incdist, azimuth)
          _, azimuth, _ = vdist(latptnew, lonptnew, lat2, lon2, ell=ell)
          lats.append(latptnew)
          lons.append(lonptnew)
          latpt = latptnew
          lonpt = lonptnew
      lons.append(lon2)
      lats.append(lat2)
  
      if not deg:
          lats = np.radians(lats)
          lons = np.radians(lons)
  
      return lats, lons