code: def dBinaryRochedr (r, D, q, F):
      
  
      r2 = (r*r).sum()
      r1 = np.sqrt(r2)
  
      return -1./r2 - q*(r1-r[0]/r1*D)/((r[0]-D)*(r[0]-D)+r[1]*r[1]+r[2]*r[2])**1.5 - q*r[0]/r1/D/D + F*F*(1+q)*(1-r[2]*r[2]/r2)*r1
docs: Computes a derivative of the potential with respect to r.
  
      @param r:      relative radius vector (3 components)
      @param D:      instantaneous separation
      @param q:      mass ratio
      @param F:      synchronicity parameter