docs: def LOS_get_sample(D, u, dL, DL=None, dLMode='abs', method='sum', Test=True):
      """ Return the sampled line, with the specified method
  
      'linspace': return the N+1 edges, including the first and last point
      'sum' : return the N middle of the segments
      'simps': return the N+1 egdes, where N has to be even (scipy.simpson requires an even number of intervals)
      'romb' : return the N+1 edges, where N+1 = 2**k+1 (fed to scipy.romb for integration)
      """
      if Test:
          assert all([type(dd) is np.ndarray and dd.shape==(3,) for dd in [D,u]])
          assert not hasattr(dL,'__iter__')
          assert DL is None or all([hasattr(DL,'__iter__'), len(DL)==2, all([not hasattr(dd,'__iter__') for dd in DL])])
          assert dLMode in ['abs','rel']
          assert type(method) is str and method in ['linspace','sum','simps','romb']
      # Compute the minimum number of intervals to satisfy the specified resolution
      N = int(np.ceil((DL[1]-DL[0])/dL)) if dLMode=='abs' else int(np.ceil(1./dL))
      # Modify N according to the desired method
      if method=='simps':
          N = N if N%2==0 else N+1
      elif method=='romb':
          N = 2**int(np.ceil(np.log(N)/np.log(2.)))
  
      # Derive k and dLr
      if method=='sum':
          dLr = (DL[1]-DL[0])/N
          k = DL[0] + (0.5+np.arange(0,N))*dLr
      else:
          k, dLr = np.linspace(DL[0], DL[1], N+1, endpoint=True, retstep=True, dtype=float)
  
      Pts = D[:,np.newaxis] + k[np.newaxis,:]*u[:,np.newaxis]
      return Pts, k, dLr