R/pulse_sequences.R

In spant: MR Spectroscopy Analysis Tools

#' Simple pulse and acquire sequence with ideal pulses.
#' @param spin_params spin system definition.
#' @param ft transmitter frequency in Hz.
#' @param ref reference value for ppm scale.
#' @param nuc acquisition nucleus.
#' @param acq_delay delay between excitation and acquisition.
#' @return list of resonance amplitudes and frequencies.
#' @export
seq_pulse_acquire <- function(spin_params, ft, ref, nuc = "1H",
                              acq_delay = 0) {
  
  sys <- spin_sys(spin_params, ft, ref)
  
  sys$rho <- gen_F(sys, "z", nuc)
   
  angle <- 90
  Fx <- gen_F(sys, "x",  nuc)
  lhs_pulse <- matexp(-Fx * 1i * angle * pi / 180)
  rhs_pulse <- matexp( Fx * 1i * angle * pi / 180)
  sys$rho <- lhs_pulse %*% sys$rho %*% rhs_pulse
  
  if (acq_delay > 0) {
    # apply delay
    t <- acq_delay
    # find the inverse of the eigenvector matrix
    eig_vec_inv <- solve(sys$H_eig_vecs)
    lhs <- sys$H_eig_vecs %*% diag(exp(sys$H_eig_vals * 2i * pi * t)) %*%
      eig_vec_inv
    rhs <- sys$H_eig_vecs %*% diag(exp(-sys$H_eig_vals * 2i * pi * t)) %*%
      eig_vec_inv
    sys$rho <- lhs %*% sys$rho %*% rhs
  }
  
  # acquire
  acquire(sys, detect = nuc)
}

#' PRESS sequence with ideal pulses.
#' @param spin_params spin system definition.
#' @param ft transmitter frequency in Hz.
#' @param ref reference value for ppm scale.
#' @param TE1 TE1 sequence parameter in seconds (TE=TE1+TE2).
#' @param TE2 TE2 sequence parameter in seconds.
#' @return list of resonance amplitudes and frequencies.
#' @export
seq_press_ideal <- function(spin_params, ft, ref, TE1 = 0.01, TE2 = 0.02) {
  
  sys <- spin_sys(spin_params, ft, ref)
  sys$rho <- -gen_F(sys, "y", "1H")
  
  # apply delay
  t <- TE1 / 2
  
  # find the inverse of the eigenvector matrix
  eig_vec_inv <- solve(sys$H_eig_vecs)
  lhs <- sys$H_eig_vecs %*% diag(exp(sys$H_eig_vals * 2i * pi * t)) %*%
         eig_vec_inv
  
  rhs <- sys$H_eig_vecs %*% diag(exp(-sys$H_eig_vals * 2i * pi * t)) %*%
         eig_vec_inv
  
  sys$rho <- lhs %*% sys$rho %*% rhs
  
  # apply pulse
  angle <- 180
  Fy <- gen_F(sys, "y", "1H")
  lhs_pulse <- matexp(-Fy * 1i * angle * pi / 180)
  rhs_pulse <- matexp( Fy * 1i * angle * pi / 180)
  sys$rho <- lhs_pulse %*% sys$rho %*% rhs_pulse
  
  # apply delay
  t <- (TE1 + TE2) / 2
  lhs <- sys$H_eig_vecs %*% diag(exp(sys$H_eig_vals * 2i * pi * t)) %*% 
         eig_vec_inv
  
  rhs <- sys$H_eig_vecs %*% diag(exp(-sys$H_eig_vals * 2i * pi * t)) %*% 
         eig_vec_inv
  
  sys$rho <- lhs %*% sys$rho %*% rhs
  
  # apply pulse
  sys$rho <- lhs_pulse %*% sys$rho %*% rhs_pulse
  
  # apply delay
  t <- TE2 / 2
  lhs <- sys$H_eig_vecs %*% diag(exp(sys$H_eig_vals * 2i * pi * t)) %*%
         eig_vec_inv
  
  rhs <- sys$H_eig_vecs %*% diag(exp(-sys$H_eig_vals * 2i * pi * t)) %*%
         eig_vec_inv
  
  sys$rho <- lhs %*% sys$rho %*% rhs
  
  # acquire
  acquire(sys, detect = "1H")
}

#' Spin echo sequence with ideal pulses.
#' @param spin_params spin system definition.
#' @param ft transmitter frequency in Hz.
#' @param ref reference value for ppm scale.
#' @param nuc acquisition nucleus.
#' @param TE echo time in seconds.
#' @return list of resonance amplitudes and frequencies.
#' @export
seq_spin_echo_ideal <- function(spin_params, ft, ref, nuc = "1H", TE = 0.03) {
  sys <- spin_sys(spin_params, ft, ref)
  
  sys$rho <- gen_F(sys, "z")
   
  angle <- 90
  Fx <- gen_F(sys, "x", nuc)
  lhs_pulse <- matexp(-Fx * 1i * angle * pi / 180)
  rhs_pulse <- matexp( Fx * 1i * angle * pi / 180)
  sys$rho <- lhs_pulse %*% sys$rho %*% rhs_pulse
  
  # apply delay
  t = TE / 2
  # find the inverse of the eigenvector matrix
  eig_vec_inv <- solve(sys$H_eig_vecs)
  lhs <- sys$H_eig_vecs %*% diag(exp(sys$H_eig_vals * 2i * pi * t)) %*% 
         eig_vec_inv
  
  rhs <- sys$H_eig_vecs %*% diag(exp(-sys$H_eig_vals * 2i * pi * t)) %*%
         eig_vec_inv
  
  sys$rho <- lhs %*% sys$rho %*% rhs
  
  # apply pulse
  angle <- 180
  Fy <- gen_F(sys, "y", nuc)
  lhs_pulse <- matexp(-Fy * 1i * angle * pi / 180)
  rhs_pulse <- matexp( Fy * 1i * angle * pi / 180)
  sys$rho <- lhs_pulse %*% sys$rho %*% rhs_pulse
  
  # apply delay
  sys$rho <- lhs %*% sys$rho %*% rhs
  
  # acquire
  acquire(sys, detect = nuc)
}

#' CPMG style sequence with ideal pulses.
#' @param spin_params spin system definition.
#' @param ft transmitter frequency in Hz.
#' @param ref reference value for ppm scale.
#' @param TE echo time in seconds.
#' @param echoes number of echoes.
#' @return list of resonance amplitudes and frequencies.
#' @export
seq_cpmg_ideal <- function(spin_params, ft, ref, TE = 0.03, echoes = 4) {
  sys <- spin_sys(spin_params, ft, ref)
  sys$rho <- -gen_F(sys, "y", "1H")
  
  # calc delay operator
  t <- TE / (echoes * 2)
  # find the inverse of the eigenvector matrix
  eig_vec_inv <- solve(sys$H_eig_vecs)
  lhs <- sys$H_eig_vecs %*% diag(exp(sys$H_eig_vals * 2i * pi * t)) %*%
         eig_vec_inv
  
  rhs <- sys$H_eig_vecs %*% diag(exp(-sys$H_eig_vals * 2i * pi * t)) %*%
         eig_vec_inv
  
  # calc pulse operator
  angle <- 180
  Fy <- gen_F(sys, "y", "1H")
  lhs_pulse <- matexp(-Fy * 1i * angle * pi / 180)
  rhs_pulse <- matexp( Fy * 1i * angle * pi / 180)
    
  sys$rho <- lhs %*% sys$rho %*% rhs
  sys$rho <- lhs_pulse %*% sys$rho %*% rhs_pulse
  if (echoes > 1) {
    for (n in 1:(echoes - 1)) {
      sys$rho <- lhs %*% sys$rho %*% rhs
      sys$rho <- lhs %*% sys$rho %*% rhs
      sys$rho <- lhs_pulse %*% sys$rho %*% rhs_pulse
    }
  }
    
  sys$rho <- lhs %*% sys$rho %*% rhs
  
  # acquire
  acq_res <- acquire(sys, detect = "1H")
}

#' MEGA-PRESS sequence with ideal localisation pulses and Gaussian shaped
#' editing pulse.
#' @param spin_params spin system definition.
#' @param ft transmitter frequency in Hz.
#' @param ref reference value for ppm scale.
#' @param ed_freq editing pulse frequency in ppm.
#' @param TE1 TE1 sequence parameter in seconds (TE=TE1+TE2).
#' @param TE2 TE2 sequence parameter in seconds.
#' @param BW editing pulse bandwidth in Hz.
#' @param steps number of hard pulses used to approximate the editing pulse.
#' @return list of resonance amplitudes and frequencies.
#' @export
seq_mega_press_ideal <- function(spin_params, ft, ref, ed_freq = 1.89,
                                 TE1 = 0.015, TE2 = 0.053, BW = 110,
                                 steps = 50) {
  
  # ed pulse duration 
  duration <- 1.53 / BW # 180 deg 1% Guass (p245 de Graff)
  
  sys <- spin_sys(spin_params, ft, ed_freq)
  sys$rho <- -gen_F(sys, "y", "1H")
  
  # filter -1
  sys$rho <- coherence_filter(sys, sys$rho, -1)
  
  # apply delay
  t <- TE1 / 2
  
  # find the inverse of the eigenvector matrix
  eig_vec_inv <- solve(sys$H_eig_vecs)
  lhs <- sys$H_eig_vecs %*% diag(exp(sys$H_eig_vals * 2i * pi * t)) %*%
    eig_vec_inv
  
  rhs <- sys$H_eig_vecs %*% diag(exp(-sys$H_eig_vals * 2i * pi * t)) %*%
    eig_vec_inv
  
  sys$rho <- lhs %*% sys$rho %*% rhs
  
  # apply pulse
  angle <- 180
  Fy <- gen_F(sys, "y", "1H")
  lhs_pulse <- matexp(-Fy * 1i * angle * pi / 180)
  rhs_pulse <- matexp( Fy * 1i * angle * pi / 180)
  sys$rho <- lhs_pulse %*% sys$rho %*% rhs_pulse
  
  # filter 1
  sys$rho <- coherence_filter(sys, sys$rho, 1)
  
  # apply delay
  t <- (TE1 + TE2 / 2) / 2 - duration / 2
  if (t < 0) stop("Error, negative delay duration required.")
  lhs <- sys$H_eig_vecs %*% diag(exp(sys$H_eig_vals * 2i * pi * t)) %*% 
    eig_vec_inv
  
  rhs <- sys$H_eig_vecs %*% diag(exp(-sys$H_eig_vals * 2i * pi * t)) %*% 
    eig_vec_inv
  
  sys$rho <- lhs %*% sys$rho %*% rhs
  
  # apply selective 180
  pulse <- get_guassian_pulse(180, steps)
  dt <- duration / steps
  
  eig_vec_inv <- solve(sys$H_eig_vecs)
  lhs_dt <- sys$H_eig_vecs %*% diag(exp(sys$H_eig_vals * 2i * pi * dt)) %*%
    eig_vec_inv
  
  rhs_dt <- sys$H_eig_vecs %*% diag(exp(-sys$H_eig_vals * 2i * pi * dt)) %*%
         eig_vec_inv
  
  for (x in pulse) {
    angle <- x
    Fx <- gen_F(sys, "y", "1H")
    lhs_pulse_gaus <- matexp(-Fx * 1i * angle)
    rhs_pulse_gaus <- matexp( Fx * 1i * angle)
    sys$rho <- lhs_pulse_gaus %*% sys$rho %*% rhs_pulse_gaus
    sys$rho <- lhs_dt %*% sys$rho %*% rhs_dt
  }
  
  # filter 1
  sys$rho <- coherence_filter(sys, sys$rho, 1)
  
  # apply TE2/4 delay
  t <- TE2 / 4 - duration / 2
  if (t < 0) stop("Error, negative delay duration required.")
  lhs <- sys$H_eig_vecs %*% diag(exp(sys$H_eig_vals * 2i * pi * t)) %*%
         eig_vec_inv
  
  rhs <- sys$H_eig_vecs %*% diag(exp(-sys$H_eig_vals * 2i * pi * t)) %*%
         eig_vec_inv
  
  sys$rho <- lhs %*% sys$rho %*% rhs
  
  # apply pulse
  sys$rho <- lhs_pulse %*% sys$rho %*% rhs_pulse
  
  # filter -1
  sys$rho <- coherence_filter(sys, sys$rho, -1)
  
  # apply TE2/4 delay
  sys$rho <- lhs %*% sys$rho %*% rhs
  
  # apply selective 180 
  for (x in pulse) {
    angle <- x
    Fx <- gen_F(sys, "x", "1H")
    lhs_pulse_gaus <- matexp(-Fx * 1i * angle)
    rhs_pulse_gaus <- matexp( Fx * 1i * angle)
    sys$rho <- lhs_pulse_gaus %*% sys$rho %*% rhs_pulse_gaus
    sys$rho <- lhs_dt %*% sys$rho %*% rhs_dt
  }
  
  # filter -1
  sys$rho <- coherence_filter(sys, sys$rho, -1)
  
  # apply TE2/4 delay
  sys$rho <- lhs %*% sys$rho %*% rhs

  # acquire
  acq_res <- acquire(sys, detect = "1H")
  
  # shift back to requested ref
  acq_res$freqs <- acq_res$freqs + (ref - ed_freq) * ft / 1e6
  acq_res
}

#' STEAM sequence with ideal pulses.
#' @param spin_params spin system definition.
#' @param ft transmitter frequency in Hz.
#' @param ref reference value for ppm scale.
#' @param TE sequence parameter in seconds.
#' @param TM sequence parameter in seconds.
#' @param amp_scale amplitude scaling factor. Set to 2 (default) to ensure
#' correct scaling for water reference scaling. Set to 1 to maintain the
#' inherent loss of signal associated with STEAM.
#' @return list of resonance amplitudes and frequencies.
#' @export
seq_steam_ideal <- function(spin_params, ft, ref, TE = 0.03, TM = 0.02,
                            amp_scale = 2) {
  
  sys <- spin_sys(spin_params, ft, ref)
  sys$rho <- gen_F(sys, "z")
  
  # TE/2 delay operator 
  eig_vec_inv <- solve(sys$H_eig_vecs)
  t <- (TE / 2)
  lhs_half_te <- sys$H_eig_vecs %*% diag(exp(sys$H_eig_vals * 2i * pi * t)) %*% 
    eig_vec_inv
  
  rhs_half_te <- sys$H_eig_vecs %*% diag(exp(-sys$H_eig_vals * 2i * pi * t)) %*% 
    eig_vec_inv
  
  # TM delay operator 
  t <- TM
  lhs_tm <- sys$H_eig_vecs %*% diag(exp(sys$H_eig_vals * 2i * pi * t)) %*% 
    eig_vec_inv
  
  rhs_tm <- sys$H_eig_vecs %*% diag(exp(-sys$H_eig_vals * 2i * pi * t)) %*% 
    eig_vec_inv
  
  # 90 degree pulse operator
  angle <- 90
  Fx <- gen_F(sys, "x", "1H")
  lhs_x_pulse <- matexp(-Fx * 1i * angle * pi / 180)
  rhs_x_pulse <- matexp( Fx * 1i * angle * pi / 180)
  
  basis_size   <- prod(sys$spin_n * 2 + 1)
  rho_combined <- matrix(0, basis_size, basis_size)
  
  # phase cycling loop
  for (n in 0:3) {
    phase <- n * 90
    
    # first and third 90 pulse operator
    angle <- 90
    Fxy <- gen_F_xy(sys, phase, "1H")
    
    # XY rotation operator
    lhs_xy_pulse <- matexp(-Fxy * 1i * angle * pi / 180)
    rhs_xy_pulse <- matexp( Fxy * 1i * angle * pi / 180)
    
    # first 90 plus rotation
    rho <- lhs_xy_pulse %*% sys$rho %*% rhs_xy_pulse
    
    # evolve TE/2
    rho <- lhs_half_te %*% rho %*% rhs_half_te
    
    # second 90
    rho <- lhs_x_pulse %*% rho %*% rhs_x_pulse
    
    # zero non-zero-order coherences
    rho <- coherence_filter(sys, rho)
    
    # evolve TM
    rho <- lhs_tm %*% rho %*% rhs_tm
    
    # third 90 plus rotation
    rho <- lhs_xy_pulse %*% rho %*% rhs_xy_pulse
    
    # evolve TE/2
    rho <- lhs_half_te %*% rho %*% rhs_half_te
    
    rho_combined <- rho_combined + rho / 4
  }
  
  sys$rho <- rho_combined
  
  # acquire, and double the output intensity to ensure consistent concentration
  # scaling
  acquire(sys, rec_phase = 180, detect = "1H", amp_scale = amp_scale)
}

#' STEAM sequence with ideal pulses using the z-rotation gradient simulation
#' method described by Young et al JMR 140, 146-152 (1999).
#' @param spin_params spin system definition.
#' @param ft transmitter frequency in Hz.
#' @param ref reference value for ppm scale.
#' @param TE sequence parameter in seconds.
#' @param TM sequence parameter in seconds.
#' @param amp_scale amplitude scaling factor. Set to 2 (default) to ensure
#' correct scaling for water reference scaling. Set to 1 to maintain the
#' inherent loss of signal associated with STEAM.
#' @return list of resonance amplitudes and frequencies.
#' @export
seq_steam_ideal_young <- function(spin_params, ft, ref, TE = 0.03, TM = 0.02,
                                  amp_scale = 2) {
  
  sys <- spin_sys(spin_params, ft, ref)
  Fz  <- gen_F(sys, "z")
  sys$rho <- Fz
  
  # TE/2 delay operator 
  eig_vec_inv <- solve(sys$H_eig_vecs)
  t <- (TE / 2)
  lhs_half_te <- sys$H_eig_vecs %*% diag(exp(sys$H_eig_vals * 2i * pi * t)) %*% 
    eig_vec_inv
  
  rhs_half_te <- sys$H_eig_vecs %*% diag(exp(-sys$H_eig_vals * 2i * pi * t)) %*% 
    eig_vec_inv
  
  # TM delay operator 
  t <- TM
  lhs_tm <- sys$H_eig_vecs %*% diag(exp(sys$H_eig_vals * 2i * pi * t)) %*% 
    eig_vec_inv
  
  rhs_tm <- sys$H_eig_vecs %*% diag(exp(-sys$H_eig_vals * 2i * pi * t)) %*% 
    eig_vec_inv
  
  # 90 degree x pulse operator
  angle <- 90
  Fx <- gen_F(sys, "x", "1H")
  lhs_x_pulse <- matexp(-Fx * 1i * angle * pi / 180)
  rhs_x_pulse <- matexp( Fx * 1i * angle * pi / 180)
  
  basis_size   <- prod(sys$spin_n * 2 + 1)
  rho_combined <- matrix(0, basis_size, basis_size)
  
  # apply first 90
  sys$rho <- lhs_x_pulse %*% sys$rho %*% rhs_x_pulse
  
  # evolve TE/2
  sys$rho <- lhs_half_te %*% sys$rho %*% rhs_half_te
  
  # phase cycling loop
  for (n in 0:3) {
    phase_ang <- n * 90
    
    # rotate about z
    lhs_z_rot <- matexp(-Fz * 1i * phase_ang * pi / 180)
    rhs_z_rot <- matexp( Fz * 1i * phase_ang * pi / 180)
    rho <- lhs_z_rot %*% sys$rho %*% rhs_z_rot
    
    # second 90
    rho <- lhs_x_pulse %*% rho %*% rhs_x_pulse
    
    # zero non-zero-order coherences
    rho <- coherence_filter(sys, rho)
    
    # evolve TM
    rho <- lhs_tm %*% rho %*% rhs_tm
    
    # third 90
    rho <- lhs_x_pulse %*% rho %*% rhs_x_pulse
    
    # rotate about z
    rho <- lhs_z_rot %*% rho %*% rhs_z_rot

    rho_combined <- rho_combined + rho / 4
  }
  
  # evolve TE/2
  sys$rho <- lhs_half_te %*% rho_combined %*% rhs_half_te
  
  # acquire, and double the output intensity to ensure consistent concentration
  # scaling
  acquire(sys, detect = "1H", rec_phase = 180, amp_scale = amp_scale)
}

#' STEAM sequence with ideal pulses and coherence order filtering to simulate
#' gradient crushers.
#' 
#' See Landheer et al NMR Biomed 2021 34(5):e4129 and Landheer et al MRM 2019
#' Apr;81(4):2209-2222 for more details on the coherence order filtering method.
#' 
#' @param spin_params spin system definition.
#' @param ft transmitter frequency in Hz.
#' @param ref reference value for ppm scale.
#' @param TE sequence parameter in seconds.
#' @param TM sequence parameter in seconds.
#' @param amp_scale amplitude scaling factor. Set to 2 (default) to ensure
#' correct scaling for water reference scaling. Set to 1 to maintain the
#' inherent loss of signal associated with STEAM.
#' @return list of resonance amplitudes and frequencies.
#' @export
seq_steam_ideal_cof <- function(spin_params, ft, ref, TE = 0.03, TM = 0.02,
                                amp_scale = 2) {
  
  sys <- spin_sys(spin_params, ft, ref)
  Fz  <- gen_F(sys, "z")
  sys$rho <- Fz
  
  # TE/2 delay operator 
  eig_vec_inv <- solve(sys$H_eig_vecs)
  t <- (TE / 2)
  lhs_half_te <- sys$H_eig_vecs %*% diag(exp(sys$H_eig_vals * 2i * pi * t)) %*% 
    eig_vec_inv
  
  rhs_half_te <- sys$H_eig_vecs %*% diag(exp(-sys$H_eig_vals * 2i * pi * t)) %*% 
    eig_vec_inv
  
  # TM delay operator 
  t <- TM
  lhs_tm <- sys$H_eig_vecs %*% diag(exp(sys$H_eig_vals * 2i * pi * t)) %*% 
    eig_vec_inv
  
  rhs_tm <- sys$H_eig_vecs %*% diag(exp(-sys$H_eig_vals * 2i * pi * t)) %*% 
    eig_vec_inv
  
  # 90 degree x pulse operator
  angle <- 90
  Fx <- gen_F(sys, "x", "1H")
  lhs_x_pulse <- matexp(-Fx * 1i * angle * pi / 180)
  rhs_x_pulse <- matexp( Fx * 1i * angle * pi / 180)
  
  basis_size   <- prod(sys$spin_n * 2 + 1)
  rho_combined <- matrix(0, basis_size, basis_size)
  
  # apply first 90
  sys$rho <- lhs_x_pulse %*% sys$rho %*% rhs_x_pulse
  
  # filter +1
  sys$rho <- coherence_filter(sys, sys$rho, 1)
  
  # evolve TE/2
  sys$rho <- lhs_half_te %*% sys$rho %*% rhs_half_te
  
  # second 90
  sys$rho <- lhs_x_pulse %*% sys$rho %*% rhs_x_pulse
  
  # filter 0
  sys$rho <- coherence_filter(sys, sys$rho, 0)
  
  # evolve TM
  sys$rho <- lhs_tm %*% sys$rho %*% rhs_tm
  
  # third 90
  sys$rho <- lhs_x_pulse %*% sys$rho %*% rhs_x_pulse
  
  # filter -1
  sys$rho <- coherence_filter(sys, sys$rho, -1)
  
  # evolve TE/2
  sys$rho <- lhs_half_te %*% sys$rho %*% rhs_half_te
  
  # acquire, and double the output intensity to ensure consistent concentration
  # scaling
  acquire(sys, rec_phase = 180, detect = "1H", amp_scale = amp_scale)
}

#' sLASER sequence with ideal pulses.
#' @param spin_params spin system definition.
#' @param ft transmitter frequency in Hz.
#' @param ref reference value for ppm scale.
#' @param TE1 first echo time (between exc. and 1st echo) in seconds.
#' @param TE2 second echo time (between 2nd echo and 4th echo) in seconds.
#' @param TE3 third echo time (between 4th echo and 5th echo) in seconds.
#' @return list of resonance amplitudes and frequencies.
#' @export
seq_slaser_ideal <- function(spin_params, ft, ref, TE1 = 0.008, TE2 = 0.011,
                             TE3 = 0.009) {
  
  sys <- spin_sys(spin_params, ft, ref)
  
  sys$rho <- gen_F(sys, "z")
  
  angle <- 90
  Fx <- gen_F(sys, "x", "1H")
  lhs_pulse <- matexp(-Fx * 1i * angle * pi / 180)
  rhs_pulse <- matexp( Fx * 1i * angle * pi / 180)
  sys$rho <- lhs_pulse %*% sys$rho %*% rhs_pulse
  
  # apply delay
  t <- TE1 / 2
  # find the inverse of the eigenvector matrix
  eig_vec_inv <- solve(sys$H_eig_vecs)
  lhs <- sys$H_eig_vecs %*% diag(exp(sys$H_eig_vals * 2i * pi * t)) %*%
    eig_vec_inv
  
  rhs <- sys$H_eig_vecs %*% diag(exp(-sys$H_eig_vals * 2i * pi * t)) %*%
    eig_vec_inv
  sys$rho <- lhs %*% sys$rho %*% rhs
  
  # apply pulse
  angle <- 180
  Fy <- gen_F(sys, "y", "1H")
  lhs_pulse <- matexp(-Fy * 1i * angle * pi / 180)
  rhs_pulse <- matexp( Fy * 1i * angle * pi / 180)
  sys$rho <- lhs_pulse %*% sys$rho %*% rhs_pulse
  
  # apply delay
  t = TE1 / 2 + TE2 / 4
  # find the inverse of the eigenvector matrix
  eig_vec_inv <- solve(sys$H_eig_vecs)
  lhs <- sys$H_eig_vecs %*% diag(exp(sys$H_eig_vals * 2i * pi * t)) %*%
    eig_vec_inv
  
  rhs <- sys$H_eig_vecs %*% diag(exp(-sys$H_eig_vals * 2i * pi * t)) %*%
    eig_vec_inv
  sys$rho <- lhs %*% sys$rho %*% rhs
  
  # apply pulse
  angle <- 180
  Fy <- gen_F(sys, "y", "1H")
  lhs_pulse <- matexp(-Fy * 1i * angle * pi / 180)
  rhs_pulse <- matexp( Fy * 1i * angle * pi / 180)
  sys$rho <- lhs_pulse %*% sys$rho %*% rhs_pulse
  
  # apply delay
  t <- 2 * TE2 / 4
  # find the inverse of the eigenvector matrix
  eig_vec_inv <- solve(sys$H_eig_vecs)
  lhs <- sys$H_eig_vecs %*% diag(exp(sys$H_eig_vals * 2i * pi * t)) %*%
    eig_vec_inv
  
  rhs <- sys$H_eig_vecs %*% diag(exp(-sys$H_eig_vals * 2i * pi * t)) %*%
    eig_vec_inv
  sys$rho <- lhs %*% sys$rho %*% rhs
  
  # apply pulse
  angle <- 180
  Fy <- gen_F(sys, "y", "1H")
  lhs_pulse <- matexp(-Fy * 1i * angle * pi / 180)
  rhs_pulse <- matexp( Fy * 1i * angle * pi / 180)
  sys$rho <- lhs_pulse %*% sys$rho %*% rhs_pulse
  
  # apply delay
  t <- TE2 / 4 + TE3 /2
  # find the inverse of the eigenvector matrix
  eig_vec_inv <- solve(sys$H_eig_vecs)
  lhs <- sys$H_eig_vecs %*% diag(exp(sys$H_eig_vals * 2i * pi * t)) %*%
    eig_vec_inv
  
  rhs <- sys$H_eig_vecs %*% diag(exp(-sys$H_eig_vals * 2i * pi * t)) %*%
    eig_vec_inv
  sys$rho <- lhs %*% sys$rho %*% rhs
  
  # apply pulse
  angle <- 180
  Fy <- gen_F(sys, "y", "1H")
  lhs_pulse <- matexp(-Fy * 1i * angle * pi / 180)
  rhs_pulse <- matexp( Fy * 1i * angle * pi / 180)
  sys$rho <- lhs_pulse %*% sys$rho %*% rhs_pulse
  
  # apply delay
  t <- TE3 /2
  # find the inverse of the eigenvector matrix
  eig_vec_inv <- solve(sys$H_eig_vecs)
  lhs <- sys$H_eig_vecs %*% diag(exp(sys$H_eig_vals * 2i * pi * t)) %*%
    eig_vec_inv
  
  rhs <- sys$H_eig_vecs %*% diag(exp(-sys$H_eig_vals * 2i * pi * t)) %*%
    eig_vec_inv
  sys$rho <- lhs %*% sys$rho %*% rhs
  
  # acquire
  acquire(sys, detect = "1H")
}