qval and bval should be in units of rad

src/gradients/changing_gradient_blocks.jl

@@ -38,7 +38,7 @@ grad_start(cgb::AbstractChangingGradientBlock) = cgb.gradient_strength_start
 slew_rate(cgb::AbstractChangingGradientBlock) = cgb.slew_rate
 grad_end(cgb::AbstractChangingGradientBlock) = grad_start(cgb) .+ slew_rate(cgb) .* duration(cgb)
 gradient_strength(cgb::AbstractChangingGradientBlock) = max.(grad_start(cgb), grad_end(cgb))
-qvec(cgb::AbstractChangingGradientBlock) = (grad_start(cgb) .+ grad_end(cgb)) .* (duration(cgb)/2)
+qvec(cgb::AbstractChangingGradientBlock) = (grad_start(cgb) .+ grad_end(cgb)) .* (duration(cgb) * π)
 
 function bmat_gradient(cgb::AbstractChangingGradientBlock, qstart)
     # grad = (g1 * (duration - t) + g2 * t) / duration
@@ -55,14 +55,14 @@ function bmat_gradient(cgb::AbstractChangingGradientBlock, qstart)
         qstart .* permutedims(qstart) .* duration(cgb) .+
         duration(cgb)^2 .* qstart .* permutedims(
             2 .* grad_start(cgb) .+
-            grad_end(cgb)) ./ 3 .+
+            grad_end(cgb)) .* 2π ./ 3 .+
         bmat_gradient(cgb)
     )
 end
 
 function bmat_gradient(cgb::AbstractChangingGradientBlock)
     diff = slew_rate(cgb) .* duration(cgb)
-    return (
+    return (2π)^2 .* (
         grad_start(cgb) .* permutedims(grad_start(cgb)) ./ 3 .+
         grad_start(cgb) .* permutedims(diff) ./ 4 .+
         diff .* permutedims(diff) ./ 10

src/gradients/constant_gradient_blocks.jl

@@ -35,10 +35,10 @@ const FixedConstantGradientBlock = AbstractConstantGradientBlock{Float64}
 duration(cgb::AbstractConstantGradientBlock) = cgb.duration
 gradient_strength(cgb::AbstractConstantGradientBlock) = cgb.gradient_strength
 slew_rate(::AbstractConstantGradientBlock) = zero(SVector{3, Float64})
-qvec(cgb::AbstractConstantGradientBlock) = duration(cgb) .* gradient_strength(cgb)
+qvec(cgb::AbstractConstantGradientBlock) = duration(cgb) .* gradient_strength(cgb) .* 2π
 
 function bmat_gradient(cgb::AbstractConstantGradientBlock)
-    grad = gradient_strength(cgb)
+    grad = 2π .* gradient_strength(cgb)
     return (grad .* permutedims(grad)) .* duration(cgb)^3 ./3
 end
 
@@ -47,7 +47,7 @@ function bmat_gradient(cgb::AbstractConstantGradientBlock, qstart)
     #   qstart^2 * duration +
     #   qstart * grad * duration^2 +
     #   grad * grad * duration^3 / 3 +
-    grad = gradient_strength(cgb)
+    grad = 2π .* gradient_strength(cgb)
     return (
         qstart .* permutedims(qstart) .* duration(cgb) .+
         qstart .* permutedims(grad) .* duration(cgb)^2 .+

src/gradients/pulsed_gradients.jl

@@ -107,7 +107,7 @@ gradient_strength(g::PulsedGradient) = gradient_strength(g.flat)
 slew_rate(g::PulsedGradient) = g.slew_rate
 δ(g::PulsedGradient) = rise_time(g) + flat_time(g)
 duration(g::PulsedGradient) = 2 * rise_time(g) + flat_time(g)
-qvec(g::PulsedGradient) = δ(g) .* gradient_strength(g)
+qvec(g::PulsedGradient) = δ(g) .* gradient_strength(g) .* 2π
 
 get_children_indices(::PulsedGradient) = (:rise, :flat, :fall)
 Base.getindex(pg::PulsedGradient, symbol::Symbol) = pg[Val(symbol)]

src/pathways.jl

@@ -143,7 +143,7 @@ area_under_curve(pathway::Pathway; scale=nothing, rotate=nothing) = qval(pathway
 """
     bmat(pathway::Pathway; scale=nothing, rotate=nothing)
 
-Return 3x3 diffusion-weighted matrix experienced by the spins following a specific [`Pathway`](@ref).
+Return 3x3 diffusion-weighted matrix experienced by the spins following a specific [`Pathway`](@ref) in rad^2 ms/um^2.
 
 Only gradients active while the spins are in the transverse plane are considered.
 
@@ -155,12 +155,12 @@ bmat(pathway::Pathway; scale=nothing, rotate=nothing)  = get(pathway.bmat, (scal
 """
     bval(pathway::Pathway; scale=nothing, rotate=nothing)
 
-Return size of diffusion-weighting experienced by the spins following a specific [`Pathway`](@ref).
+Return size of diffusion-weighting experienced by the spins following a specific [`Pathway`](@ref) in rad^2 ms/um^2.
 
-Only gradients active while the spins are in the transverse plane are considered.
+Only gradients active while the spins are in the transverse plane will contribute to the diffusion weighting.
 
 By default gradients that are affected by user-provided `scale` or `rotate` parameters (e.g., bvals/bvecs) are ignored.
-You can set `scale` and/or `rotate` to specific symbols to only consider gradients that are affected by speficic `scale`/`rotate` parameters
+You can set `scale` and/or `rotate` to specific symbols to only consider gradients that are affected by specific `scale`/`rotate` parameters
 """
 bval(pathway::Pathway; scale=nothing, rotate=nothing) = tr(bmat(pathway; scale, rotate))
 

src/variables.jl

@@ -19,8 +19,8 @@ all_variables_symbols = [
     :slice_thickness => (:pulse, "Slice thickness of an RF pulse that is active during a gradient."),
 
     # gradients
-    :qvec => (:gradient, "The spatial range and orientation on which the displacements can be detected due to this gradient in 1/um."),
-    :qval => (:gradient, "The spatial range on which the displacements can be detected due to this gradient in 1/um."),
+    :qvec => (:gradient, "The spatial range and orientation on which the displacements can be detected due to this gradient in rad/um."),
+    :qval => (:gradient, "The spatial range on which the displacements can be detected due to this gradient in rad/um (i.e., norm of [`qvec`](@ref))."),
     :δ => (:gradient, "Effective duration of a gradient pulse ([`rise_time`](@ref) + [`flat_time`](@ref)) in ms."),
     :rise_time => (:gradient, "Time for gradient pulse to reach its maximum value in ms."),
     :flat_time => (:gradient, "Time of gradient pulse at maximum value in ms."),
@@ -101,7 +101,7 @@ end
 """
     bmat_gradient(gradient::GradientBlock, qstart=(0, 0, 0))
 
-Computes the diffusion-weighting matrix due to a single gradient block.
+Computes the diffusion-weighting matrix due to a single gradient block in rad^2 ms/um^2.
 
 This should be defined for every `GradientBlock`, but not be called directly.
 Instead, the `bmat` and `bval` should be constrained for specific `Pathways`