code: def rdf(coords_a, coords_b, binsize=0.002,
          cutoff=1.5, periodic=None, normalize=True):
      
      
      period = periodic[0, 0], periodic[1,1], periodic[2,2]
      distances = distances_within(coords_a, coords_b, cutoff,
                                   np.array(period, dtype=np.double))
      
      n_a = len(coords_a)
      n_b = len(coords_b)
  
      volume = periodic[0, 0] * periodic[1, 1] * periodic[2, 2]
  
      int_distances = np.rint(distances/binsize).astype(int)
      hist = np.bincount(int_distances)
      
      bin_edges = np.arange(len(hist)+1) * binsize
          
      
      if normalize:
          dr  = binsize
          normfac = volume/(n_a*n_b)
  
          # Normalize this by a sphere shell
          for i, r in enumerate(bin_edges[1:]):
              hist[i] /= ((4.0/3.0 * np.pi * (r + 0.5*dr)**3)
                          - (4.0/3.0 * np.pi * (r- 0.5*dr)**3))
  
          # Normalize by density
          hist = hist * normfac
  
      # Cutting up to rmax value
          
      width = cutoff/binsize + 1
      return bin_edges[0:width], hist[0:width]
docs: Calculate the radial distribution function of *coords_a* against
      *coords_b*.
  
      **Parameters**
      - coords_a: np.ndarray((3, NA))
          first set of coordinates
      - coords_b: np.ndarray((3, NB))
          coordinates to calculate the RDF against
      - periodic: np.ndarray((3, 3)) or None
          Wether or not include periodic images in the calculation
      - normalize: True or False
          gromacs-like normalization
      - cutoff: 
          where to cutoff the RDF