Added sequence examples.

examples/COMPFID.py

@@ -0,0 +1,59 @@
+"""Composite FID example.
+
+This example demonstrates how to simulate a composite FID signal using the quackseq simulator.
+
+See the paper: 
+Sauer, K.L., Klug, C.A., Miller, J.B. et al. Using quaternions to design composite pulses for spin-1 NQR. Appl. Magn. Reson. 25, 485–500 (2004). https://doi.org/10.1007/BF03166543
+
+This also works for Samples with spin > 1.
+"""
+
+import logging
+
+from quackseq_simulator.simulator import Simulator
+from quackseq.pulsesequence import QuackSequence
+from quackseq.functions import RectFunction
+from matplotlib import pyplot as plt
+
+if __name__ == "__main__":
+    logging.basicConfig(level=logging.INFO)
+
+    logger = logging.getLogger(__name__)
+
+    seq = QuackSequence("COMPFID")
+    seq.add_pulse_event("tx1", "3u", 100, 0, RectFunction())
+    seq.add_pulse_event("tx2", "6u", 100, 45, RectFunction())
+    # This makes the phase 45, 135, 225, 315
+    seq.set_tx_n_phase_cycles("tx2", 4)
+    seq.add_blank_event("blank", "5u")
+    
+    seq.add_readout_event("rx", "100u")
+
+    # No phase shifiting of the receive data but weighting of -1 for the 45 degree pulse, +1 for the 135 degree pulse, -1 for the 225 degree pulse and +1 for the 315 degree pulse
+    readout_scheme = [[1, 0], [-1, 0], [1, 0], [-1, 0]]
+
+    sim = Simulator()
+    sim.set_averages(100)
+
+    sim.settings.noise = 1 # microvolts
+
+    result = sim.run_sequence(seq)
+    # Plot time and frequency domain next to each other
+    plt.subplot(1, 2, 1)
+    plt.title("Time domain Simulation of BiPh3 COMPFID")
+    plt.xlabel("Time (µs)")
+    plt.ylabel("Signal (a.u.)")
+    plt.plot(result.tdx[-1], result.tdy[-1].imag, label="imaginary")
+    plt.plot(result.tdx[-1], result.tdy[-1].real, label="real")
+    plt.plot(result.tdx[-1], abs(result.tdy[-1]), label="abs")
+
+    plt.subplot(1, 2, 2)
+    plt.title("Frequency domain Simulation of BiPh3 COMPFID")
+    plt.xlabel("Frequency (kHz)")
+    plt.ylabel("Signal (a.u.)")
+    plt.plot(result.fdx[-1], result.fdy[-1].imag, label="imaginary")
+    plt.plot(result.fdx[-1], result.fdy[-1].real, label="real")
+    plt.plot(result.fdx[-1], abs(result.fdy[-1]), label="abs")
+
+    plt.legend()
+    plt.show()

examples/FID.py

@@ -0,0 +1,48 @@
+"""Example Free Induction Decay (FID) simulation using the quackseq simulator.
+
+The sample is the default BiPh3 NQR sample.
+"""
+
+import logging
+
+from quackseq_simulator.simulator import Simulator
+from quackseq.pulsesequence import QuackSequence
+from quackseq.functions import RectFunction
+from matplotlib import pyplot as plt
+
+if __name__ == "__main__":
+    logging.basicConfig(level=logging.INFO)
+
+    logger = logging.getLogger(__name__)
+
+    seq = QuackSequence("FID")
+    seq.add_pulse_event("tx", "3u", 100, 0, RectFunction())
+    seq.add_blank_event("blank", "5u")
+    seq.add_readout_event("rx", "100u")
+    seq.add_blank_event("TR", "1m")
+
+    sim = Simulator()
+    sim.set_averages(100)
+
+    sim.settings.noise = 1 # microvolts
+
+    result = sim.run_sequence(seq)
+    # Plot time and frequency domain next to each other
+    plt.subplot(1, 2, 1)
+    plt.title("Time domain Simulation of BiPh3 FID")
+    plt.xlabel("Time (µs)")
+    plt.ylabel("Signal (a.u.)")
+    plt.plot(result.tdx[0], result.tdy[0].imag, label="imaginary")
+    plt.plot(result.tdx[0], result.tdy[0].real, label="real")
+    plt.plot(result.tdx[0], abs(result.tdy[0]), label="abs")
+
+    plt.subplot(1, 2, 2)
+    plt.title("Frequency domain Simulation of BiPh3 FID")
+    plt.xlabel("Frequency (kHz)")
+    plt.ylabel("Signal (a.u.)")
+    plt.plot(result.fdx[0], result.fdy[0].imag, label="imaginary")
+    plt.plot(result.fdx[0], result.fdy[0].real, label="real")
+    plt.plot(result.fdx[0], abs(result.fdy[0]), label="abs")
+
+    plt.legend()
+    plt.show()

examples/SE.py

@@ -0,0 +1,49 @@
+"""Example Spin Echo (SE) simulation using the quackseq simulator.
+
+The sample is the default BiPh3 NQR sample.
+"""
+
+import logging
+
+from quackseq_simulator.simulator import Simulator
+from quackseq.pulsesequence import QuackSequence
+from quackseq.functions import RectFunction
+from matplotlib import pyplot as plt
+
+if __name__ == "__main__":
+    logging.basicConfig(level=logging.INFO)
+
+    logger = logging.getLogger(__name__)
+
+    seq = QuackSequence("SE")
+    seq.add_pulse_event("pi-half", "3u", 100, 0, RectFunction())
+    seq.add_blank_event("te-half", "150u")
+    seq.add_pulse_event("pi", "6u", 100, 0, RectFunction())
+    seq.add_blank_event("blank", "50u")
+    seq.add_readout_event("rx", "200u")
+
+    sim = Simulator()
+    sim.set_averages(100)
+
+    sim.settings.noise = 1 # microvolts
+
+    result = sim.run_sequence(seq)
+    # Plot time and frequency domain next to each other
+    plt.subplot(1, 2, 1)
+    plt.title("Time domain Simulation of BiPh3 SE")
+    plt.xlabel("Time (µs)")
+    plt.ylabel("Signal (a.u.)")
+    plt.plot(result.tdx[0], result.tdy[0].imag, label="imaginary")
+    plt.plot(result.tdx[0], result.tdy[0].real, label="real")
+    plt.plot(result.tdx[0], abs(result.tdy[0]), label="abs")
+
+    plt.subplot(1, 2, 2)
+    plt.title("Frequency domain Simulation of BiPh3 SE")
+    plt.xlabel("Frequency (kHz)")
+    plt.ylabel("Signal (a.u.)")
+    plt.plot(result.fdx[0], result.fdy[0].imag, label="imaginary")
+    plt.plot(result.fdx[0], result.fdy[0].real, label="real")
+    plt.plot(result.fdx[0], abs(result.fdy[0]), label="abs")
+
+    plt.legend()
+    plt.show()

examples/SEPC.py

@@ -0,0 +1,58 @@
+"""Spin Echo with Phase Cycling (SEPC) simulation using the quackseq simulator.
+
+The sample is the default BiPh3 NQR sample.
+"""
+
+import logging
+
+from quackseq_simulator.simulator import Simulator
+from quackseq.pulsesequence import QuackSequence
+from quackseq.functions import RectFunction
+from matplotlib import pyplot as plt
+
+if __name__ == "__main__":
+    logging.basicConfig(level=logging.INFO)
+
+    logger = logging.getLogger(__name__)
+
+    seq = QuackSequence("SEPC")
+    seq.add_pulse_event("pi-half", "3u", 100, 0, RectFunction())
+    # This causes the phase to cycle through 0, 90, 180, 270
+    seq.set_tx_n_phase_cycles("pi-half", 4)
+
+    seq.add_blank_event("te-half", "150u")
+    # For the second pulse we just need a phase of 180
+    seq.add_pulse_event("pi", "6u", 100, 180, RectFunction())
+    seq.add_blank_event("blank", "50u")
+
+    seq.add_readout_event("rx", "200u")
+    # Readout scheme for phase cycling TX pulses have the scheme 0 90 180 270 for the first pulse and 180 always for the second pulse
+    readout_scheme = [[1, 0], [1, 90], [1, 180], [1, 270]]
+
+    seq.set_rx_readout_scheme("rx", readout_scheme)
+
+    sim = Simulator()
+    sim.set_averages(100)
+
+    sim.settings.noise = 1 # microvolts
+
+    result = sim.run_sequence(seq)
+    # Plot time and frequency domain next to each other
+    plt.subplot(1, 2, 1)
+    plt.title("Time domain Simulation of BiPh3 SEPC")
+    plt.xlabel("Time (µs)")
+    plt.ylabel("Signal (a.u.)")
+    plt.plot(result.tdx[-1], result.tdy[-1].imag, label="imaginary")
+    plt.plot(result.tdx[-1], result.tdy[-1].real, label="real")
+    plt.plot(result.tdx[-1], abs(result.tdy[-1]), label="abs")
+
+    plt.subplot(1, 2, 2)
+    plt.title("Frequency domain Simulation of BiPh3 SEPC")
+    plt.xlabel("Frequency (kHz)")
+    plt.ylabel("Signal (a.u.)")
+    plt.plot(result.fdx[-1], result.fdy[-1].imag, label="imaginary")
+    plt.plot(result.fdx[-1], result.fdy[-1].real, label="real")
+    plt.plot(result.fdx[-1], abs(result.fdy[-1]), label="abs")
+
+    plt.legend()
+    plt.show()

src/quackseq_simulator/simulator_model.py

@@ -38,7 +38,7 @@ class SimulatorModel(SpectrometerModel):
 
     # Sample settings, this will  be done in a separate module later on
     SAMPLE_NAME = "Name"
-    NUMBER_ATOMS = "Number of atoms per unit volume (1/m^3)"
+    NUMBER_ATOMS = "N. atoms (1/m^3)"
     DENSITY = "Density (g/cm^3)"
     MOLAR_MASS = "Molar mass (g/mol)"
     RESONANT_FREQUENCY = "Resonant freq. (MHz)"

tests/simulator.py

@@ -170,10 +170,22 @@ class TestQuackSequence(unittest.TestCase):
         # Plotting the results
         plt.title("Phase cycling")
         logger.info(f"Number of data sets {len(result.tdy)}")
-        plt.plot(result.tdx[0], result.tdy[0].real, label="pc 1")
+        plt.plot(result.tdx[0], result.tdy[0].real, label="pc 1 real")
+        plt.plot(result.tdx[0], result.tdy[0].imag, label="pc 2 imag")
+        plt.legend()
+        plt.show()
         plt.plot(result.tdx[1], result.tdy[1].real, label="pc 2")
+        plt.plot(result.tdx[1], result.tdy[1].imag, label="pc 2")
+        plt.legend()
+        plt.show()
         plt.plot(result.tdx[2], result.tdy[2].real, label="pc 3")
+        plt.plot(result.tdx[2], result.tdy[2].imag, label="pc 3")
+        plt.legend()
+        plt.show()
         plt.plot(result.tdx[3], result.tdy[3].real, label="pc 4")
+        plt.plot(result.tdx[3], result.tdy[3].imag, label="pc 4")
+        plt.legend()
+        plt.show()
         plt.plot(result.tdx[4], abs(result.tdy[4]), label="Phase Cycling")
         plt.legend()
         plt.show()