function sig = liouville_von_neumann_fun_2(T,hfc,J,B0)
% Calculation options
opt.fixtrace = false; % Fix trace of rho to 1 after each timestep
opt.check = false; % Do sanity checks (for advanced users only)

%% Setup
% 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); % Identity matrix

% Successive Kronecker product of three and four elements
tkron = @(x,y,z) kron(kron(x,y),z);
qk = @(x,y,z,q) kron(kron(kron(x,y),z),q);

tol = 10*eps; % Tolerance for asserts
 
% Similarity transform to Singlet-Triplet basis
s2 = 1/sqrt(2); Pe = [1 0 0 0; 0 s2 s2 0; 0 -s2 s2 0; 0 0 0 1];
P = tkron(Pe,Ie,Ie); % Nuclei stay in alpha/beta basis
ST = @(A) P*A*P'; % Similarity transform to the ST basis
 
% Check that P is unitary
if opt.check
    condition = ST(eye(16)) - eye(16) < tol;
    assert(all(condition(:)),'P is not unitary');
end
 
%% Spin matrices
% Components of SA,B,C,D vector operators
SAx=ST(qk(Ix,Ie,Ie,Ie)); SAy=ST(qk(Iy,Ie,Ie,Ie)); SAz=ST(qk(Iz,Ie,Ie,Ie));
SBx=ST(qk(Ie,Ix,Ie,Ie)); SBy=ST(qk(Ie,Iy,Ie,Ie)); SBz=ST(qk(Ie,Iz,Ie,Ie));
SCx=ST(qk(Ie,Ie,Ix,Ie)); SCy=ST(qk(Ie,Ie,Iy,Ie)); SCz=ST(qk(Ie,Ie,Iz,Ie));
SDx=ST(qk(Ie,Ie,Ie,Ix)); SDy=ST(qk(Ie,Ie,Ie,Iy)); SDz=ST(qk(Ie,Ie,Ie,Iz));
SASC = SAx*SCx + SAy*SCy + SAz*SCz; % A-C product operator
SASD = SAx*SDx + SAy*SDy + SAz*SDz; % A-D product operator
SASB = SAx*SBx + SAy*SBy + SAz*SBz; % A-C product operator

%% Make the Hamiltonian
gamma_n = 42.577e2; % Hydrogen nucleus gyromagnetic ratio in rad/s/G
gamma_e_1 = 28024.951e2; % Electron gyromagnetic ratio in MHz/G
gamma_e_2 = (2.01/2.00)*28024.951e2; % Electron gyromagnetic ratio in MHz/G
J = J.*1e6; % convert to MHz
% Hyperfine frequencies of spin A
a1 = gamma_e_1 * hfc(1); a2 = gamma_e_1 * hfc(2);
omega_1 = gamma_e_1*B0; % Electron 1 Larmor frequency
omega_2 = gamma_e_2*B0; % Electron 2 Larmor frequency
omega_n = gamma_n*B0; % Nuclear Larmor frequency
fprintf('B0 = %.1f G\n',B0)
fprintf('a1 = %.1f MHz | a2 = %.1f MHz \n',a1/1e6,a2/1e6);
fprintf('omega1 = %.1f MHz | omega2 = %.1f MHz\n',omega_1/1e6,omega_2/1e6);
fprintf('omega_n = %.1f MHz\n',omega_n/1e6);
fprintf('J = %.1f MHz\n',J/1e6);
% Electron and nuclear Zeeman
H_mag = omega_1.*SAz + omega_2.*SBz + omega_n.*SCz + omega_n.*SDz;
H_hyperfine = a1.*SASC + a2.*SASD;
H_J = J.*SASB;
H = H_mag + H_hyperfine + H_J;
 
%% Define the projection operators, propagators and the observable
% Singlet projection operator in the aba basis
PS = 0.25*eye(16)-qk(Ix,Ix,Ie,Ie)-qk(Iy,Iy,Ie,Ie)-qk(Iz,Iz,Ie,Ie);
PS = ST(PS); % Project into ST basis
rh0 = PS./sum(trace(PS)); % Starts as singlet with no nuclear polarisation
PS = sparse(PS);

% 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));
dt = min(1./[a1 a2 omega_1 omega_n])/2;
dt = min(dt,1e-6);
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
observable = @(rh) real(sum(diag(PS*rh))); % Singlet character

%% Iterate through all the time points
Tsig = (0:ceil(max(T)./dt))*dt; % Timesteps at which the 'sig' is simulated
fprintf('N = %.0f\n',numel(Tsig));
sig = NaN(1,length(Tsig)); % Preallocate the matrix for the signal
sig(1) = observable(rh0); % Observable of the initial density matrix
for i = 2:length(Tsig)
    % Propagate the density matrix by dt
    rht = Um*rhp*Up;
    if opt.fixtrace && opt.check; rht = rht/trace(rht); end % Fix roundoff
    rhp = rht;

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