code: def _direct_nbody_force(q,m,t,pot,softening,softening_args):
      
      #First do the particles
      #Calculate all the distances
      nq= len(q)
      dim= len(q[0])
      dist_vec= nu.zeros((nq,nq,dim))
      dist= nu.zeros((nq,nq))
      for ii in range(nq):
          for jj in range(ii+1,nq):
              dist_vec[ii,jj,:]= q[jj]-q[ii]
              dist_vec[jj,ii,:]= -dist_vec[ii,jj,:]
              dist[ii,jj]= linalg.norm(dist_vec[ii,jj,:])
              dist[jj,ii]= dist[ii,jj]
      #Calculate all the forces
      force= []
      for ii in range(nq):
          thisforce= nu.zeros(dim)
          for jj in range(nq):
              if ii == jj: continue
              thisforce+= m[jj]*softening(dist[ii,jj],*softening_args)\
                  /dist[ii,jj]*dist_vec[ii,jj,:]
          force.append(thisforce)
      #Then add the external force
      if pot is None: return force
      for ii in range(nq):
          force[ii]+= _external_force(q[ii],t,pot)
      return force
docs: Calculate the force