Fix for #162 #85

KomaMRICore/src/datatypes/Sequence.jl

@@ -331,9 +331,9 @@ Generates a trapezoidal waveform vector.
 			B = A
 		end
 		#Trapezoidal waveform
-		aux = (ζ1    .<= t .- delay .< ζ1+ΔT) .* B
+		aux = (ζ1    .< t .- delay .< ζ1+ΔT) .* B
 		if ζ1 != 0
-			aux .+= (0     .<= t .- delay .< ζ1) .* A[1] .* (t .- delay)./ζ1
+			aux .+= (0     .< t .- delay .<= ζ1) .* A[1] .* (t .- delay)./ζ1
 		end
 		if ζ2 !=0
 			aux .+= (ζ1+ΔT .<= t .- delay .< ζ1+ΔT+ζ2) .* A[end] .* (1 .- (t.-delay.-ζ1.-ΔT)./ζ2)
@@ -813,4 +813,93 @@ get_M2(seq::Sequence; Δt=1) = begin
 	M2_adc = [M2x_adc M2y_adc M2z_adc]
 	#Final
 	M2, M2_adc
+end
+
+"""
+	SR, SR_adc = get_slew_rate(seq::Sequence; Δt=1)
+
+Outputs the designed slew rate of the Sequence `seq`.
+
+# Arguments
+- `seq`: (`::Sequence`) Sequence struct
+- `Δt`: (`::Real`, `=1`, `[s]`) nominal delta time separation between two time samples
+    for ADC acquisition and Gradients
+
+# Returns
+- `SR`: (`3-column ::Matrix{Float64}`) Slew rate
+- `SR_adc`: (`3-column ::Matrix{Float64}`) Slew rate sampled at ADC points
+"""
+get_slew_rate(seq::Sequence; Δt=1) = begin
+	t, Δt = get_uniform_times(seq, Δt)
+	Gx, Gy, Gz = get_grads(seq, t)
+	G = Dict(1=>[Gx; 0], 2=>[Gy; 0], 3=>[Gz; 0])
+	#kspace
+	Nt = length(t)
+	m2 = zeros(Nt,3)
+	#get_RF_center_breaks
+	idx_rf, rf_type = get_RF_types(seq, t)
+	parts = kfoldperm(Nt,1,type="ordered", breaks=idx_rf)
+	for i = 1:3
+		m2f = 0
+		for (rf, p) in enumerate(parts)
+			m2[p,i] = (G[i][p[1]+1:p[end]+1] .- G[i][p[1]:p[end]])[:] ./ Δt[p]
+		end
+	end
+	M2 = m2 #[m^-1]
+	#Interp, as Gradients are generally piece-wise linear, the integral is piece-wise cubic
+	#Nevertheless, the integral is sampled at the ADC times so a linear interp is sufficient
+	#TODO: check if this interpolation is necessary
+	ts = t .+ Δt
+	t_adc =  get_adc_sampling_times(seq)
+	M2x_adc = linear_interpolation(ts,M2[:,1],extrapolation_bc=0)(t_adc)
+	M2y_adc = linear_interpolation(ts,M2[:,2],extrapolation_bc=0)(t_adc)
+	M2z_adc = linear_interpolation(ts,M2[:,3],extrapolation_bc=0)(t_adc)
+	M2_adc = [M2x_adc M2y_adc M2z_adc]
+	#Final
+	M2, M2_adc
+end
+
+"""
+    EC, EC_adc = get_eddy_currents(seq::Sequence; Δt=1)
+
+Outputs the designed eddy currents of the Sequence `seq`.
+
+# Arguments
+- `seq`: (`::Sequence`) Sequence struct
+- `Δt`: (`::Real`, `=1`, `[s]`) nominal delta time separation between two time samples
+    for ADC acquisition and Gradients
+
+# Returns
+- `EC`: (`3-column ::Matrix{Float64}`) Eddy currents
+- `EC_adc`: (`3-column ::Matrix{Float64}`) Eddy currents sampled at ADC points
+"""
+get_eddy_currents(seq::Sequence; Δt=1, λ=80e-3) = begin
+	t, Δt = get_uniform_times(seq, Δt)
+	Gx, Gy, Gz = get_grads(seq, t)
+	G = Dict(1=>[Gx; 0], 2=>[Gy; 0], 3=>[Gz; 0])
+	#kspace
+	Nt = length(t)
+	m2 = zeros(Nt,3)
+	#get_RF_center_breaks
+	idx_rf, rf_type = get_RF_types(seq, t)
+	parts = kfoldperm(Nt,1,type="ordered", breaks=idx_rf)
+	for i = 1:3
+		m2f = 0
+		for (rf, p) in enumerate(parts)
+			m2[p,i] = (G[i][p[1]+1:p[end]+1] .- G[i][p[1]:p[end]])[:] ./ Δt[p]
+		end
+	end
+	ec(t, λ) = exp.(-t ./ λ) .* (t .>= 0)
+	M2 = [sum( m2[:, j] .* ec(t[i] .- t, λ) .* Δt ) for i = 1:Nt, j = 1:3] #[m^-1]
+	#Interp, as Gradients are generally piece-wise linear, the integral is piece-wise cubic
+	#Nevertheless, the integral is sampled at the ADC times so a linear interp is sufficient
+	#TODO: check if this interpolation is necessary
+	ts = t .+ Δt
+	t_adc =  get_adc_sampling_times(seq)
+	M2x_adc = linear_interpolation(ts,M2[:,1],extrapolation_bc=0)(t_adc)
+	M2y_adc = linear_interpolation(ts,M2[:,2],extrapolation_bc=0)(t_adc)
+	M2z_adc = linear_interpolation(ts,M2[:,3],extrapolation_bc=0)(t_adc)
+	M2_adc = [M2x_adc M2y_adc M2z_adc]
+	#Final
+	M2*1_000, M2_adc*1_000
 end

KomaMRICore/src/simulation/TimeStepCalculation.jl

@@ -1,3 +1,4 @@
+const EPS = eps(1.0)
 """
     array_of_ranges = kfoldperm(N, k; type="random", breaks=[])
 
@@ -178,23 +179,19 @@ This function returns non-uniform time points that are relevant in the sequence
 """
 function get_variable_times(seq; dt=1e-3, dt_rf=1e-5)
 	t = Float64[]
+	ϵ = EPS #Smallest Float64
 	ΔT = durs(seq) #Duration of sequence block
 	T0 = cumsum([0; ΔT[:]]) #Start time of each block
 	for i = 1:length(seq)
 		s = seq[i] #Current sequence block
 		t0 = T0[i]
-		if is_ADC_on(s)
-			ts = get_adc_sampling_times(s) .+ t0
-			taux = points_from_key_times(ts; dt) # ADC sampling
-			append!(t, taux)
-		end
 		if is_RF_on(s)
 			y = s.RF[1]
 			delay, T = y.delay, y.T
 			t1 = t0 + delay
 			t2 = t1 + sum(T)
 			rf0 = t0 + get_RF_center(y) #get_RF_center includes delays
-			taux = points_from_key_times([t1,rf0,t2]; dt=dt_rf) # Arbitrary RF
+			taux = points_from_key_times([t1,t1+ϵ,rf0,t2-ϵ,t2]; dt=dt_rf) # Arbitrary RF. Points (t1+ϵ, t2-ϵ) added to fix bug with ADCs
 			append!(t, taux)
 		end
 		if is_GR_on(s)
@@ -204,13 +201,19 @@ function get_variable_times(seq; dt=1e-3, dt_rf=1e-5)
 			if is_Gz_on(s) append!(active_gradients, s.GR.z) end
 			for y = active_gradients
 				ts = get_theo_t(y) .+ t0
-				taux = points_from_key_times(ts; dt)
+				taux = points_from_key_times([ts[1]+ϵ; ts; ts[end]-ϵ]; dt) #The ±ϵ fixes #
 				append!(t, taux)
 			end
 		end
 	end
+	#Adding ADC samples, and removing repeated points
 	tadc = get_adc_sampling_times(seq)
 	t = sort(unique([t; tadc])) #Removing repeated points
+	#Fixes a problem with ADC at the start and end of the seq
+	t0 = t[1]   - ϵ
+	tf = t[end] + ϵ
+	t = [t0; t; tf]
+	#Final time points
 	Δt = t[2:end] .- t[1:end-1]
 	t = t[1:end-1]
 	t, Δt
@@ -237,15 +240,16 @@ Thus, it is possible to split the simulation into excitation and preccesion comp
 function get_breaks_in_RF_key_points(seq::Sequence, t)
 	T0 = cumsum([0; durs(seq)[:]])
 	# Identifying RF key points
+	ϵ = EPS #Smallest Float64
 	key_points = Float64[]
 	key_idxs = Int[]
 	for (i, s) = enumerate(seq)
 		if is_RF_on(s)
 			t0 = T0[i] + s.RF.delay[1]	#start of RF waverform
 			tf = T0[i] + s.RF.dur[1]	#end of RF waveform
 			append!(key_points, [t0; tf])
-			idx0 = argmin(abs.(t.-t0))
-			idxf = argmin(abs.(t.-tf))
+			idx0 = argmin(abs.(t.-(t0+ϵ)))
+			idxf = argmin(abs.(t.-(tf-ϵ)))
 			append!(key_idxs,   [idx0; idxf])
 		end
 	end