hsqc_si.m

function hsqc_si(x,cycle,J)
% Density matrix simulation of sensitivity-enhanced HSQC (not dependent on Spinach).
% Takes two parameters:
%    x, which is t1 in seconds.
%    cycle, which is the phase cycle number (from 1 to 8),
%    J, the C-H coupling in Hz.

% Includes helder functions which plot the trajectory of a density matrix
% on a Bloch sphere. Set J(CH) to 150 Hz to visualise what happens with a
% 13C-1H pair, or 0 Hz for a 12C-1H pair.
% The number of slices (for spatial averaging) can be fairly low, ca. 31.

% Start parallel pool if not already running
p = gcp('nocreate');
if isempty(p)
    parpool();
end

    function plotBlochSphere()
        % Sets up an empty Bloch sphere.
        [X, Y, Z] = sphere(100);
        s = surf(X, Y, Z);
        axis equal
        shading interp
        set(s, 'FaceAlpha', 0.1, 'FaceColor', '#88ee88', 'EdgeColor', 'none')

        line([-1.2 1.2], [0 0], [0 0], 'LineWidth', 1.2, 'Color', '#888888')
        line([0 0], [-1.2 1.2], [0 0], 'LineWidth', 1.2, 'Color', '#888888')
        line([0 0], [0 0], [-1.2 1.2], 'LineWidth', 1.2, 'Color', '#888888')


        text(0, 0, 1.3, '$z$', 'Interpreter', 'latex', 'FontSize', 20, 'HorizontalAlignment', 'Center')
        text(1.3, 0, 0, '$x$', 'Interpreter', 'latex', 'FontSize', 20, 'HorizontalAlignment', 'Center')
        text(0, 1.3, 0, '$y$', 'Interpreter', 'latex', 'FontSize', 20, 'HorizontalAlignment', 'Center')
        view([60 15])
        set(gca, 'visible', 'off')
    end

    function plotDensityMatrixI(rho)
        % Plots a density matrix rho as a red vector.
        cIx = trace(rho * Ix);
        cIy = trace(rho * Iy);
        cIz = trace(rho * Iz);
        quiver3(0, 0, 0, real(cIx), real(cIy), real(cIz), 'LineWidth', 2, 'Color', '#ee2222', 'MaxHeadSize', 0.4);
    end

    function plotReceiver(ph)
        % Plots the word 'rec' along the receiver axis, as determined by ph30.
        xc = ((mod(ph,4) == 0) - (mod(ph,4) == 2)) * 2;
        yc = ((mod(ph,4) == 1) - (mod(ph,4) == 3)) * 2;
        text(xc, yc, 0, 'rec', 'FontSize', 20, 'HorizontalAlignment', 'Center')
    end

    function decompose(rho)
        % Decomposes a given density matrix into its constituent elements.
        threshold = 0.05;
        if abs(trace(rho * Ix)) > threshold
        fprintf("\namount of Ix: %.3f;  ", trace(rho * Ix))
        end
        if abs(trace(rho * Iy)) > threshold
        fprintf("\namount of Iy: %.3f;  ", trace(rho * Iy))
        end
        if abs(trace(rho * Iz)) > threshold
        fprintf("\namount of Iz: %.3f;  ", trace(rho * Iz))
        end
        if abs(trace(rho * Sx)) > threshold
        fprintf("\namount of Sx: %.3f;  ", trace(rho * Sx))
        end
        if abs(trace(rho * Sy)) > threshold
        fprintf("\namount of Sy: %.3f;  ", trace(rho * Sy))
        end
        if abs(trace(rho * Sz)) > threshold
        fprintf("\namount of Sz: %.3f;  ", trace(rho * Sz))
        end
        if abs(trace(rho * (2*Ix*Sx))) > threshold
        fprintf("\namount of 2IxSx: %.3f;  ", trace(rho * (2*Ix*Sx)))
        end
        if abs(trace(rho * (2*Ix*Sy))) > threshold
        fprintf("\namount of 2IxSy: %.3f;  ", trace(rho * (2*Ix*Sy)))
        end
        if abs(trace(rho * (2*Ix*Sz))) > threshold
        fprintf("\namount of 2IxSz: %.3f;  ", trace(rho * (2*Ix*Sz)))
        end
        if abs(trace(rho * (2*Iy*Sx))) > threshold
        fprintf("\namount of 2IySx: %.3f;  ", trace(rho * (2*Iy*Sx)))
        end
        if abs(trace(rho * (2*Iy*Sy))) > threshold
        fprintf("\namount of 2IySy: %.3f;  ", trace(rho * (2*Iy*Sy)))
        end
        if abs(trace(rho * (2*Iy*Sz))) > threshold
        fprintf("\namount of 2IySz: %.3f;  ", trace(rho * (2*Iy*Sz)))
        end
        if abs(trace(rho * (2*Iz*Sx))) > threshold
        fprintf("\namount of 2IzSx: %.3f;  ", trace(rho * (2*Iz*Sx)))
        end
        if abs(trace(rho * (2*Iz*Sy))) > threshold
        fprintf("\namount of 2IzSy: %.3f;  ", trace(rho * (2*Iz*Sy)))
        end
        if abs(trace(rho * (2*Iz*Sz))) > threshold
        fprintf("\namount of 2IzSz: %.3f;  ", trace(rho * (2*Iz*Sz)))
        end
        fprintf("\n\n")
    end


% tic

% Spin system settings
freqs = [67 456789]; % H and C offsets respectively (in Hz)
% J = 150; % manually set 1J(CH) (Hz)
J_IS = J;
Omega_I = 2*pi*freqs(1);
Omega_S = 2*pi*freqs(2);

% Pauli matrices
Sigmax = 0.5*[0,1;1,0];Sigmay = 0.5*[0,-1i;1i,0];Sigmaz = 0.5*[1,0;0,-1];Unity_mat = [1,0;0,1];
Ix = kron(Sigmax,Unity_mat);
Iy = kron(Sigmay,Unity_mat);
Iz = kron(Sigmaz,Unity_mat);
Sx = kron(Unity_mat,Sigmax);
Sy = kron(Unity_mat,Sigmay);
Sz = kron(Unity_mat,Sigmaz);
all_x = Ix + Sx;
all_y = Iy + Sy;
all_z = Iz + Sz;
all_plus = all_x + 1i*all_y;
all_minus = all_x - 1i*all_y;
Ipul = @(ph) ((mod(ph,4) == 0) - (mod(ph,4) == 2))*Ix + ((mod(ph,4) == 1) - (mod(ph,4) ==3))*Iy;
Spul = @(ph) ((mod(ph,4) == 0) - (mod(ph,4) == 2))*Sx + ((mod(ph,4) == 1) - (mod(ph,4) ==3))*Sy;

% GRADIENT SETUP
grad_max = 67.5 ; % maximum gradient on the instrument (G/cm)
gamma_H = 4257.747892; % gamma of proton (Hz/Gauss)
L_sample = 1.8; % length of the active volume (cm)
np_slices = 100;

gpz1 = 40; % gradient z-amplitude
grad_used1 = grad_max * gpz1 /100;
sw_g1 = gamma_H*grad_used1*L_sample;

gpz2 = 20; % gradient z-amplitude
grad_used2 = grad_max * gpz2 /100;
sw_g2 = gamma_H*grad_used2*L_sample;

gpz6 = 31; % gradient z-amplitude
grad_used6 = grad_max * gpz6 /100;
sw_g6 = gamma_H*grad_used6*L_sample;

gpz7 = 46; % gradient z-amplitude
grad_used7 = grad_max * gpz7 /100;
sw_g7 = gamma_H*grad_used7*L_sample;

if np_slices > 1
    Omega_g1 = 2*pi*linspace(-sw_g1/2,sw_g1/2,np_slices); % frequency distribution induced by the gradient
    Omega_g2 = 2*pi*linspace(-sw_g2/2,sw_g2/2,np_slices); % frequency distribution induced by the gradient
    Omega_g6 = 2*pi*linspace(-sw_g6/2,sw_g6/2,np_slices); % frequency distribution induced by the gradient
    Omega_g7 = 2*pi*linspace(-sw_g7/2,sw_g7/2,np_slices); % frequency distribution induced by the gradient
else
    Omega_g1 = (0);
    Omega_g2 = (0);
    Omega_g6 = (0);
    Omega_g7 = (0);
end


% EXPERIMENT SETUP
cnst2 = 150; % in theory should set to J
td = 16384;
sw = 10000;
dw = 1/sw;
fid(td) = 0;

% HAMILTONIAN SETUP
pw_90 = 1e-9; % duration of hard 90-degree pulse - can be the same for I and S
rf_hard = 1/(4*pw_90);
% Anonymous function for unitary evolution
evolve = @(rho_init, H, t) (expm(-1i*H*t) * rho_init * expm(1i*H*t));
% Common Hamiltonians
H_free = Omega_I*Iz + Omega_S*Sz + 2*pi*J_IS*(Iz*Sz); % free precession, weak coupling
H_hard_I = @(ph) 2*pi*rf_hard*Ipul(ph);  % proton hard pulse with phase [0,1,2,3]
H_hard_S = @(ph) 2*pi*rf_hard*Spul(ph);  % carbon hard pulse with phase [0,1,2,3]

figure(1);
plotBlochSphere();

%% PHASE CYCLING

ph1  = [0 0 0 0 0 0 0 0];
ph2  = [1 1 1 1 1 1 1 1];
ph3  = [0 2 0 2 0 2 0 2];
ph4  = [0 0 0 0 2 2 2 2];
ph5  = [0 0 2 2 0 0 2 2];
ph6  = [0 0 0 0 0 0 0 0];
ph7  = [0 0 0 0 2 2 2 2];
ph8  = [1 1 1 1 1 1 1 1];
ph10 = [0 0 0 0 0 0 0 0];
ph30 = [0 2 2 0 0 2 2 0];
% cycle = 2;  % manually select which phase cycle to use - from 1 through 8

%% EXPERIMENT

rho_0 = Iz; % initial state
EA = 1; % -1 for antiecho experiment

% INEPT
rho = evolve(rho_0, H_free + H_hard_I(ph1(cycle)), pw_90);                                      % 90(x) H
rho = evolve(rho, H_free, 1/(4*cnst2));                                                         % 1/(4J) delay
rho = evolve(rho, H_free + H_hard_I(ph1(cycle)) + H_hard_S(ph6(cycle)), pw_90*2);               % 180(x) HC
rho = evolve(rho, H_free, 1/(4*cnst2));                                                         % 1/(4J) delay
rho = evolve(rho, H_free + H_hard_I(ph2(cycle)) + H_hard_S(ph3(cycle)), pw_90);                 % 90(y) H + 90(x) C

% Spatial averaging begins
rho_sumoverslices = zeros(size(rho));
for k=1:np_slices
    rho_k = rho;
    H_gradient1 = Omega_g1(k)*Iz + Omega_g1(k)*Sz/4;                                    % accounting for gyromagnetic ratio of C
    H_gradient2 = Omega_g2(k)*Iz + Omega_g2(k)*Sz/4;
    H_gradient6 = Omega_g6(k)*Iz + Omega_g6(k)*Sz/4;
    H_gradient7 = Omega_g7(k)*Iz + Omega_g7(k)*Sz/4;
    
    % Multiplicity editing, part 1
    rho_k = evolve(rho_k, H_free, 0.01);                                                        % to refocus CS evolution during Gz1
    rho_k = evolve(rho_k, H_free + H_hard_S(ph1(cycle)), pw_90*2);                              % 180(x) C
    rho_k = evolve(rho_k, H_free + H_gradient1*EA, 0.01);                                       % 1 ms gradient Gz1*EA
    
    % t1 period
    rho_k = evolve(rho_k, H_free, x/2);                                                 % t1/2
    rho_k = evolve(rho_k, H_free + H_hard_I(ph4(cycle)), pw_90*2);                      % 180(x) H
    rho_k = evolve(rho_k, H_free, x/2);                                                 % t1/2

    % Multiplicity editing, part 2
    rho_k = evolve(rho_k, H_free, 1/(2*cnst2) - 0.01);                                          % 1/(2J) delay - 1 ms
    rho_k = evolve(rho_k, H_free + H_gradient1*EA, 0.01);                                       % 1 ms gradient Gz1*EA
    rho_k = evolve(rho_k, H_free + H_hard_I(ph2(cycle)) + H_hard_S(ph6(cycle)), pw_90*2);       % 180(y) H + 180(x) C
    rho_k = evolve(rho_k, H_free, 1/(2*cnst2));                                                 % 1/(2J) delay

    % Reverse INEPT
    rho_k = evolve(rho_k, H_free + H_hard_I(ph5(cycle)) + H_hard_S(ph1(cycle)), pw_90);         % 90(x) HC
    rho_k = evolve(rho_k, H_free, 1/(4*cnst2) - 0.01);                                          % 1/(4J) delay - 1 ms
    rho_k = evolve(rho_k, H_free + H_gradient6, 0.01);                                          % 1 ms gradient Gz6
    rho_k = evolve(rho_k, H_free + H_hard_I(ph1(cycle)) + H_hard_S(ph6(cycle)), pw_90*2);       % 180(x) HC
    rho_k = evolve(rho_k, H_free, 1/(4*cnst2) - 0.01);                                          % 1/(4J) delay - 1 ms
    rho_k = evolve(rho_k, H_free + H_gradient6, 0.01);                                          % 1 ms gradient Gz6

    % Sensitivity enhancement block from hsqcedetgpsisp2.3
    rho_k = evolve(rho_k, H_free + H_hard_I(ph2(cycle)) + H_hard_S(ph8(cycle) + 1 - EA), pw_90);    % 90(y)H, 90(+-y)C depending on EA
    rho_k = evolve(rho_k, H_free, 1/(4*cnst2) - 0.01);                                              % 1/(4J) delay - 1 ms
    rho_k = evolve(rho_k, H_free + H_gradient7, 0.01);                                              % 1 ms gradient Gz7
    rho_k = evolve(rho_k, H_free + H_hard_I(ph1(cycle)) + H_hard_S(ph1(cycle)), pw_90*2);           % 180(x) HC
    rho_k = evolve(rho_k, H_free, 1/(4*cnst2) - 0.01);                                              % 1/(4J) delay - 1 ms
    rho_k = evolve(rho_k, H_free + H_gradient7, 0.01);                                              % 1 ms gradient Gz

    % Convert to detectable magnetisation & refocusing gradient
    rho_k = evolve(rho_k, H_free + H_hard_I(ph1(cycle)), pw_90);                    % 90(x) H
    rho_k = evolve(rho_k, H_free, 0.01);                                            % 1 ms delay
    rho_k = evolve(rho_k, H_free + H_hard_I(ph1(cycle)), pw_90*2);                  % 180(x) H
    rho_k = evolve(rho_k, H_free + H_gradient2, 0.01);                              % 1 ms gradient Gz2
    rho_sumoverslices = rho_sumoverslices + rho_k;
    
    hold on
    plotDensityMatrixI(rho_k);
    hold off
end

% Rescaling density matrix
rho_avgoverslices = rho_sumoverslices / np_slices;
fprintf("final magnetisation ");
decompose(rho_avgoverslices);

% Detection
plotReceiver(ph30(cycle));
text(1, 1, 1, sprintf("%d", cycle), 'FontSize', 20, 'HorizontalAlignment', 'Center')
det = (Ix - 1i*Iy) * exp(1i*(pi/2)*ph30(cycle));
rho_det = rho_avgoverslices;
for j=1:td
    fid(j) = trace(det' * rho_det);
    rho_det = evolve(rho_det, H_free, dw);
end

% faux-relaxation (+ exponential line-broadening)
win = exp(-6 * linspace(0, 1, length(fid)));
fid = fid .* win;

% Fourier transform and plotting
si = 32768;
frq = linspace(-sw/2,sw/2,si);
spec = fftshift(fft(fid,si));

% toc


figure(2);
plot(frq,real(spec));

end