code: def funTransChica(imagl, idril, ibfield, gamma0, xoy='x'):
      
      m0 = 9.10938215e-31
      e0 = 1.602176487e-19
      c0 = 299792458
      rho = np.sqrt(gamma0 ** 2 - 1) * m0 * c0 / ibfield / e0
      theta = np.arcsin(imagl / rho)
      ld = idril
      mx = reduce(np.dot, [funTransDrift(idril),
                           funTransSectX(theta, rho), funTransEdgeX(theta, rho),
                           funTransDrift(ld),
                           funTransEdgeX(-theta, -rho), funTransSectX(-theta, -rho),
                           funTransDrift(ld),
                           funTransSectX(-theta, -rho), funTransEdgeX(-theta, -rho),
                           funTransDrift(ld),
                           funTransEdgeX(theta, rho), funTransSectX(theta, rho),
                           funTransDrift(idril)])
      my = reduce(np.dot, [funTransDrift(idril),
                           funTransSectY(theta, rho), funTransEdgeY(theta, rho),
                           funTransDrift(ld),
                           funTransEdgeY(-theta, -rho), funTransSectY(-theta, -rho),
                           funTransDrift(ld),
                           funTransSectY(-theta, -rho), funTransEdgeY(-theta, -rho),
                           funTransDrift(ld),
                           funTransEdgeY(theta, rho), funTransSectY(theta, rho),
                           funTransDrift(idril)])
      if xoy == 'x':
          m = mx
      else:
          m = my
  
      return m
docs: Chicane matrix, composed of four rbends, seperated by drifts
  
      :param imagl: rbend width, in [m]
      :param idril: drift length between two adjacent rbends, in [m]
      :param ibfield: rbend magnetic strength, in [T]
      :param gamma0: electron energy, gamma
      :param xoy: ``'x'`` or ``'y'``, matrix in X or Y direction, ``'x'`` by default
      :return: 2x2 numpy array