docs: def spin_gen(particles, index, gauge=1):
      """Generates the generic spin operator in z basis for a system of
          N=particles and for the selected spin index name. where index=0..N-1
          The gauge term sets the behavoir for a system away from half-filling"""
      mat = np.zeros((2**particles, 2**particles))
  
      flipper = 2**index
      for i in range(2**particles):
          ispin = btest(i, index)
          if ispin == 1:
              mat[i ^ flipper, i] = 1
          else:
              mat[i ^ flipper, i] = gauge
      return mat