code: def transQuad(length=0.0, k1=0.0, gamma=None):
      
      m = np.eye(6, 6, dtype=np.float64)
      if length == 0.0:
          print("warning: 'length' should be a positive float number.")
      elif gamma is not None and gamma != 0.0:
          if k1 == 0:
              print("warning: 'k1' should be a positive float number.")
              m[0, 1] = m[2, 3] = 1.0
              m[4, 5] = float(length) / gamma / gamma
          else:
              sqrtk = np.sqrt(complex(k1))
              sqrtkl = sqrtk * length
              m[0, 0] = m[1, 1] = (np.cos(sqrtkl)).real
              m[0, 1] = (np.sin(sqrtkl) / sqrtk).real
              m[1, 0] = (-np.sin(sqrtkl) * sqrtk).real
              m[2, 2] = m[3, 3] = (np.cosh(sqrtkl)).real
              m[2, 3] = (np.sinh(sqrtkl) / sqrtk).real
              m[3, 2] = (-np.sinh(sqrtkl) * sqrtk).real
              m[4, 5] = float(length) / gamma / gamma
      else:
          print("warning: 'gamma' should be a positive float number.")
      return m
docs: Transport matrix of quadrupole
  
      :param length: quad width in [m]
      :param k1: quad k1 strength in [T/m]
      :param gamma: electron energy, gamma value
      :return: 6x6 numpy array