spidyan

Simulation of spin dynamics with arbitrary wave forms

Syntax
spidyan(Sys,Exp)
spidyan(Sys,Exp,Opt)
sig = spidyan(...)
[t,sig] = spidyan(...)
[t,sig,out] = spidyan(...)

See also the examples on how to use spidyan and the userguide.

Description

This function simulates spin dynamics of pulse EPR experiments with arbitrary waveforms. The focus of spidyan lies on investigation of spin dynamics and on effects of non-ideal or frequency-swept pulses on the outcome of experiments. Hence, the outputs of spidyan are a time domain signal and density matrices of a single spin packet. For simulating transients or spectra of realistic samples, use the function saffron, which is based one the same propagation engine. Time-domain signals or state trajectories simulated by spidyan are typically further analyzed by (Matlab) software home-written by the user. At the expense of requiring more knowledge on spin quantum mechanics than higher-level EasySpin functions, spidyan provides a lot of freedom in specifying input and output (e.g., custom (non-physical) initial states, transition selective excitation operators, freely-choosable detection operators,...). This freedom should help you to investigate your spin dynamics problem in depth.

Outputs

There are up to three possible output arguments.

A detailed discussion of the structure of the output can be found further below.

Input: Spin system

The three input arguments to spidyan are

Sys is a spin system structure. Most of the regular fields of Sys can be used for the construction of the spin Hamiltonian. Line broadening parameters used by other simulation functions (lw, lwpp, gStrain, etc.) are ignored. Some additional fields can be used with spidyan: To facilitate working with frequency-swept pulses, the field ZeemanFreq is available. For simulations with relaxation, at least generic relaxation times T1 or T2 should be provided. If only one is provided, the other one is set by default to a very large value (1010 microseconds), essentially switching that type of relaxation off.

ZeemanFreq
Array with Zeeman frequencies (in GHz) of electron spins defined in S of the spin system structure. It is possible to use ZeemanFreq in combination with g. If a non-zero Zeeman frequency is found, any g-values given for this electron spin are ignored.
T1
Longitudinal relaxation time in microseconds. If it is a scalar, the relaxation time is applied to all transitions. Transition specific relaxation times can be provided in the form of a matrix. See section Relaxation times below for more details. Relaxation times for transitions that are not defined (or 0) are automatically set to 1010 microseconds. In order for spidyan to simulate relaxation effects, it is necessary to switch relaxation on with Opt.Relaxation (see below).
T2
Transverse relaxation time in microseconds. If it is a scalar, the relaxation time is applied to all transitions. Transition specific relaxation times can be provided in the form of a matrix. See section Relaxation times more details. Relaxation times for transitions that are not defined (or 0) are automatically set to 1010 microseconds. In order for spidyan to simulate relaxation effects, it is necessary to switch relaxation on with Opt.Relaxation (see below).
initState
It is possible to run spidyan simulations with a custom initial state, which must be provided in matrix form. By default, the initial state assumes the high temperature approximation for all electron spins. The structure of density matrices (ordering of operators in terms of product Zeeman basis states) is explained at sop.
Sys.S = 1/2;				% Creates an isolated spin 1/2
Sys.initState = [0 1/2; 1/2 0]; 	% Start simulation from Sx
eqState
For simulations with relaxation it is possible to provide an equilibrium state eqState that the system relaxes to. This state has to be a density operator matrix form. By default the initial state is used as equilibrium state. The equilibrium state is only used if relaxation is active, else it is ignored.
Input: Experimental parameters

Exp contains experimental parameters.

Field
Magnetic field (in mT) at which the experiment is performed, needs to be provided only if the spin system is defined with Sys.g. If Sys.ZeemanFreq is used, Field can be omitted.
Sequence
A cell array that contains the pulse definition and the inter-pulse delays. Available pulse types and their parameters are described at pulse. For a pulse it is necessary to insert a structure containing the pulse definition, free evolution events are declared with a scalar that corresponds to their length in microsecond:
P90.Type = 'rectangular';  	% Define Pulse, see examples for more details
P90.tp = 0.032;		  	% Pulse Length in microseconds
P90.Flip = pi/2;		  	% Flip angle of the pulse

P180.Type = 'rectangular';  	% Define Pulse, see examples for more details
P180.tp = 0.032;		% Pulse Length in microseconds
P180.Flip = pi;		  	% Flip angle of the pulse

Exp.Sequence = {P90 0.2 P180 0.4}; % Creates a two pulse (echo) sequence, with an inter-pulse 
				   % delay of 0.2 microseconds and a final free-evolution period 
				   % of 0.4 microseconds
The delays are defined to go from the end of one previous element in Sequence to the beginning of the next (unlike in Bruker spectrometers).
mwFreq
EPR spectrometer frequency in GHz. All frequencies in the pulse definition (e.g. Pulse.Frequency) need to be defined relative to it:
Exp.mwFreq = 33.5;

Pulse.Frequency = [-250 250];  % 500 MHz sweep
PhaseCycle
Cell array that contains the individual phase cycles, with the indexing corresponding to the pulse indices. The phase cycle itself needs to be given in array form, where rows correspond to the phase cycling steps. In each row the first element is the phase, and the second element the detection phase.
PC = [0, 1; pi, -1];   		% phase cycle
Exp.PhaseCycle = {[] PC};	% phase cycles the second pulse
Phase cycling that steps more than one pulse simultaneously is not available through the PhaseCycle structure, but can be achieved through adding an additional indirect dimension to the experiment (see Dim).
mwPolarization
Optional, can be 'linear' or 'circular', default is 'linear'. If mwPolarization is 'circular', the excitation operator also uses the imaginary part of the waveform. The default excitation then becomes real(IQ) Ŝx + imag(IQ) Ŝy.
Exp.mwPolarization = 'circular';     % mwPolarziation set to circular for all pulses
mwPolarization can also be used in combination with custom excitation operators (Opt.ExcOperator), if they have a imaginary component. Keep in mind, that using 'circular' does not allow spidyan to use some of its speed up tricks and can lead to significantly longer simulation times.

Detection:

DetWindow
The field DetWindow controls detection and is given in microseconds: The zero-time of the detection is taken as the end of the last element in Exp.Sequence.

To do single point detection, it is sufficient to tell DetWindow the time of the detection point. E.g. to detect at the center of the echo after a two pulse experiment:
tau = 0.5
Exp.Sequence = {pulse90 tau pulse180 tau};   	% a two-pulse echo

Exp.DetWindow = 0;	% detect a single point at the end of the second delay (center of the echo)
The detection position can be set to anywhere at the end of the sequence, as long as it does not overlap with a pulse:
Exp.DetWindow =  0.2;	% detect a single point 200 ns after the second delay
Exp.DetWindow = -0.2;	% detect a single point 200 ns before the end of the second delay
In order to detect a transient, DetWindow has to be provided with a start and end value:
tau = 0.5
Exp.Sequence = {pulse90 tau pulse180 tau};   	% a two-pulse echo

Exp.DetWindow = [-0.1 0.1];	% detect a from -100 ns to 100 ns centered around of the end of tau
Of course, it is also possible to detect a transient anywhere after the last pulse, as long is it does not overlap with it:
Exp.DetWindow = [0 0.3];	% start detection at the end of the last event and detect for 300 ns
Exp.DetWindow = [-0.3 -0.1];	% start detection 300 ns before the end of the last event and detect for 200 ns
If you want to detect an FID or have a fixed experiment length (e.g. 3p ESEEM) you can also write:
Exp.Sequence = {pulse90 0.5 pulse90 0.8 pulse90};   	% a three-pulse ESEEM sequence (the last element is a pulse)

Exp.DetWindow = 0.5;	% detect a single point after 500ns after the last pi/2 pulse (center of echo)

Exp.DetWindow = [0.250 0.750] % detect a transient, centered around the echo

Exp.DetWindow = [0 0.100]    % detect an FID after the last pulse
DetSequence

A vector for advanced control of the detection, e.g. detect during pulses. Detection can be set for the entire simulation (meaning detection is on for all events, even pulses) or for specific events:
Exp.DetSequence = true;   	% the entire sequence is detected
Exp.DetSequence = [0 0 1 1];	% expectation values are only returned the third and fourth element in Exp.Sequence
Exp.DetSequence = 0;		% no detection
Besides booleans, Exp.DetSequence also accepts character arrays:
Exp.DetSequence = 'all';   	% the entire sequence is detected
Exp.DetSequence = 'last';	% the last element in Exp.Sequence is detected
If the field DetSequence is not provided, detection is switched off by default. Detection operators can be defined in Exp.DetOperator. If DetSequence is active, but no detection operator is defined, Ŝ+ is used for all electron spins. The length of DetSequence has to be either 1 or the same length as Exp.Sequence.

If Exp.DetWindow is provided, any input from DetSequence is ignored.

DetOperator
The field DetOperator is a cell array that contains detection operators. If no detection operator is defined, but detection is active, Ŝ+ is used for all electron spins. Detection operators can be defined by using the same syntax as in sop (see there):
Exp.DetOperator = {'+1' 'z1'};   % detects Ŝ+ and Ŝz of the first electron spin
Detection operators that can not be defined using the sop syntax, can be provided in matrix form:
Exp.DetOperator = {[0 1; 0 0] [1/2 0; 0 -1/2]};   % the same operators as a above for S = 1/2 in matrix form
For convenience it is possible write a single detection operator without the curly brackets:
Exp.DetOperator = '+1'; % same as {'+1'}

Exp.DetOperator = [1/2 0; 0 -1/2];   % identical to {[1/2 0; 0 -1/2]}

If you are interested in the expectation values of Ŝx, it is usually beneficial to use Ŝ+ or Ŝ- as detection operator and then take the real part of the obtained signal. This removes artifacts at the beginning and end of the time traces that are introduced when translating a purely real signal during the signal processing.

Ŝ+ and Ŝ- are not Hermitian operators. Hence, if you use Ŝ+ as detection operator, the signal that is returned will in fact correspond to ⟨Ŝ-⟩. The real part of the signal (Ŝx) is nof affected by this. But if your data processing involves the imaginary part of the signal and you encounter frequencies with a wrong sign, you might want to consider using Ŝ- instead.

DetFreq
The signal is detected in the simulation frame, which causes all signals from detection operators that have non-zero off-diagonal elements to have an oscillating component. By translating/shifting the frequency it is possible to center the signal around 0 (baseband) or any other desired frequency. The frequency shift has to given as absolute frequency (even if used in conjunction with FrameShift) in GHz and for each detection operator separately. Indexing in DetFreq corresponds to the ordering in DetOperator:
Exp.DetOperator = {'+1' 'z1'};   % detects  Ŝ+ and  Ŝz of the first electron spin
Exp.DetFreq = [33.5 0]; % shifts Ŝ+ by -33.5 GHz
				 % Ŝz does not contain a rotating component and does not need to shifted
For counter-rotating detection operators (e.g. Ŝ-) the sign of the corresponding element in DetFreq has to reflect this:
Exp.DetOperator = {'+1' '-1'};	% detects Ŝ+ and Ŝ- of the first electron spin
Exp.DetFreq = [33.5 -33.5]; 	% shifts Ŝ+ by -33.5 GHz and Ŝ- by 33.5 GHz
DetPhase
Phase of detection in radians, optional. Has to be specified for each detection operator separately:
Exp.DetOperator = {'z1' '+1'};
Exp.DetPhase = [0 pi]; 		% change detection phase of second detection operator by π
The default value is 0. The phase provided in DetPhase is added to the signal, regardless the type of operator.

By defining nPoints and DimX it is possible to vary parameters of the pulse sequence and create one or multidimensional experiments. This also changes the structure of the output compared to a single acquisition. More details on how the output changes can be found below.

nPoints
A vector that contains the number of points in each dimension, e. g. [10 150] corresponds to 10 points in the first and 150 points in the second dimension.
Exp.nPoints = [10 150]; 	% 10 points along the first indirect dimension, 150 along the second
Dim1, Dim2,...
Dim1, Dim2,... provide the fields that are to be changed along the indirect dimensions (field nPoints). The first data point always uses the values defined initially in the experiment definition. All fields that appear in the pulse definition can be changed, e.g:
Exp.Dim1 = {'p2.Flip' pi/8};        % increments the flip angle of the second pulse by pi/8 each step
For free evolution events only the length can be changed:
Exp.Dim1 = {'d3' -0.1};             % decrements the length of the third delay by 100 ns each step
Several parameters can be changed in one dimension:
Exp.Dim1 = {'p2.Flip' pi/8; 'd3' -0.1}; 	% changes flip angle of pulse and duration of free evolution
Exp.Dim1 = {'p2.Flip,p3.Flip' pi/8};  		% flip angles of 2nd and 3rd pulse are simultaneously stepped  
For experiments that involve one or several moving pulses, the identifier Position can be used. This is only possible for pulses that are not the first or last event in the sequence. Pulses are allowed to cross, but must not overlap.
Exp.Dim2 = {'p2.Position' 0.1};   % moves the second pulse 100 ns back each step in the 2nd dimension
A list of increments can be used by providing vectors with the precomputed increments. All increments are always applied to the initial value of the field (as defined in the Exp) and are not related to each other. Hence, if the value in the experiment definition is desired as first data point of the indirect dimension, the fist element has to be zero:
Exp.nPoints = 4;

% this
Exp.Dim1 = {'d2' [0 0.1 0.2 0.3]};
% is equal to: 
Exp.Dim1 = {'d2' 0.1}; 
Complete freedom is given when it comes to providing the list of increments.
Exp.nPoints = 5;
Exp.Dim1 = {'d2' [0 0.1 0.3 0.65 -0.2]};   % increment of the second delay
However, the program checks that changing lengths of events do not lead to overlapping pulses at any acquisition point.

A special case is the field Par.Frequency in the pulse definition of a frequency-swept or of a monochromatic pulse that is defined by identical initial and final frequency. If one value is provided in Exp.DimX, this is going to be added to both elements in the Par.Frequency field - the pulse is moved in the frequency domain:

Exp.Dim1 = {'p1.Frequency' 5};   % changes both elements in Par.Frequency by 5 MHz per step
An equivalent input as the above is the following, where each the increment for each field in Par.Frequency is given:
Exp.Dim1 = {'p1.Frequency' [5 5]};   % changes both elements in Par.Frequency by 5 MHz per step
This syntax can also be used to change the sweep width:
Exp.Dim1 = {'p1.Frequency' [-5 5]};   			% increases sweep width by 10 MHz each step
Exp.Dim1 = {'p1.Frequency' [0 5]};    			% increment only the final frequency of the frequency sweep
A list of increments can be provided as well, by using ';' as separator:
Exp.Dim1 = {'p1.Frequency' [0 0; 10 10; 30 30; 80 80]};  % vector increment for shifting the frequency range of the pulse
Exp.Dim1 = {'p1.Frequency' [0 0; -10 0; 0 30; -80 80]};  % vector increment changing the excitation range of the pulse

For other pulse parameters that are defined by a vector (e.g., the order of an asymmetric hyperbolic secant (HS) pulse Pulse.n or the list of relatives amplitudes of a Gaussian cascade Pulse.A0), selected elements can be incremented by adding an index to the field name:

Exp.Dim1 = {'p1.A0(3)' 0.1; 'p1.A0(4)' -0.1};   % changes relative amplitudes of the third and fourth pulse in a Gaussian cascade
Exp.Dim1 = {'p1.n(2)' 2};    % increases order of the falling flank of a hyperbolic secant pulse by 2 each step

% Also possible for 'Frequency'
Exp.Dim1 = {'p1.Frequency(2)' 5};    % changes only the final frequency of the pulse
Exp.Dim1 = {'p1.Frequency' [0 5]};   % identical to the above and not a list of increments!
Adding an indirect dimension also allows for simultaneous phase cycling of two or more pulses (something that can not be achieved through the Exp.PhaseCycle structure):
Exp.nPoints = 4;
Exp.Dim1 = {'p2.Phase,p3.Phase' pi/4};   % changes the phase of the 2nd and 3rd pulse by pi/4 each step
In the above example, the output of spidyan will contain the individual transients from each phase cycling step and manual merging of the dimensions (with proper detection phase/sign) is required to obtain the phase-cycled signal.

This is an example for a two-dimensional experiment (e.g. for a two-pulse echo) to optimize two pulse parameters:

Exp.nPoints = [20 15];
Exp.Dim1 = {'p1.Frequency,p2.Frequency' [-2.5 2.5]}; % increase excitation band by 5 MHz each step of both pulses
Exp.Dim1 = {'p1.tp' 2; 'p2.tp' 1}; % step pulse lengths, in µs

Bandwidth limiting effects of a resonator on a pulse can be incorporated in two ways:

The resonator profile is defined using the resonator center frequency ResonatorFrequency and the loaded Q-value ResonatorQL. The resonator frequency response is then computed from the ideal transfer function for an RLC series circuit)

Exp.ResonatorFrequency = 33.5; 	% resonator center frequency in GHz
Exp.ResonatorQL = 300;	  	% loaded Q-value

Or the transfer function in combination with the frequency axis are given with FrequencyResponse.

Exp.FrequencyResponse = [Frequency; TransferFunction] 	% transfer function and its correspond frequency axis

By setting ResonatorMode it is possible to compensate for the resonator. Compensation aims to provide a waveform that excites spin packets in the resonator with constant critical adiabaticity if the originally specified waveform would lead to excitation with constant critical adiabaticity in the absence of a resonator (mainly for chirp and hyperbolic secant pulses). For more details see pulse and resonator.

ResonatorFrequency
Center frequency of resonator in GHz.
ResonatorQL
The loaded Q-value.
FrequencyResponse
This gives the frequency response of the resonator, in the form Par.FrequencyResponse = [Frequency; TransferFunction] with the (possibly experimental) resonator transfer function in TransferFunction and the corresponding frequency axis in Frequency (in GHz). A complex transfer function input in Par.FrequencyResponse is used directly in the bandwidth compensation. A real transfer function input is assumed to correspond to the magnitude response, and the associated phase response is estimated.

If FrequencyResponse is given, the alternative input fields ResonatorFrequency and ResonatorQL are ignored.

ResonatorMode
Optional, can be 'simulate' or 'compensate'. By default the effect of the resonator on the signal is simulated. If set to 'compensate' the pulse shape is adapted such that it compensates for the resonator profile by using a uniform critical adiabaticity criterion.
Input: Simulation options

Opt contains additional simulation parameters.

Relaxation
Relaxation can be controlled for the entire sequence or specific events:
Opt.Relaxation = true;   	% switches on relaxation for the entire simulation
Opt.Relaxation = [0 0 1];	% relaxation is active only during the third event
If relaxation is switched on, relaxation times (Sys.T1 and Sys.T2) have to be provided and the simulation is then performed in Liouville space, which slows down the simulation, for larger spin systems tremendously so. In order to avoid unexpected behavior, if the Opt.Relaxation is not used to switch relaxation on/off for the entire sequence, the length of it has to match the length of Exp.Sequence (all events need to be defined).
SimFreq
Optional, frequency of the (rotating) simulation frame in GHz. If SimFreq is not provided, a suitable simulation frame is automatically selected. No matter the choice of the simulation frame, all other frequencies must still be given in the lab frame - the program handles all required frequency shifts. For example, for running a Q-band simulation without using the default simulation frame, but one that rotates at 30 GHz:
Exp.mwFreq = 34.4;	        % Set carrier frequency to 34.4 GHz
Pulse.Frequency = [-50 50]; % Pulse with bandwidth of 100 MHz, sweeps from 34.35 to 34.45 GHz 

Opt.SimFreq = 30;  	      % Changes into a rotating frame at 30 GHz
To force a simulation in the labframe SimFreq has to be set to zero:
Opt.SimFreq = 0;  	      % Propagation in the labframe
Changing the simulation frame allows for using a larger time step, which in turn can strongly reduce computation time. Keep in mind, that if a time step is provided with Opt.IntTimeStep it is not adapted automatically. To have the program adapt the time step, remove the field Opt.IntTimeStep .
IntTimeStep
Integration time step in microseconds, optional. If no time step is provided, the program computes a suitable time step (taking into account the highest frequency that is being used for the pulse definition Exp.Frequency). Accurate results are obtained for an integration time step that corresponds to about 50 data points per oscillation or 1/50th of the largest time step size to fulfil the Nyquist criterion.

User provided integration time steps must at least fulfil the Nyquist criterion, or an error is thrown. For integration time steps larger than the recommended ones, a warning is printed to the command line.

StateTrajectories
If StateTrajectories is active, the output contains a cell array with density matrices at each time step during events that are detected. StateTrajectories can be switched on for the entire sequence as well as for specific events:
Opt.StateTrajectories = true;   	% switches on state trajectories for the entire simulation
Opt.StateTrajectories = [0 0 1];	% state trajectories is active only during the third event
State trajectories are only fully recorded for events where detection is active. Otherwise, the density matrices are stored only at the beginning and end of the event. In order to avoid unexpected behavior, if the Opt.StateTrajectories is not control the entire sequence, the length of it has to match the length of Exp.Sequence (all events need to be defined).
ExcOperator
By default Ŝx is used as excitation operator. With the cell array ExcOperator it is possible to use custom excitation operators. The indexing of ExcOperator corresponds to the pulse index. It is possible to use the syntax from sop as well as matrices:
Opt.ExcOperator = {[] 'x(1|2)'};     % transition selective excitation operator for the second pulse
Opt.ExcOperator = {[0 1/2; 1/2 0]};  % custom excitation operator for the first pulse in matrix form
The imaginary part of a custom excitation operator is only considered if ComplexExcitation is active for the respective pulse.
Output Structure
Time Traces

The data type of the returned time-domain signal depends on experiment parameters. If you run a single experiment (no Exp.nPoints), the output sig is a two-dimensional numeric array, with the first dimension corresponding being the detected transient and the second dimension specified detection operators. If one the specified detection operators is a ladder operator, say Ŝ-, sig is a complex-valued vector. The time axis is returned in t and will be a vector with the same length as sig.

If you run a multidimensional experiment, sig can be an (n+2)-dimensional numeric array (if each acquisition point has the same length) or an n-dimensional cell array (if the detection length changes), where n is the number of indirect dimensions of your simulation. The indexing of n corresponds to how Exp.nPoints was defined:

If the total length of all detected events is identical for each acquisition point, sig will be a (n+2) dimensional array. For Exp.nPoints = [3 4]; a total of 12 data sets will be returned in a four dimensional array:

Exp.nPoints = [3 4];
Exp.DetOperator = {'z1','+1'};
size(sig) = 

           3	  4	  11003		2 	
The first dimension of your output has a length of 3 and the second dimension 4. These correspond to the indirect dimensions. The second-to-last dimension corresponds are the expectation values at each point of time during which was detected (the transient) and the last dimension to the detection operators. To get the time trace that corresponds to the second point in the first dimension and the last point in the second dimension you could use:
trace = sig(2,4,:,:);
% or if you want to remove the singular dimensions (e.g. for plotting)
trace = squeeze(sig(2,4,:,:));
If all acquisition points have the same time axis, t is a vector:
size(t) = 

           1 	11003
If not, t is an (n+1) dimensional numeric array, with the same indexing for n as for sig. This can be encountered when detecting a moving echo, with a detection window that is always centered around the echo, but moves in terms of absolute time (see example on two-pulse ESEEM).
size(t) = 

           3        4 	   11003

When the detection length changes (e.g. if you vary the length of a pulse or a free evolution event that is being detected) between acquisition points, t and sig become n-dimensional cell arrays. The indexing of the elements in each cell array again corresponds to Exp.nPoints:

size(sig) = 

           3	  4
Each element of sig is a two-dimensional numeric array with the first dimension corresponding to the detection operators and the second to the time axis:
size(sig{1,1}) = 

           2 	11003
		   
size(sig{2,4}) = 

           2 	21003
The time axes for each acquisition point are stored as vectors in the elements of t:
size(t) = 

           3	  4
		   
size(t{1,1}) = 

           2 	11003
		   
size(t{2,4}) = 

           2 	21003

Final States

For a single acquisition, the output out.FinalState is the final density matrix of size nH x nH, where nH is the dimensionality of your spin system in Hilbert space. If a simulation with more than one acquisition point is run, s is an (n+2)-dimensional numeric array, with the same rules for indexing as demonstrated above for sig. E. g. to get the final state of the second data point in the first dimension and the fourth data point in the second dimension, write:

InterestingState = squeeze(out.FinalState(2,4,:,:)); % squeeze removes the singleton dimensions
InterestingState = 

		0.0592 - 0.0000i   0.4559 + 0.1965i
		0.4559 - 0.1965i  -0.0592 + 0.0000i

State Trajectories

If Opt.StateTrajectories is active during at least one event, out.StateTrajectories contains the state trajectories. This means, it contains the state (density matrix) of each propagation step during the events selected in Opt.StateTrajectories. out.StateTrajectories is an n-dimensional cell array where the dimensions once again correspond to Exp.nPoints. To get the state trajectories for the second data point in the first dimension and for the fourth datapoint in the second dimension of the example above, write:

InterestingStateTrajectory = out.StateTrajectories{2,4}
size(out.StateTrajectories)

     3		4
size(InterestingStateTrajectory)

     1		21003
For full recording of state trajectories, detection needs to be active during those events. Otherwise only the states at the beginning and the end of that event are recorded.

See also

saffron, resonator, signalprocessing, pulse, sop