content: import numpy as np
  import qutip
  
  def qubit_integrate(w, theta, gamma1, gamma2, psi0, tlist):
      # operators and the hamiltonian
      sx = qutip.sigmax()
      sy = qutip.sigmay()
      sz = qutip.sigmaz()
      sm = qutip.sigmam()
      H = w * (np.cos(theta) * sz + np.sin(theta) * sx)
      # collapse operators
      c_op_list = []
      n_th = 0.5 # temperature
      rate = gamma1 * (n_th + 1)
      if rate > 0.0: c_op_list.append(np.sqrt(rate) * sm)
      rate = gamma1 * n_th
      if rate > 0.0: c_op_list.append(np.sqrt(rate) * sm.dag())
      rate = gamma2
      if rate > 0.0: c_op_list.append(np.sqrt(rate) * sz)
      # evolve and calculate expectation values
      output = qutip.mesolve(H, psi0, tlist, c_op_list, [sx, sy, sz])
      return output.expect[0], output.expect[1], output.expect[2]
  
  ## calculate the dynamics
  w     = 1.0 * 2 * np.pi  # qubit angular frequency
  theta = 0.2 * np.pi      # qubit angle from sigma_z axis (toward sigma_x axis)
  gamma1 = 0.5             # qubit relaxation rate
  gamma2 = 0.2             # qubit dephasing rate
  # initial state
  a = 1.0
  psi0 = (a*qutip.basis(2, 0) + (1-a)*qutip.basis(2, 1))/np.sqrt(a**2 + (1-a)**2)
  tlist = np.linspace(0, 4, 250)
  #expectation values for ploting
  sx, sy, sz = qubit_integrate(w, theta, gamma1, gamma2, psi0, tlist)
  