function sig = liouville_von_neumann_fun_1(T,B0)
hfc = 1.0; % Hyperfine coupling
B0 = B0./1e4; % Convert B0 to Tesla

% Calculation options
opt.fixtrace = true; % Fix trace of the density matrix to 1 after each timestep
opt.check = true; % Do sanity checks (for advanced users only)
 
%% Basic spin matrices (in alpha/beta basis)
Ix = 0.5 * [0  1; 1   0];
Iy = 0.5 * [0 -1i;1i  0];
Iz = 0.5 * [1  0; 0  -1];
Ie = eye(2);
 
tkron = @(x,y,z) kron(kron(x,y),z); % Successive Kronecker product of three elements
tol = 10*eps; % Tolerance for asserts
 
%% Transform to Singlet Triplet basis
s2 = 1/sqrt(2);
%      T+ S  T0   T-
Pe = [1  0   0   0;
      0 s2  s2   0;
      0 -s2 s2   0;
      0  0   0   1];
P = kron(Pe,eye(2)); % Nucleus stays in alpha/beta basis
ST = @(A) P*A*P'; % Transformation to the ST basis
 
% Check that P is unitary
if opt.check
    condition = ST(eye(8)) - eye(8) < tol;
    assert(all(condition(:)),'P is not unitary');
end
 
%% Product spin matrices
% Note: In Matlab one cannot (easily) work on the vector operators.
% Therefore, we work with the components thereof.
 
% Components of SA vector operator
SAx = ST(tkron(Ix,Ie,Ie));
SAy = ST(tkron(Iy,Ie,Ie));
SAz = ST(tkron(Iz,Ie,Ie));
% Raising and lowering operators
SAp = SAx + 1i*SAy; SAm = SAx - 1i*SAy;
 
% Components of SB vector operator
SBx = ST(tkron(Ie,Ix,Ie));
SBy = ST(tkron(Ie,Iy,Ie));
SBz = ST(tkron(Ie,Iz,Ie));
SBp = SBx + 1i*SBy;
SBm = SBx - 1i*SBy;
 
% Components of SC vector operator
SCx = ST(tkron(Ie,Ie,Ix));
SCy = ST(tkron(Ie,Ie,Iy));
SCz = ST(tkron(Ie,Ie,Iz));
 
% Total spin angular momentum (squared) of each species
SA2 = SAx^2 + SAy^2 + SAz^2;
SB2 = SBx^2 + SBy^2 + SBz^2;
SC2 = SCx^2 + SCy^2 + SCz^2;
 
% Product operators
SASB = SAx*SBx + SAy*SBy + SAz*SBz;
SASC = SAx*SCx + SAy*SCy + SAz*SCz;
SBSC = SBx*SCx + SBy*SCy + SBz*SCz;
 
% Total electron spin operator (squared)
SAB2 = SA2 + SB2 + 2*SASB; % Method #1
%SAB2 = SA2 + SB2 + 2*SAz*SBz + SAp*SBm + SBp*SAm; % Method #2
 
% Components of the SAB vector operator
SABx = SAx + SBx;
SABy = SAy + SBy;
SABz = SAz + SBz;
 
% Product of the vector operator SAB SC
SABSC = SABx*SCx + SABy*SCy + SABz*SCz;
 
% Total spin operator F2 = (SA + SB + SC)^2 <- SA, SB, SC are vector operators
%F2 = SA2 + SB2 + SC2 + 2*SASB + 2*SASC + 2*SBSC; % Method #1
F2 = SAB2 + SC2 + 2*SABSC; % Method 2
 
gn = 267.513e3; % Hydrogen nucleus gyromagnetic ratio in rad/s/mT
ge = 1.760859644e8; % Electron gyromagnetic ratio in rad/s/mT
 
%% Make the Hamiltonian
a = ge * hfc; % Hyperfine coupling constant
omega0 = ge*B0; % Electron Larmor frequency
omega0n = -gn*B0; % Nuclear Larmor frequency
%disp2(a,omega0);

H_mag = omega0.*(SAz + SBz) + omega0n.*SCz; % Electron and nuclear Zeeman
H_hyperfine = a*(SAx*SCx + SAy*SCy + SAz*SCz);
 
H = H_mag + H_hyperfine;
 
%% Define the projection operators, propagators and the observable
% Singlet projection operator in the aba basis
PS = 0.25*eye(8) - tkron(Ix,Ix,Ie) - tkron(Iy,Iy,Ie) - tkron(Iz,Iz,Ie);
PS = ST(PS); % Project into ST basis

%Create the initial density matrix
% Starts as singlet with no nuclear polarisation
rh0 = zeros(8); rh0(5,5) = 0.5; rh0(6,6) = 0.5;
 
% Sanity checks
if opt.check
    assert(abs(trace(rh0)-1) < tol,'Invalid starting point trace')
    condition = diag(rh0) >= -tol;
    assert(all(condition(:)),'Invalid sign of the starting point')
end
    
% Calculate the propagator matrices
dt = mean(diff(T));
Um = expm(-1i*H*dt); Up = expm(1i*H*dt);
rhp = rh0; % Spin density at previous timestep
 
% Define the function to calculate the observable on the density matrix
observable = @(rh) real(trace(PS*rh)); % Singlet character
 
%% Iterate through all the time points
sig = NaN(1,length(T)); % Preallocate the matrix for the signal
sig(1) = observable(rh0); % Calculate the observable of the initial density matrix
Tsig = (0:length(T)-1)*dt; % Timesteps at which the 'sig' is simulated
for i = 2:length(T)
    % Propagate the density matrix by dt
    rht = Um*rhp*Up;
    if opt.fixtrace && opt.check; rht = rht/trace(rht); end % Fix the roundoff errors on trace
    rhp = rht;

    % Sanity checks
    if opt.check
        % The trace of the density matrix must be equal to one (within machine precision)
        % Error accumulates therefore multiply `tol` by `i`
        assert(abs(trace(rht) - 1) < i*tol,'Invalid trace %d', abs(trace(rht) - 1)/tol)
 
        % Diagonal elements of the density matrix must be positive (within machine precision)
        condition = diag(rht) >= -tol;
        assert(all(condition(:)),'Invalid sign')
    end
    
    % Calculate the observable at time t
    sig(i) = observable(rht);
end
 
% Interpolate from the simulation points to the desired time points
sig = interp1(Tsig,sig,T,'pchip'); % Cubic interpolate the desired T points
 
end