Tutorial

This tutorial will teach you the basics to write your first scripts. It does show you

  • How to set up a simple spin system

  • How to use the classes provided by teacups

  • How to start a teacups simulation

  • Where to find the outputs

  • A possible workflow for determining a plausible relaxation mechanism

It doesn’t show you how to

  • program using Python3 (here are Python3 tutorials)

  • use all features of the routine. See the sections about the Spin System, Experimental and Simulation Options for further information about possible attributes.

The main simulation function

The simulation is done by the main simulation function. You can use the function by importing the simulations module from the teacups package:

import teacups.simulations as sim

Then it can be called in your script:

spec = sim.teacups(Sys, Exp, SimOpt)

Input classes provided by teacups

The main simulation function needs three classes as input. In your script theses classes need to be set up. Their attributes define the spin system parameters (Sys), the experimental parameters (Exp) and the simulation options (SimOpt). You can import the classes form the classes module by:

import teacups.classes as cl

Sys = cl.Sys()
Exp = cl.Exp()
SimOpt = cl.SimOpt()

Some important parameters are predefined but can be easily changed by setting the attribute in your skript. You can look at all attributes by typing:

vars(Sys)
vars(Exp)
vars(SimOpt)

To get a full overview of all possible attributes please consult the specific section in the documentation.

A first simulation

Define necessary inputs

After importing the classes and simulation modules and setting up the (nearly) empty classes you can fill your skript with live by defining the necessary inputs. The default syntax is:

Class.Attribute = Value

In the following some important attributes are listed. They are either obligatory (o) or it is at least recommended (r) to set a good value. The explanation for the attributes and possible values can be found in the specific section of the documentation or you can consult the quick reference.

  • Sys.spin_system (o)

  • Sys.precursor (o)

  • Sys.population (o)

  • Sys.width_gauss (r)

  • Sys.sigma_time (r)

  • Sys.”spin system parameters” (o)

  • Exp.B_z (r)

  • Exp.t_scale (r)

  • Exp.t_points (r)

  • SimOpt.grid_points (r)

Using these parameters it is possible to run a simulation, nevertheless you should definitly make yourself familiar with further details to get the right things calculated.

Getting the outputs

You can now end your skript by running the simulation. Dependent on the input four different outputs can be expected:

  1. If you set SimOpt.eigval_mode = True the eigenvalues dependent on the magnetic field are given back in a numpy array → eigvals = sim.teacups(Sys, Exp, SimOpt)

  2. If you you run any simulation in hilbert or liouville space a two-dimensional numpy array with the intensities dependent on the magnetic field an the time is returned. → spc = sim.teacups(Sys, Exp, SimOpt)

  3. If you run a simulation in liouville space (SimOpt.space = 'liouville') and set SimOpt.pop_evolution = True the spectrum is returned as well as the evolution of the population of all eigenstates dependent on the time (as a numpy array). → spc, pop_evolution = sim.teacups(Sys, Exp, SimOpt)

Example Skripts

In the following example scripts can be used as a starting point for your own simulations.

Doublet (transient nutations)

import teacups.simulations as sim
import teacups.classes as cl
import numpy as np
import matplotlib.pyplot as plt

plt.style.use("stylesheet.mplstyle")

# initialize classes with default parameters
Sys = cl.Sys()
Exp = cl.Exp()
SimOpt = cl.SimOpt()

# set up spin System parameters
Sys.g = [2, 2, 2]
Sys.width_gauss = 3
Sys.decay = 5e-6


Sys.spin_system = 'doub'
Sys.precursor = 'eigen'
Sys.population = [1, 0]

# set up Experimental parameters
Exp.B_z = np.linspace(320, 380, 3000)
Exp.t_scale = [0, 10e-6]
Exp.t_points = 100
Exp.B_mw = 0.2

# set up simulation SimOption parameters
SimOpt.grid_points = 1
SimOpt.space = 'hilbert'

# do simulation
spec = sim.teacups(Sys, Exp, SimOpt)
spec = spec/abs(spec).real.max()


# %%

plt.figure()
plt.xlabel("$B_z / \mathrm{mT}$")
plt.ylabel("Int./a.u.")
plt.plot(Exp.B_z, spec[25].real)

# plt.savefig("doublet_transient_nutations_B.pdf")


plt.figure()
plt.ylabel("Int./a.u.")
plt.xlabel("$t / \mu\mathrm{s}$")
plt.plot(np.linspace(
    Exp.t_scale[0], Exp.t_scale[1], Exp.t_points), spec[:, 1415].real)

# plt.savefig("doublet_transient_nutations_t.pdf")

spectrum (time)

spectrum (magnetic field)

Doublet (with hyperfines)

import teacups.classes as cl
import matplotlib.pyplot as plt
import numpy as np
import teacups.simulations as sim

plt.style.use("stylesheet.mplstyle")

# initialize classes with default parameters
sys = cl.Sys()
exp = cl.Exp()
opt = cl.SimOpt()

# set up spin system parameters
sys.g = [2, 2, 2]
sys.width_gauss = 3

sys.A1 = [150, 150, 300]
sys.I1 = 1
sys.n1 = 1
sys.A1_frame = [0, 1, 0]

# alternative hyperfine input (2x I=1/2)
# sys.A = [[[250, 250, 250], [250, 250, 250]]]
# sys.I = [[1/2, 1/2]]

sys.spin_system = 'doub'
sys.precursor = 'eigen'
sys.population = [1, 0]

# set up experimental parameters
exp.B_z = np.linspace(320, 380, 600)
exp.t_scale = [0, 2e-6]
exp.t_points = 2

# set up simulation option parameters
opt.grid_points = 20

# do simulation
spec = sim.teacups(sys, exp, opt)
spec = spec[1].real/max(abs(spec[1].real))

# %%
# plot result
plt.figure()
plt.plot(exp.B_z, spec)
plt.xlabel("$B_z$ / mT")
plt.ylabel("norm. intensity")

# plt.savefig("doublet_with_hyperfines.pdf")

spectrum

Triplet (2D)

import teacups.classes as cl
import teacups.simulations as sim
import numpy as np
import matplotlib.pyplot as plt

plt.style.use("stylesheet.mplstyle")

Sys = cl.Sys()
Exp = cl.Exp()
SimOpt = cl.SimOpt()

# choose a triplet spin system
Sys.spin_system = "trip"

# set up g-tensors of the triplet
Sys.g_tri = [2.003, 2.003, 2.003]

# define triplet ZFS
Sys.D_tri = 700
Sys.E_tri = -100

# define initial state
Sys.precursor = "zf"
Sys.population = [1, 0, 0]

# define decay time and line width
Sys.decay = 1e-6
Sys.width_gauss = 1

# experimental setup
Exp.B_z = np.linspace(300, 400, 1024)
Exp.t_scale = [0, 2e-6]
Exp.t_points = 2
Exp.B_mw = 0.001
Exp.freq_mw = 9.75e9

# simulation options
SimOpt.grid_points = 40

# do the simulation
spec_sim = sim.teacups(Sys, Exp, SimOpt)
spec_sim = spec_sim.real[1]/max(abs(spec_sim.real[1]))



plt.figure()
plt.plot(Exp.B_z, spec_sim)

plt.xlabel("$B_z$ / mT")
plt.ylabel("norm. intensity")

# plt.savefig("triplet_2D_hilbert.pdf")

spectrum

Triplet (asymmetric relaxation)

import numpy as np
import matplotlib.pyplot as plt
import teacups.simulations as sim
import teacups.classes as cl

plt.style.use("stylesheet.mplstyle")


# %%
# initialize classes with default parameters
Sys = cl.Sys()
Exp = cl.Exp()
SimOpt = cl.SimOpt()

# set up spin System parameters
Sys.g_tri = [2.008, 2.008, 2.008]
Sys.D_tri = 898
Sys.E_tri = -161

Sys.spin_system = 'trip'
Sys.precursor = 'zf'
Sys.population = [0.95, 0, 0.05]

Sys.decay = 1e-6
Sys.dynamics = np.array([[0, 0.25e6, 0],
                         [0.25e6, 0, 0.01e6],
                         [0, 0.01e6, 0]])

Sys.width_gauss = 2

# set up Experimental parameters
Exp.B_z = np.linspace(295, 395, 256*4)
Exp.t_scale = [1.9e-6, 3.8e-6]
Exp.t_points = 476
Exp.B_mw = 0.0001

# set up simulation SimOption parameters
SimOpt.grid_points = 40
SimOpt.grid = 'sophe'
SimOpt.sym = "D2h"
SimOpt.space = 'liouville'
SimOpt.pop_evolution = True
SimOpt.cpu_cores = 0



# do simulation
spec, pop = sim.teacups(Sys, Exp, SimOpt)

#%%
sim = spec.real/max(abs(spec[25].real))
t = np.linspace(1.9e-6, 3.8e-6, 476)

plt.figure()
plt.plot(t, pop.real)
plt.xlabel("$t$")
plt.ylabel("population of states")

# plt.savefig("znp_triplet_asymmetric_pop.pdf")


#%%
plt.figure()
for i in range(0, 451, 150):
    plt.plot(Exp.B_z, sim[25+i], label=str(t[25+i]))

plt.legend()
plt.xlabel("$B_z$ / mT")
plt.ylabel("norm. intensity")
plt.yticks([])

# plt.savefig("znp_triplet_asymmetric_spec.pdf")

spectrum population evolution

Radical pair (quantum beats)

import teacups.simulations as sim
import teacups.classes as cl
import numpy as np
import matplotlib.pyplot as plt

plt.style.use("stylesheet.mplstyle")

# initialize classes with default parameters
Sys = cl.Sys()
Exp = cl.Exp()
SimOpt = cl.SimOpt()


Sys.spin_system = 'rp'
Sys.g1 = [2.00304, 2.00262, 2.00232]
Sys.g1_frame = [-0.262, 0.489, 0.471]
Sys.g2 = [2.00564, 2.00494, 2.00217]

# Sys.precursor = 'singlet'
# Triplet precursor
Sys.precursor = 'triplet-zf'
Sys.g_tri = [2.00370, 2.00285, 2.00246]
Sys.population = [1, 0, 0]
Sys.D_tri = -1.9217e+03
Sys.E_tri = 525.4678

Sys.decay = 1e-6
Sys.T_relax_1 = 1e-6
Sys.T_relax_2 = 1e-6

Sys.D = 1/3*3.3630
Sys.D_frame = [0, 1.012, -0.017]
Sys.J_ex = 0

Sys.width_gauss = 0.08

Exp.B_z = np.linspace(350.5, 352.3, 800)
Exp.t_scale = [0, 2e-6]
Exp.t_points = 300
Exp.B_mw = 0.03
Exp.freq_mw = 9.8562*1e9

SimOpt.grid_points = 20
SimOpt.space = 'hilbert'
SimOpt.cpu_cores = 0

# do simulation
spec = sim.teacups(Sys, Exp, SimOpt)
spec = spec/abs(spec).real.max()

# %%
plt.figure()
plt.xlabel("$B_z / \mathrm{mT}$")
plt.ylabel("Int./a.u.")
plt.plot(Exp.B_z, spec[150].real)
plt.savefig("rp_quantum_beats_B.pdf")

plt.figure()
plt.ylabel("norm. intensity")
plt.xlabel("$t$ / $\mu$s")
plt.plot(np.linspace(
    Exp.t_scale[0], Exp.t_scale[1], Exp.t_points)*1e6, spec[:, 269].real)
plt.savefig("rp_quantum_beats_t.pdf")

# %% Reproduction of a worse time resolution

Sys.sigma_time = 0.1e-6
SimOpt.extend_t = True

# do simulation
spec = sim.teacups(Sys, Exp, SimOpt)
spec = spec/abs(spec).real.max()

plt.figure()
plt.ylabel("norm. intensity")
plt.xlabel("$t$ / $\mu$s")
plt.plot(Exp.t*1e6, spec[:, 269].real)
plt.savefig("rp_no_quantum_beats_t.pdf")

spectrum (time)

spectrum (magnetic field)

spectrum (bad time resolution)

Radical pair (early dynamics)

import numpy as np
import matplotlib.pyplot as plt
import teacups.simulations as sim
import teacups.classes as cl

plt.style.use("stylesheet.mplstyle")

# %%
Sys = cl.Sys()
Exp = cl.Exp()
SimOpt = cl.SimOpt()

Sys.spin_system = 'rp'
Sys.g1 = [2.00304, 2.00262, 2.00232]
Sys.g2 = [2.00564, 2.00494, 2.00217]
Sys.g1_frame = np.array([-15, 28, 27])*(np.pi/180)

Sys.precursor = 'singlet'

Sys.decay = 0.75e-7
Sys.T_relax_1 = 2e-6
Sys.T_relax_2 = 500e-9

Sys.D = -3.3630
Sys.D_frame = np.array([0, 58, -1])*(np.pi/180)
Sys.J_ex = 0

Sys.width_gauss = 0.125

Exp.B_z = np.linspace(350.55, 352.33, 500)
Exp.t_scale = [0, 100e-9]
Exp.t_points = 11
Exp.B_mw = 0.03
Exp.freq_mw = 9.8562*1e9

SimOpt.grid_points = 10
SimOpt.space = 'hilbert'

# do the simulation
spec = sim.teacups(Sys, Exp, SimOpt)
spec = spec.real/(max(abs(spec[10].real)))



# %%
m = 1
t = [0, 40, 60, 80, 100]

fig, ax1 = plt.subplots(1, 1, figsize=(4, 6))
for m in range(1, 5):
    n = (m+1)*2
    ax1.plot(Exp.B_z, (spec[n].real)-(1.5*m-1), label=t[m])


plt.legend(loc="upper left", ncol=1)
plt.xlabel("$B_z$ / mT")
plt.ylabel("norm. intensity")
plt.yticks([])

# plt.savefig("psi_rp_early_dynamics.pdf")

spectrum

Triplet-doublet pair (2D)

import teacups.classes as cl
import teacups.simulations as sim
import numpy as np
import matplotlib.pyplot as plt

plt.style.use("stylesheet.mplstyle")

Sys = cl.Sys()
Exp = cl.Exp()
SimOpt = cl.SimOpt()

# choose a triplet doublet pair spin system
Sys.spin_system = "tdp"

# set up g-tensors of the radical and the triplet
Sys.g = [2.0059, 2.0059, 2.0059]
Sys.g_tri = [2.003, 2.003, 2.003]

# define triplet ZFS
Sys.D_tri = -500
Sys.E_tri = 50

# define couplings of the radical and the triplet
Sys.J_ex = -400000

Sys.precursor = "triplet-zf"
Sys.population = [0.2, 0.205, 0.2, 0.1, 0.1]

# define decay time and line width
Sys.decay = 1e-6
Sys.width_gauss = 1

# experimental setup
Exp.B_z = np.linspace(333, 363, 550)
Exp.t_scale = [0, 2e-6]
Exp.t_points = 2
Exp.freq_mw = 9.75e9

# simulation options
SimOpt.grid = "sophe"
SimOpt.sym = "D2h"
SimOpt.grid_points = 25
SimOpt.CUPY = False


# do the simulation
spec_sim = sim.teacups(Sys, Exp, SimOpt)
spec_sim = spec_sim.real[1]/max(abs(spec_sim.real[1]))

plt.figure()
plt.plot(Exp.B_z, spec_sim)
plt.xlabel("$B_z$ / mT")
plt.ylabel("norm. intensity")

# plt.savefig("tdp_2D_hilbert.pdf")

spectrum

Triplet-doublet pair (reverse quartet mechanism)

import teacups.simulations as sim
import teacups.classes as cl
import numpy as np
import matplotlib.pyplot as plt

plt.style.use("stylesheet.mplstyle")

Sys = cl.Sys()
Exp = cl.Exp()
SimOpt = cl.SimOpt()

# choose a triplet doublet pair spin system
Sys.spin_system = "tdp"

# set up g-tensors of the radical and the triplet
Sys.g = [2.0059, 2.0059, 2.0059]
Sys.g_tri = [2.003, 2.003, 2.003]

# define triplet ZFS
Sys.D_tri = 700

# define couplings of the radical and the triplet
Sys.J_ex = -20000

# define initial state
Sys.precursor = "eigen"
Sys.population = [0.3, 0.225, 0.2, 0.3, 0.5, 0.5]

# define line width
Sys.width_gauss = 1

# experimental setup
Exp.B_z = np.linspace(320, 380, 700)
Exp.t_scale = [0, 2e-6]
Exp.t_points = 60
Exp.B_mw = 0.001
Exp.freq_mw = 9.75e9

# simulation options
SimOpt.grid_points = 7
SimOpt.space = "liouville"
SimOpt.pop_evolution = True
SimOpt.cpu_cores = 0

# %% set up dynamics-matrix

# determine eigenvalues
SimOpt.eigval_mode = True
eigval = sim.teacups(Sys, Exp, SimOpt)
e = np.mean(eigval[350], axis=0)

# build delta E between trip-doublet and trip-quartet states
de_53 = e[5]-e[3]
de_51 = e[5]-e[1]
de_50 = e[5]-e[0]
de_43 = e[4]-e[3]
de_42 = e[4]-e[2]
de_40 = e[4]-e[0]

# Dipolar coupling squared
D = (Sys.D_tri*1e6)**2

# set isc and doublet decay rates
k_isc = 0.3/1e-11
k_d = 0.25/1e-6
k_e = 0

# set up dynamics matrix
R = np.zeros((6, 6))
R[5, 5] = -(k_d + k_e)
R[4, 4] = -(k_d + k_e)
R[3, 3] = -k_e
R[2, 2] = -k_e
R[1, 1] = -k_e
R[0, 0] = -k_e

R[5, 3] = k_isc/45*(D/(de_53)**2)
R[5, 1] = k_isc/135*(D/(de_51)**2)
R[5, 0] = k_isc/45*(D/(de_50)**2)
R[3, 5] = k_isc/45*(D/(de_53)**2)
R[1, 5] = k_isc/135*(D/(de_51)**2)
R[0, 5] = k_isc/45*(D/(de_50)**2)

R[4, 3] = k_isc/45*(D/(de_42)**2)
R[4, 2] = k_isc/135*(D/(de_42)**2)
R[4, 0] = k_isc/45*(D/(de_40)**2)
R[3, 4] = k_isc/45*(D/(de_43)**2)
R[2, 4] = k_isc/135*(D/(de_42)**2)
R[0, 4] = k_isc/45*(D/(de_40)**2)

Sys.dynamics = R

# %%
SimOpt.eigval_mode = False

# do simulation
spec, pop_evolution = sim.teacups(Sys, Exp, SimOpt)
spec = spec/abs(spec).real.max()

# %%
plt.figure()
plt.xlabel("$B_z$ / mT")
plt.ylabel("norm. intensity")
plt.plot(Exp.B_z, spec[10].real, label="$t=0.3$ $\mu$s")
plt.plot(Exp.B_z, spec[30].real, label="$t=1.0$ $\mu$s")
plt.legend()
# plt.savefig("tdp_rqm_B.pdf")


plt.figure()
plt.ylabel("population")
plt.xlabel("$t$ / $\mu$s")
a = [0, 1, 2, 3, 4, 5]
labels = ["Q$_{-3/2}$", "Q$_{-1/2}$", "Q$_{+1/2}$",
          "Q$_{+3/2}$", "D$_{-1/2}$", "D$_{+1/2}$"]

time = np.linspace(Exp.t_scale[0], Exp.t_scale[1], Exp.t_points)*1e6
for val in a:
    plt.plot(time, pop_evolution[:, val].real /
             pop_evolution.real.max()*1/6, label=labels[val])
plt.legend(ncol=2)
# plt.savefig("tdp_rqm_pop.pdf")

plt.figure()
plt.plot(time, spec[:, 320].real, label="$B_z$ = 347 mT")
plt.legend()
plt.xlabel("$t$ / $\mu$s")
plt.ylabel("norm. intensity")
# plt.savefig("tdp_rqm_t.pdf")

spectrum (time)

spectrum (magnetic field)

population evolution

Workflow using Teacups

When running a simulation in order to calculate a timeresolved spectrum you should first think about what you excactly want to simulate. Think about possible relaxation models and try to determine all parameters of the spin system before starting the simulation as the calculation may be very expensive. In the following a potential workflow is described and illustrated with code snippets how you could run a TEACUPS simulation from the start.

Determine the spin system

First you have to determine the spin system for your transient simulation. Therefore a static spectrum of the spin system can be calculated. You may use some other programms like e.g. EasySpin to simulate any spin system more quickly. To simulate a (nearly) static spectrum in TEACUPS choose as SimOpt.space = 'hilbert'. Further you can calculate only two time points by setting Exp.t_points = 2. If you look at the second time trace you can see your initial static spectrum. Choose carefully the precursor state by setting Sys.precursor:

zf:

only for triplet, population of the zero-field levels X, Y and Z

eigen:

population of the eigenstates of the system

singlet:

only for a radical pair, only the singlet state gets populated

triplet-zf:

for radical pairs or triplet-doublet pairs, population of the triplet precursors zf-levels

As an example in the following file the calculation of a static spectrum of a triplet doublet pair with the initial population of the eigenvalues is shown:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
import teacups.simulations as sim
import teacups.classes as cl
import numpy as np
import matplotlib.pyplot as plt

Sys = cl.Sys()
Exp = cl.Exp()
SimOpt = cl.SimOpt()

# choose a triplet doublet pair spin system
Sys.spin_system = "tdp"

# set up g-tensors of the radical and the triplet
Sys.g = [2.0059, 2.0059, 2.0059]
Sys.g_tri = [2.003, 2.003, 2.003]

# define triplet ZFS
Sys.D_tri = -700

# define couplings of the radical and the triplet
Sys.J_ex = -20000

# define initial state
Sys.precursor = "eigen"
Sys.population = [0.3, 0.23, 0.2, 0.3, 0.5, 0.49]

# define decay time and line width
Sys.decay = 1e-6
Sys.width_gauss = 1

# experimental setup
Exp.B_z = np.linspace(320, 380, 700)
Exp.t_scale = [0, 2e-6]
Exp.t_points = 2
Exp.B_mw = 0.001
Exp.freq_mw = 9.75e9

# simulation options
SimOpt.grid_points = 7
SimOpt.space = "hilbert"

# do the simulation
spec = sim.teacups(Sys, Exp, SimOpt)
plt.plot(Exp.B_z, spec.real[1])
_images/tutorial_teacups_1.png

Analyze the eigenvalues

If you are happy with the initial spin system you should have a look at the eigenvalues of the system. Therefore you have to set SimOpt.eigval_mode = True. If you add the following lines to the code above you’ll get the eigenvalues in ascending order dependent on the magnetic field.

1
2
3
4
SimOpt.eigval_mode = True
eigvals = sim.teacups(Sys, Exp, SimOpt)
plt.figure()
plt.plot(Exp.B_z, eigvals[:, 0])
_images/tutorial_teacups_2.png

Now it is possible to analyze the eigenstates and think about a suitable relaxation mechanism.

Calculate the time evolution

To calculate a certain time evolution you can define Sys.dynamics as a relaxation matrix containing rate constants in 1/s. For example, if you want the eigenstate 5 to decay with the rate constant k_d, set Sys.dynamics[5, 5] = -k_d. For transitions between the eigenstates you can use the off-diagonal elements. E.g. the transition between eigenstates 5 and 2 is described by the elements [5, 2] and [2, 5]. After setting up the relaxation operator, you have to change the simulation options. Set SimOpt.eigval_mode = False in order to calculate a spectrum instead of eigenvalues. Further SimOpt.space = "liouville" activates the calculation of the time evolution in liouville space. If you set SimOpt.pop_evolution = True the matrix Exp.pop_evolution will be calculated additionally, which contains the population of each eigenstate dependent on the time. This is a good poissbility to check your relaxation model and to compare the dependence of the spectrum on the populations.

As an example you can add the following lines to your code to a time resolved spectrum:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
# Define relaxation operator for RQM, Rateconstants in 1/s
k_d = 0.25 / 1e-6

R = np.zeros((6, 6))
R[5, 5] = -k_d
R[4, 4] = -k_d
R[5, 2] = k_d
R[2, 5] = k_d

Sys.dynamics = R

# Change simulationoptions
SimOpt.eigval_mode = False
SimOpt.space = "liouville"
SimOpt.pop_evolution = True
Exp.t_points = 60

# Run the simulation
spec, pop_evolution = sim.teacups(Sys, Exp, SimOpt)

t = np.linspace(Exp.t_scale[0], Exp.t_scale[1], Exp.t_points)*1e6
# %% Plot population evolution
plt.figure()
plt.xlabel("$t / \mu\mathrm{s}$")
plt.ylabel("population of states")
plt.plot(t, np.array(pop_evolution).real)

# %% Plot single 2D
plt.figure()
plt.xlabel("$B_z / \mathrm{mT}$")
plt.plot(Exp.B_z, spec[30].real)
plt.yticks([])

# %% Plot multiple 2D
plt.figure()
plt.xlabel("$B_z / \mathrm{mT}$")
for i in range(14, 30, 5):
    plt.plot(Exp.B_z, spec[i].real, label=str(np.round(t[i], 2))+" $\mu$s")

plt.yticks([])
plt.legend()

# %% Plot surface
plt.figure()
ax = plt.axes()
plt.xlabel("$B_z / \mathrm{mT}$")
plt.ylabel("$t / \mu\mathrm{s}$")
x, y = np.meshgrid(Exp.B_z, t)
ax.pcolormesh(x, y, spec.real, cmap="RdBu", shading="auto")
_images/tutorial_teacups_3.png _images/tutorial_teacups_4.png _images/tutorial_teacups_5.png