From 60208b60c1377bbde758ee07a7c1142dd32cf5c7 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Thu, 12 Aug 2021 14:45:49 -0400 Subject: [PATCH 01/13] update type promotions and flux --- .../Semi-infinite/DAsemiinf.jl | 80 +++++++++++++------ test/DAsemiinfTests/runtests.jl | 4 + 2 files changed, 60 insertions(+), 24 deletions(-) diff --git a/src/forwardmodels/Diffusion Approximation/Semi-infinite/DAsemiinf.jl b/src/forwardmodels/Diffusion Approximation/Semi-infinite/DAsemiinf.jl index 240d877..5d91769 100644 --- a/src/forwardmodels/Diffusion Approximation/Semi-infinite/DAsemiinf.jl +++ b/src/forwardmodels/Diffusion Approximation/Semi-infinite/DAsemiinf.jl @@ -75,8 +75,10 @@ function fluence_DA_semiinf_TD(t, ρ, μa, μsp; n_ext = 1.0, n_med = 1.0, z = 0 if isa(t, AbstractFloat) return _kernel_fluence_DA_semiinf_TD(D, ν, t, μa, z, z0, ρ, zb) - elseif isa(t, AbstractArray) - ϕ = zeros(eltype(ρ), length(t)) + elseif isa(t, AbstractVector) + T = promote_type(eltype(ρ), eltype(μa), eltype(μsp), eltype(z)) + ϕ = zeros(T, length(t)) + @inbounds Threads.@threads for ind in eachindex(t) ϕ[ind] = _kernel_fluence_DA_semiinf_TD(D, ν, t[ind], μa, z, z0, ρ, zb) end @@ -119,32 +121,27 @@ function fluence_DA_semiinf_FD(ρ, μa, μsp; n_ext = 1.0, n_med = 1.0, z = 0.0, end ##################################### -# Time-Domain Reflectance (Flux) +# Time-Domain Flux ##################################### """ - refl_DA_semiinf_TD(t, ρ, μa, μsp, n_ext, n_med) + flux_DA_semiinf_TD(t, ρ, μa, μsp; n_ext = 1.0, n_med = 1.0) -Compute the time-domain reflectance from a semi-infinite medium from Eqn. 36 Contini 97 (m = 0). +Compute the time-domain flux (D*∂ϕ(t)/∂z @ z = 0) from a semi-infinite medium from Eqn. 36 Contini 97 (m = 0). # Arguments - `t`: the time vector (ns). - `ρ`: the source detector separation (cm⁻¹) - `μa`: absorption coefficient (cm⁻¹) - `μsp`: reduced scattering coefficient (cm⁻¹) -- `ndet`: the boundary's index of refraction (air or detector) -- `nmed`: the sample medium's index of refraction +- `n_det`: the boundary's index of refraction (air or detector) +- `n_med`: the sample medium's index of refraction # Examples ```jldoctest -julia> `refl_DA_semiinf_TD(0.2:0.6:2.0, 1.0, 0.1, 10.0, 1.0, 1.0)` -4-element Vector{Float64}: - 0.03126641311563058 - 0.00042932742937583005 - 2.016489075323064e-5 - 1.4467359376429123e-6 +julia> `flux_DA_semiinf_TD(1.0, 1.0, 0.1, 10.0, n_ext = 1.0, n_med = 1.0)` ``` """ -function refl_DA_semiinf_TD(t, ρ, μa, μsp, n_ext, n_med) +function flux_DA_semiinf_TD(t, ρ, μa, μsp; n_ext = 1.0, n_med = 1.0) D = D_coeff(μsp, μa) A = A_coeff(n_med / n_ext) ν = ν_coeff(n_med) @@ -152,18 +149,53 @@ function refl_DA_semiinf_TD(t, ρ, μa, μsp, n_ext, n_med) z0 = z0_coeff(μsp) zb = zb_coeff(A, D) - Rt = Array{eltype(ρ)}(undef, length(t)) - z3m = - z0 z4m = 2 * zb + z0 - Threads.@threads for n in eachindex(t) - tmp1 = 4 * D * ν * t[n] - Rt[n] = -(ρ^2 / (tmp1)) - Rt[n] -= μa * ν * t[n] - Rt[n] = -exp(Rt[n]) / (2 * (4 * π * D * ν)^(3/2) * sqrt(t[n] * t[n] * t[n] * t[n] * t[n])) - Rt[n] *= (z3m * exp(-(z3m^2 / (tmp1))) - z4m * exp(-(z4m^2 / (tmp1)))) + if isa(t, AbstractFloat) + return _kernel_flux_DA_semiinf_TD(D, ν, t, μa, z3m, z4m, ρ) + elseif isa(t, AbstractArray) + T = promote_type(eltype(ρ), eltype(μa), eltype(μsp)) + Rt = zeros(T, length(t)) + @inbounds Threads.@threads for ind in eachindex(t) + Rt[ind] = _kernel_flux_DA_semiinf_TD(D, ν, t[ind], μa, z3m, z4m, ρ) + end + return Rt end - - return Rt + end +@inline function _kernel_flux_DA_semiinf_TD(D, ν, t, μa, z3m, z4m, ρ) + @assert t > zero(eltype(D)) + tmp1 = 4 * D * ν * t + Rt = -(ρ^2 / (tmp1)) + Rt -= μa * ν * t + Rt = -exp(Rt) / (2 * (4 * π * D * ν)^(3/2) * sqrt(t * t * t * t * t)) + Rt *= (z3m * exp(-(z3m^2 / (tmp1))) - z4m * exp(-(z4m^2 / (tmp1)))) + + return Rt +end + + +##################################### +# Steady-State Flux +##################################### +""" + flux_DA_semiinf_CW(ρ, μa, μsp; n_ext = 1.0, n_med = 1.0) + +Compute the steady-state flux (D*∂ϕ(ρ)/∂z @ z = 0) from a semi-infinite medium. + +# Arguments +- `ρ`: the source detector separation (cm⁻¹) +- `μa`: absorption coefficient (cm⁻¹) +- `μsp`: reduced scattering coefficient (cm⁻¹) +- `ndet`: the boundary's index of refraction (air or detector) +- `nmed`: the sample medium's index of refraction + +# Examples +```jldoctest +julia> `flux_DA_semiinf_CW(1.0, 0.1, 10.0, n_ext = 1.0, n_med = 1.0)` +``` +""" +function flux_DA_semiinf_CW(ρ, μa, μsp; n_ext = 1.0, n_med = 1.0) + D = D_coeff(μsp, μa) + return D * ForwardDiff.derivative(z -> fluence_DA_semiinf_CW(ρ, μa, μsp, n_med = n_med, n_ext = n_ext, z = z), 0.0) end diff --git a/test/DAsemiinfTests/runtests.jl b/test/DAsemiinfTests/runtests.jl index eab0a74..9a596c3 100644 --- a/test/DAsemiinfTests/runtests.jl +++ b/test/DAsemiinfTests/runtests.jl @@ -16,4 +16,8 @@ using LightPropagation: fluence_DA_semiinf_CW, fluence_DA_semiinf_TD d = [0.000288659079311426, 2.896685181887841e-6, 5.474977667638328e-8, 1.3595622802163355e-9] @test fluence_DA_semiinf_TD(1.0:4.0, 1.0, 0.1, 10.0, n_med = 1.0, n_ext = 1.0, z = 0.0) ≈ d +# this test both the ForwardDiff package and the implementations of fluence and flux in semiinfite case +D = 1 / (3 * 10.0) # hard code diffusion coefficient here but could be different +@test D*ForwardDiff.derivative(z -> fluence_DA_semiinf_TD(1.0, 1.0, 0.1, 10.0, z = z), 0.0) ≈ flux_DA_semiinf_TD(1.0, 1.0, 0.1, 10.0, n_med = 1.0, n_ext = 1.0) + end # module \ No newline at end of file From 2e75c64d4bd0a7e9bc7fd38f69719d329ee9fe2a Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Thu, 12 Aug 2021 14:55:19 -0400 Subject: [PATCH 02/13] fix space --- .../Diffusion Approximation/Semi-infinite/DAsemiinf.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/forwardmodels/Diffusion Approximation/Semi-infinite/DAsemiinf.jl b/src/forwardmodels/Diffusion Approximation/Semi-infinite/DAsemiinf.jl index 5d91769..819a5bd 100644 --- a/src/forwardmodels/Diffusion Approximation/Semi-infinite/DAsemiinf.jl +++ b/src/forwardmodels/Diffusion Approximation/Semi-infinite/DAsemiinf.jl @@ -78,7 +78,7 @@ function fluence_DA_semiinf_TD(t, ρ, μa, μsp; n_ext = 1.0, n_med = 1.0, z = 0 elseif isa(t, AbstractVector) T = promote_type(eltype(ρ), eltype(μa), eltype(μsp), eltype(z)) ϕ = zeros(T, length(t)) - + @inbounds Threads.@threads for ind in eachindex(t) ϕ[ind] = _kernel_fluence_DA_semiinf_TD(D, ν, t[ind], μa, z, z0, ρ, zb) end @@ -174,7 +174,6 @@ function flux_DA_semiinf_TD(t, ρ, μa, μsp; n_ext = 1.0, n_med = 1.0) return Rt end - ##################################### # Steady-State Flux ##################################### From de5ab860dc97e5ed097ebb790a5bdd7d9842d81b Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Thu, 12 Aug 2021 15:50:01 -0400 Subject: [PATCH 03/13] fix allocation type and flux calls --- .../Diffusion Approximation/Slab/DAslab.jl | 124 +++++++++++------- 1 file changed, 73 insertions(+), 51 deletions(-) diff --git a/src/forwardmodels/Diffusion Approximation/Slab/DAslab.jl b/src/forwardmodels/Diffusion Approximation/Slab/DAslab.jl index e3a2926..f8c927c 100644 --- a/src/forwardmodels/Diffusion Approximation/Slab/DAslab.jl +++ b/src/forwardmodels/Diffusion Approximation/Slab/DAslab.jl @@ -90,7 +90,9 @@ function fluence_DA_slab_TD(t, ρ, μa, μsp; n_ext = 1.0, n_med = 1.0, s = 1.0, if isa(t, AbstractFloat) return _kernel_fluence_DA_slab_TD(D, ν, t, ρ, μa, s, zb, z0, z, xs) elseif isa(t, AbstractArray) - ϕ = zeros(eltype(D), length(t)) + T = promote_type(eltype(ρ), eltype(D), eltype(s), eltype(z)) + ϕ = zeros(T, length(t)) + @inbounds Threads.@threads for ind in eachindex(t) ϕ[ind] = _kernel_fluence_DA_slab_TD(D, ν, t[ind], ρ, μa, s, zb, z0, z, xs) end @@ -117,26 +119,44 @@ end end ##################################### -# Time-Domain Reflectance (flux) +# Time-Domain Flux ##################################### """ - refl_DA_slab_TD(t, ρ, μa, μsp, n_ext, n_med, s; xs = -15:15) + flux_DA_slab_TD(t, ρ, μa, μsp; n_ext = 1.0, n_med = 1.0, s = 1.0, z = 0.0, xs = 15) -Compute the time-domain reflectance (flux) from a slab geometry (x,y->inf, z-> finite). + Compute the time-domain flux (D*∂ϕ(t)/∂z @ z = 0 or z = s) from a slab geometry (x,y->inf, z-> finite). # Arguments - `t`: the time vector (ns). - `ρ`: the source detector separation (cm⁻¹) - `μa`: absorption coefficient (cm⁻¹) - `μsp`: reduced scattering coefficient (cm⁻¹) -- `ndet`: the boundary's index of refraction (air or detector) -- `nmed`: the sample medium's index of refraction +- `n_ext`: the boundary's index of refraction (air or detector) +- `n_med`: the sample medium's index of refraction - `s`: the thickness (z-depth) of the slab (cm) +- `z`: the z-depth coordinate (cm) +- `xs`: the number of sources to compute in the series # Examples -julia> `refl_DA_slab_TD(1.0, 1.0, 0.1, 10.0, 1.0, 1.0, 1.0)` +julia> `flux_DA_slab_TD(1.0, 1.0, 0.1, 10.0, n_ext = 1.0, n_med = 1.0, z = 0.0, s = 1.0)` +""" +function flux_DA_slab_TD(t, ρ, μa, μsp; n_ext = 1.0, n_med = 1.0, z = 0.0, s = 1.0, xs = 15) + if z == zero(eltype(z)) + return _refl_DA_slab_TD(t, ρ, μa, μsp, n_ext = n_ext, n_med = n_med, s = s, xs = xs) + elseif z == s + return _trans_DA_slab_TD(t, ρ, μa, μsp, n_ext = n_ext, n_med = n_med, s = s, xs = xs) + else + return D * ForwardDiff.derivative(dz -> fluence_DA_slab_TD(t, ρ, μa, μsp, n_ext = n_ext, n_med = n_med, s = s, z = dz, xs = xs), z) + end +end + """ -function refl_DA_slab_TD(t, ρ, μa, μsp, n_ext, n_med, s; xs = -15:15) + _refl_DA_slab_TD(t, ρ, μa, μsp, n_ext, n_med, s; xs = -15:15) + +Compute the time-domain flux or reflectance at the top surface from a slab geometry (x,y->inf, z-> finite). + +""" +function _refl_DA_slab_TD(t, ρ, μa, μsp; n_ext = 1.0, n_med = 1.0, s = 1.0, xs = 15) D = D_coeff(μsp, μa) A = A_coeff(n_med / n_ext) ν = ν_coeff(n_med) @@ -144,63 +164,65 @@ function refl_DA_slab_TD(t, ρ, μa, μsp, n_ext, n_med, s; xs = -15:15) z0 = z0_coeff(μsp) zb = zb_coeff(A, D) - Rt1 = Array{eltype(ρ)}(undef, length(t)) - Rt2 = zeros(eltype(ρ), length(t)) - - Threads.@threads for n in eachindex(t) - Rt1[n] = -exp(-(ρ^2 / (4 * D * ν * t[n])) - μa * ν * t[n]) / (2 * (4 * π * D * ν)^(3/2) * t[n]^(5/2)) - for m in xs - z3m = -2 * m * s - 4 * m * zb - z0 - z4m = -2 * m * s - (4m - 2) * zb + z0 - Rt2[n] += z3m * exp(-(z3m^2 / (4 * D * ν * t[n]))) - z4m * exp(-(z4m^2 / (4 * D * ν * t[n]))) - end - Rt1[n] *= Rt2[n] - end - - return Rt1 + if isa(t, AbstractFloat) + return _kernel_refl_DA_slab_TD(ρ, D, ν, t, μa, xs, s, zb, z0) + elseif isa(t, AbstractVector) + T = promote_type(eltype(ρ), eltype(D), eltype(s)) + Rt = zeros(T, length(t)) + + @inbounds Threads.@threads for ind in eachindex(t) + Rt[ind] = _kernel_refl_DA_slab_TD(ρ, D, ν, t[ind], μa, xs, s, zb, z0) + end + + return Rt + end +end +@inline function _kernel_refl_DA_slab_TD(ρ, D, ν, t, μa, xs, s, zb, z0) + Rt1 = -exp(-(ρ^2 / (4 * D * ν * t)) - μa * ν * t) / (2 * (4 * π * D * ν)^(3/2) * t^(5/2)) + Rt2 = zero(eltype(Rt1)) + for m in -xs:xs + z3m = -2 * m * s - 4 * m * zb - z0 + z4m = -2 * m * s - (4m - 2) * zb + z0 + Rt2 += z3m * exp(-(z3m^2 / (4 * D * ν * t))) - z4m * exp(-(z4m^2 / (4 * D * ν * t))) + end + + return Rt1 * Rt2 end -##################################### -# Time-Domain Transmittance (flux) -##################################### """ trans_DA_slab_TD(t, ρ, μa, μsp, n_ext, n_med, s; xs = -15:15) Compute the time-domain transmittance (flux) from a slab geometry (x,y -> inf, z -> finite) with Eqn. 39 from Contini 97. - -# Arguments -- `t`: the time vector (ns). -- `ρ`: the source detector separation (cm⁻¹) -- `μa`: absorption coefficient (cm⁻¹) -- `μsp`: reduced scattering coefficient (cm⁻¹) -- `n_det`: the boundary's index of refraction (air or detector) -- `n_med`: the sample medium's index of refraction -- `s`: the thickness (z-depth) of the slab (cm) - -# Examples -julia> `trans_DA_slab_TD(1.0, 1.0, 0.1, 10.0, 1.0, 1.0, 1.0)` """ -function trans_DA_slab_TD(t, ρ, μa, μsp, n_ext, n_med, s; xs = -15:15) +function _trans_DA_slab_TD(t, ρ, μa, μsp; n_ext = 1.0, n_med = 1.0, s = 1.0, xs = 15) D = D_coeff(μsp, μa) A = A_coeff(n_med / n_ext) ν = ν_coeff(n_med) z0 = z0_coeff(μsp) zb = zb_coeff(A, D) + + if isa(t, AbstractFloat) + return _kernel_trans_DA_slab_TD(ρ, D, ν, t, μa, xs, s, zb, z0) + elseif isa(t, AbstractVector) + T = promote_type(eltype(ρ), eltype(D), eltype(s)) + Rt = zeros(T, length(t)) - Rt1 = Array{eltype(ρ)}(undef, length(t)) - Rt2 = zeros(eltype(ρ), length(t)) - - Threads.@threads for n in eachindex(t) - Rt1[n] = exp(-(ρ^2 / (4 * D * ν * t[n])) - μa * ν * t[n]) / (2 * (4 * π * D * ν)^(3/2) * t[n]^(5/2)) - for m in xs - z1m = (1 - 2 * m) * s - 4 * m * zb - z0 - z2m = (1 - 2 * m) * s - (4 * m - 2) * zb + z0 + @inbounds Threads.@threads for ind in eachindex(t) + Rt[ind] = _kernel_trans_DA_slab_TD(ρ, D, ν, t[ind], μa, xs, s, zb, z0) + end - Rt2[n] += z1m * exp(-(z1m^2 / (4 * D * ν * t[n]))) - z2m * exp(-(z2m^2 / (4 * D * ν * t[n]))) - end - Rt1[n] *= Rt2[n] + return Rt end - - return Rt1 end +@inline function _kernel_trans_DA_slab_TD(ρ, D, ν, t, μa, xs, s, zb, z0) + Rt1 = exp(-(ρ^2 / (4 * D * ν * t)) - μa * ν * t) / (2 * (4 * π * D * ν)^(3/2) * t^(5/2)) + Rt2 = zero(eltype(Rt1)) + for m in -xs:xs + z1m = (1 - 2 * m) * s - 4 * m * zb - z0 + z2m = (1 - 2 * m) * s - (4 * m - 2) * zb + z0 + Rt2 += z1m * exp(-(z1m^2 / (4 * D * ν * t))) - z2m * exp(-(z2m^2 / (4 * D * ν * t))) + end + + return Rt1 * Rt2 +end \ No newline at end of file From 25cee55260cb240cf19aa9fa7f47ac515637f3e3 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Thu, 12 Aug 2021 16:03:05 -0400 Subject: [PATCH 04/13] add slab tests --- test/DAslabTests/runtests.jl | 27 +++++++++++++++------------ 1 file changed, 15 insertions(+), 12 deletions(-) diff --git a/test/DAslabTests/runtests.jl b/test/DAslabTests/runtests.jl index e97cc80..c384de8 100644 --- a/test/DAslabTests/runtests.jl +++ b/test/DAslabTests/runtests.jl @@ -8,33 +8,36 @@ using LightPropagation: refl_DA_slab_TD, trans_DA_slab_TD # Test slabs against semi-inf for sufficiently thick slabs +fluence_DA_slab_CW(ρ, μa, μsp; n_ext = 1.0, n_med = 1.0, s = 2.0, z = 0.0, xs = 10) # CW -@test fluence_DA_slab_CW(1.0, 0.1, 10.0, 1.0, 1.0, 20.0, 0.0) == 0.024035998288332954 -@test fluence_DA_slab_CW(2.4, 0.6, 18.0, 1.3, 1.2, 15.0, 1.0) == 1.77784095787581e-7 -@test fluence_DA_semiinf_CW(1.0, 0.1, 10.0, 1.0, 1.0, 0.0) ≈ fluence_DA_slab_CW(1.0, 0.1, 10.0, 1.0, 1.0, 20.0, 0.0) -@test fluence_DA_semiinf_CW(4.0, 0.01, 50.0, 1.0, 1.0, 0.0) ≈ fluence_DA_slab_CW(4.0, 0.01, 50.0, 1.0, 1.0, 20.0, 0.0) +@test fluence_DA_slab_CW(1.0, 0.1, 10.0, n_ext = 1.0, n_med = 1.0, s = 20.0, z = 0.0) == 0.024035998288332954 +@test fluence_DA_slab_CW(2.4, 0.6, 18.0, n_ext = 1.3, n_med = 1.2, s = 15.0, z = 1.0) == 1.77784095787581e-7 +@test fluence_DA_semiinf_CW(1.0, 0.1, 10.0, n_ext = 1.0, n_med = 1.0, z = 0.0) ≈ fluence_DA_slab_CW(1.0, 0.1, 10.0, n_ext = 1.0, n_med = 1.0, z = 0.0, s = 20.0) +@test fluence_DA_semiinf_CW(4.0, 0.01, 50.0, n_ext = 1.0, n_med = 1.0, z = 0.0) ≈ fluence_DA_slab_CW(4.0, 0.01, 50.0,n_ext = 1.0, n_med = 1.0, z = 0.0, s = 20.0) # TD t = 0.1:0.5:8 -@test fluence_DA_semiinf_TD(t, 1.0, 0.1, 10.0, 1.0, 1.0, 0.0) ≈ fluence_DA_slab_TD(t, 1.0, 0.1, 10.0, 1.0, 1.0, 40.0, 0.0) -@test fluence_DA_semiinf_TD(t, 1.0, 0.01, 70.0, 1.0, 1.0, 0.0) ≈ fluence_DA_slab_TD(t, 1.0, 0.01, 70.0, 1.0, 1.0, 40.0, 0.0) +@test fluence_DA_semiinf_TD(t, 1.0, 0.1, 10.0, n_ext = 1.0, n_med = 1.0, z = 0.0) ≈ fluence_DA_slab_TD(t, 1.0, 0.1, 10.0, n_ext = 1.0, n_med = 1.0, z = 0.0, s = 40.0) +@test fluence_DA_semiinf_TD(t, 1.0, 0.01, 70.0, n_ext = 1.0, n_med = 1.0, z = 0.0) ≈ fluence_DA_slab_TD(t, 1.0, 0.01, 70.0, n_ext = 1.0, n_med = 1.0, z = 0.0, s = 40.0) -@test fluence_DA_semiinf_TD(t, 3.0, 0.1, 10.0, 1.0, 1.0, 0.0) ≈ fluence_DA_slab_TD(t, 3.0, 0.1, 10.0, 1.0, 1.0, 40.0, 0.0) -@test fluence_DA_semiinf_TD(t, 3.0, 0.01, 70.0, 1.0, 1.0, 0.0) ≈ fluence_DA_slab_TD(t, 3.0, 0.01, 70.0, 1.0, 1.0, 40.0, 0.0) +@test fluence_DA_semiinf_TD(t, 3.0, 0.1, 10.0, n_ext = 1.0, n_med = 1.0, z = 0.0) ≈ fluence_DA_slab_TD(t, 3.0, 0.1, 10.0, n_ext = 1.0, n_med = 1.0, z = 0.0, s = 40.0) +@test fluence_DA_semiinf_TD(t, 3.0, 0.01, 70.0, n_ext = 1.0, n_med = 1.0, z = 0.0) ≈ fluence_DA_slab_TD(t, 3.0, 0.01, 70.0, n_ext = 1.0, n_med = 1.0, z = 0.0, s = 40.0) +@test fluence_DA_semiinf_TD(t, 3.0, 0.01, 70.0, n_ext = 1.4, n_med = 1.2, z = 0.1) ≈ fluence_DA_slab_TD(t, 3.0, 0.01, 70.0, n_ext = 1.4, n_med = 1.2, z = 0.1, s = 40.0) ## check flux # reflectance + D = 1 / (3 * 10.0) -@test (D * ForwardDiff.derivative(z -> fluence_DA_slab_TD(1.0, 1.0, 0.1, 10.0, 1.0, 1.0, 1.0, z), 0.0)) ≈ refl_DA_slab_TD(1.0, [0.1, 10.0], 1.0, 1.0, 1.0, 1.0)[1] +@test (D * ForwardDiff.derivative(dz -> fluence_DA_slab_TD(1.0, 1.0, 0.1, 10.0, n_ext = 1.0, n_med = 1.0, s = 1.0, z = dz), 0.0)) ≈ flux_DA_slab_TD(1.0, 1.0, 0.1, 10.0, n_ext = 1.0, n_med = 1.0, s = 1.0, z = 0.0) D = 1 / (3 * 20.0) -@test (D * ForwardDiff.derivative(z -> fluence_DA_slab_TD(1.5, 1.0, 0.4, 20.0, 1.0, 1.0, 1.0, z), 0.0)) ≈ refl_DA_slab_TD(1.5, [0.4, 20.0], 1.0, 1.0, 1.0, 1.0)[1] +@test (D * ForwardDiff.derivative(dz -> fluence_DA_slab_TD(1.0, 1.3, 0.3, 20.0, n_ext = 1.0, n_med = 1.0, s = 1.0, z = dz), 0.0)) ≈ flux_DA_slab_TD(1.0, 1.3, 0.3, 20.0, n_ext = 1.0, n_med = 1.0, s = 1.0, z = 0.0) # transmittance - D = 1 / (3 * 20.0) -@test (-D * ForwardDiff.derivative(z -> fluence_DA_slab_TD(1.5, 1.0, 0.4, 20.0, 1.0, 1.0, 1.0, z), 1.0)) ≈ trans_DA_slab_TD(1.5, [0.4, 20.0], 1.0, 1.0, 1.0, 1.0)[1] +@test (-D * ForwardDiff.derivative(dz -> fluence_DA_slab_TD(1.0, 1.3, 0.3, 20.0, n_ext = 1.0, n_med = 1.0, s = 1.0, z = dz), 1.0)) ≈ flux_DA_slab_TD(1.0, 1.3, 0.3, 20.0, n_ext = 1.0, n_med = 1.0, s = 1.0, z = 1.0) + end \ No newline at end of file From 6f8cad0000ba01430ffc1c29094fd6f034b1e7b5 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Thu, 12 Aug 2021 16:36:18 -0400 Subject: [PATCH 05/13] update slab tests --- .../Diffusion Approximation/Slab/DAslab.jl | 36 +++++++++++++++++++ test/DAslabTests/runtests.jl | 2 ++ 2 files changed, 38 insertions(+) diff --git a/src/forwardmodels/Diffusion Approximation/Slab/DAslab.jl b/src/forwardmodels/Diffusion Approximation/Slab/DAslab.jl index f8c927c..da7fb47 100644 --- a/src/forwardmodels/Diffusion Approximation/Slab/DAslab.jl +++ b/src/forwardmodels/Diffusion Approximation/Slab/DAslab.jl @@ -225,4 +225,40 @@ end end return Rt1 * Rt2 +end + +##################################### +# Steady-State Flux +##################################### +""" +fluence_DA_slab_CW(ρ, μa, μsp; n_ext = 1.0, n_med = 1.0, s = 2.0, z = 0.0, xs = 10) + + flux_DA_slab_CW(ρ, μa, μsp; n_ext = 1.0, n_med = 1.0, s = 2.0, z = 0.0, xs = 10) + +Compute the steady-state flux (D*∂ϕ(ρ)/∂z for z = [0, s]) from a slab geometry (x, y -> inf, z -> finite). + +# Arguments +- `ρ`: the source detector separation (cm⁻¹) +- `μa`: absorption coefficient (cm⁻¹) +- `μsp`: reduced scattering coefficient (cm⁻¹) +- `n_ext`: the boundary's index of refraction (air or detector) +- `n_med`: the sample medium's index of refraction +- `s`: the thickness (z-depth) of the slab (cm) +- `z`: the z-depth within slab (cm) +- `xs`: the number of sources to compute in the series + +# Examples +```jldoctest +julia> `flux_DA_slab_CW(1.0, 0.1, 10.0, n_ext = 1.0, n_med = 1.0, s = 2.0, z = 0.0)` +``` +""" +function flux_DA_slab_CW(ρ, μa, μsp; n_ext = 1.0, n_med = 1.0, s = 2.0, z = 0.0, xs = 10) + D = D_coeff(μsp, μa) + if z == zero(eltype(z)) + return D * ForwardDiff.derivative(dz -> fluence_DA_slab_CW(ρ, μa, μsp, n_med = n_med, n_ext = n_ext, s = s, z = dz, xs = xs), z) + elseif z == s + return -D * ForwardDiff.derivative(dz -> fluence_DA_slab_CW(ρ, μa, μsp, n_med = n_med, n_ext = n_ext, s = s, z = dz, xs = xs), z) + else + return D * ForwardDiff.derivative(dz -> fluence_DA_slab_TD(t, ρ, μa, μsp, n_ext = n_ext, n_med = n_med, s = s, z = dz, xs = xs), z) + end end \ No newline at end of file diff --git a/test/DAslabTests/runtests.jl b/test/DAslabTests/runtests.jl index c384de8..54881bd 100644 --- a/test/DAslabTests/runtests.jl +++ b/test/DAslabTests/runtests.jl @@ -39,5 +39,7 @@ D = 1 / (3 * 20.0) D = 1 / (3 * 20.0) @test (-D * ForwardDiff.derivative(dz -> fluence_DA_slab_TD(1.0, 1.3, 0.3, 20.0, n_ext = 1.0, n_med = 1.0, s = 1.0, z = dz), 1.0)) ≈ flux_DA_slab_TD(1.0, 1.3, 0.3, 20.0, n_ext = 1.0, n_med = 1.0, s = 1.0, z = 1.0) +## CW +@test flux_DA_slab_CW(1.0, 0.1, 10.0; n_ext = 1.0, n_med = 1.0, s = 40.0, z = 0.0, xs = 15) ≈ flux_DA_semiinf_CW(1.0, 0.1, 10.0; n_ext = 1.0, n_med = 1.0) end \ No newline at end of file From c763d049c80e29d55012a866b915e6e61fdf3f48 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Thu, 12 Aug 2021 16:36:55 -0400 Subject: [PATCH 06/13] update docs --- src/forwardmodels/Diffusion Approximation/Slab/DAslab.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/forwardmodels/Diffusion Approximation/Slab/DAslab.jl b/src/forwardmodels/Diffusion Approximation/Slab/DAslab.jl index da7fb47..4e7a72a 100644 --- a/src/forwardmodels/Diffusion Approximation/Slab/DAslab.jl +++ b/src/forwardmodels/Diffusion Approximation/Slab/DAslab.jl @@ -235,7 +235,7 @@ fluence_DA_slab_CW(ρ, μa, μsp; n_ext = 1.0, n_med = 1.0, s = 2.0, z = 0.0, xs flux_DA_slab_CW(ρ, μa, μsp; n_ext = 1.0, n_med = 1.0, s = 2.0, z = 0.0, xs = 10) -Compute the steady-state flux (D*∂ϕ(ρ)/∂z for z = [0, s]) from a slab geometry (x, y -> inf, z -> finite). +Compute the steady-state flux (D*∂ϕ(ρ)/∂z for z = 0 or s from a slab geometry (x, y -> inf, z -> finite). # Arguments - `ρ`: the source detector separation (cm⁻¹) From 4f32a4bc17455407e337a64fca02fc7316d217fa Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Thu, 12 Aug 2021 16:39:14 -0400 Subject: [PATCH 07/13] docs --- src/forwardmodels/Diffusion Approximation/Slab/DAslab.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/forwardmodels/Diffusion Approximation/Slab/DAslab.jl b/src/forwardmodels/Diffusion Approximation/Slab/DAslab.jl index 4e7a72a..b167118 100644 --- a/src/forwardmodels/Diffusion Approximation/Slab/DAslab.jl +++ b/src/forwardmodels/Diffusion Approximation/Slab/DAslab.jl @@ -124,7 +124,8 @@ end """ flux_DA_slab_TD(t, ρ, μa, μsp; n_ext = 1.0, n_med = 1.0, s = 1.0, z = 0.0, xs = 15) - Compute the time-domain flux (D*∂ϕ(t)/∂z @ z = 0 or z = s) from a slab geometry (x,y->inf, z-> finite). + Compute the time-domain flux, D*∂ϕ(t)/∂z for z = 0 and -D*∂ϕ(t)/∂z for z = s from a slab geometry (x, y -> inf, z -> finite). + If z != 0 or s will default to compute D*∂ϕ(t)/∂z for the z given. # Arguments - `t`: the time vector (ns). @@ -235,7 +236,8 @@ fluence_DA_slab_CW(ρ, μa, μsp; n_ext = 1.0, n_med = 1.0, s = 2.0, z = 0.0, xs flux_DA_slab_CW(ρ, μa, μsp; n_ext = 1.0, n_med = 1.0, s = 2.0, z = 0.0, xs = 10) -Compute the steady-state flux (D*∂ϕ(ρ)/∂z for z = 0 or s from a slab geometry (x, y -> inf, z -> finite). +Compute the steady-state flux, D*∂ϕ(ρ)/∂z for z = 0 and -D*∂ϕ(ρ)/∂z for z = s from a slab geometry (x, y -> inf, z -> finite). +If z != 0 or s will default to compute D*∂ϕ(ρ)/∂z for the z given. # Arguments - `ρ`: the source detector separation (cm⁻¹) From ac461808a9dc6a314ad866f10e1b760843b1d3f0 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Thu, 12 Aug 2021 17:40:56 -0400 Subject: [PATCH 08/13] update docs --- src/forwardmodels/Diffusion Approximation/Slab/DAslab.jl | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/src/forwardmodels/Diffusion Approximation/Slab/DAslab.jl b/src/forwardmodels/Diffusion Approximation/Slab/DAslab.jl index b167118..432987b 100644 --- a/src/forwardmodels/Diffusion Approximation/Slab/DAslab.jl +++ b/src/forwardmodels/Diffusion Approximation/Slab/DAslab.jl @@ -232,12 +232,10 @@ end # Steady-State Flux ##################################### """ -fluence_DA_slab_CW(ρ, μa, μsp; n_ext = 1.0, n_med = 1.0, s = 2.0, z = 0.0, xs = 10) - - flux_DA_slab_CW(ρ, μa, μsp; n_ext = 1.0, n_med = 1.0, s = 2.0, z = 0.0, xs = 10) + fluence_DA_slab_CW(ρ, μa, μsp; n_ext = 1.0, n_med = 1.0, s = 2.0, z = 0.0, xs = 10) -Compute the steady-state flux, D*∂ϕ(ρ)/∂z for z = 0 and -D*∂ϕ(ρ)/∂z for z = s from a slab geometry (x, y -> inf, z -> finite). -If z != 0 or s will default to compute D*∂ϕ(ρ)/∂z for the z given. + Compute the steady-state flux, D*∂ϕ(ρ)/∂z for z = 0 and -D*∂ϕ(ρ)/∂z for z = s from a slab geometry (x, y -> inf, z -> finite). + If z != 0 or s will default to compute D*∂ϕ(ρ)/∂z for the z given. # Arguments - `ρ`: the source detector separation (cm⁻¹) From f8792396df0735bebc7badee0d08ca83dd87261a Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Mon, 16 Aug 2021 14:20:23 -0400 Subject: [PATCH 09/13] add flux calculations in layered cylinders --- .../DAcylinder_layered.jl | 135 ++++++++++++++++-- 1 file changed, 120 insertions(+), 15 deletions(-) diff --git a/src/forwardmodels/Diffusion Approximation/DAcylinder_layered.jl b/src/forwardmodels/Diffusion Approximation/DAcylinder_layered.jl index 7f9a5fc..c24fa29 100644 --- a/src/forwardmodels/Diffusion Approximation/DAcylinder_layered.jl +++ b/src/forwardmodels/Diffusion Approximation/DAcylinder_layered.jl @@ -68,11 +68,56 @@ function fluence_DA_Nlay_cylinder_CW(data, besselroots) return fluence_DA_Nlay_cylinder_CW(data.ρ, data.μa, data.μsp, data.n_ext, data.n_med, data.l, data.a, data.z, besselroots) end +##################################### +# Steady-State Flux +##################################### +""" + flux_DA_Nlay_cylinder_CW(ρ, μa, μsp, n_ext, n_med, l, a, z, besselroots) + +Compute the steady-state flux using Fick's law D[1]*∂ϕ(ρ)/∂z for z = 0 (reflectance) and -D[end]*∂ϕ(ρ)/∂z for z = sum(l) (transmittance) in an N-layered cylinder. +Source is assumed to be located on the top middle of the cylinder. z must be equal to 0 or the total length sum(l) of cylinder. + +# Arguments +- `ρ`: source-detector separation in cylindrical coordinates (distance from middle z-axis of cylinder) (cm⁻¹) +- `μa`: absorption coefficient (cm⁻¹) +- `μsp`: reduced scattering coefficient (cm⁻¹) +- `n_ext`: the boundary's index of refraction (air or detector) +- `n_med`: the sample medium's index of refraction +- `l`: layer thicknesses (cm) +- `a`: cylinder radius (cm) +- `z`: source depth within cylinder: must be equal to 0 or sum(l) +- `besselroots`: roots of bessel function of first kind zero order J0(x) = 0 + +# Examples +julia> `flux_DA_Nlay_cylinder_CW(1.0, [0.1, 0.1], [10.0, 10.0], 1.0, [1.0, 1.0], [4.5, 4.5], 10.0, 0.0, besselroots)` +""" +function flux_DA_Nlay_cylinder_CW(ρ, μa, μsp, n_ext, n_med, l, a, z, besselroots) + @assert z == zero(eltype(z)) || z == sum(l) + D = D_coeff.(μsp, μa) + if z == zero(eltype(z)) + return D[1] * ForwardDiff.derivative(dz -> fluence_DA_Nlay_cylinder_CW(ρ, μa, μsp, n_ext, n_med, l, a, dz, besselroots), z) + elseif z == sum(l) + return -D[end] * ForwardDiff.derivative(dz -> fluence_DA_Nlay_cylinder_CW(ρ, μa, μsp, n_ext, n_med, l, a, dz, besselroots), z) + end +end +""" + flux_DA_Nlay_cylinder_CW(data, besselroots) + +Wrapper to flux_DA_Nlay_cylinder_CW(ρ, μa, μsp, n_ext, n_med, l, a, z, besselroots) with inputs given as a structure (data). + +# Examples +julia> data = Nlayer_cylinder(a = 10.0, l = [1.0, 1.0, 1.0, 2.0], z = 5.0) +julia> `flux_DA_Nlay_cylinder_CW(data, besselroots)` +""" +function flux_DA_Nlay_cylinder_CW(data, besselroots) + return flux_DA_Nlay_cylinder_CW(data.ρ, data.μa, data.μsp, data.n_ext, data.n_med, data.l, data.a, data.z, besselroots) +end + ##################################### # Time-Domain Fluence ##################################### """ - fluence_DA_Nlay_cylinder_TD(t, ρ, μa, μsp, n_ext, n_med, l, a, z, besselroots; ; N = 24, ILT = hyper_fixed) + fluence_DA_Nlay_cylinder_TD(t, ρ, μa, μsp, n_ext, n_med, l, a, z, besselroots; N = 24, ILT = hyper_fixed) Compute the time-domain fluence in an N-layered cylinder. Source is assumed to be located on the top middle of the cylinder. It is best to use `hyper_fixed` as the inverse laplace transform if t consists of many time points. Utilize `hyperbola` for a single time point. @@ -96,17 +141,11 @@ The lowest fluence value you can compute will be no less than the machine precis # Examples julia> `fluence_DA_Nlay_cylinder_TD(0.1:0.1:2.0, [0.1, 0.1], [10.0, 10.0], 1.0, [1.0, 1.0], [4.5, 4.5], 10.0, 0.0, besselroots)` """ -function fluence_DA_Nlay_cylinder_TD(t::AbstractFloat, ρ, μa, μsp, n_ext, n_med, l, a, z, besselroots; N = 24, ILT = hyperbola) - Rt = zero(eltype(t)) - Rt = ILT(s -> _fluence_DA_Nlay_cylinder_Laplace(ρ, μa, μsp, n_ext, n_med, l, a, z, s, besselroots), t, N = N) - - return Rt +function fluence_DA_Nlay_cylinder_TD(t::AbstractFloat, ρ, μa, μsp, n_ext, n_med, l, a, z, besselroots; N = 18, ILT = hyperbola) + return ILT(s -> _fluence_DA_Nlay_cylinder_Laplace(ρ, μa, μsp, n_ext, n_med, l, a, z, s, besselroots), t, N = N) end function fluence_DA_Nlay_cylinder_TD(t::AbstractArray, ρ, μa, μsp, n_ext, n_med, l, a, z, besselroots; N = 24, ILT = hyper_fixed) - Rt = zeros(eltype(μa), length(t)) - Rt = ILT(s -> _fluence_DA_Nlay_cylinder_Laplace(ρ, μa, μsp, n_ext, n_med, l, a, z, s, besselroots), t, N = N) - - return Rt + return ILT(s -> _fluence_DA_Nlay_cylinder_Laplace(ρ, μa, μsp, n_ext, n_med, l, a, z, s, besselroots), t, N = N) end """ fluence_DA_Nlay_cylinder_TD(data; besselroots, ILT = hyper_fixed) @@ -121,12 +160,68 @@ function fluence_DA_Nlay_cylinder_TD(t, data; bessels = besselroots, N = 24) return fluence_DA_Nlay_cylinder_TD(t, data.ρ, data.μa, data.μsp, data.n_ext, data.n_med, data.l, data.a, data.z, bessels, N = N) end -# this function is the base for the laplace transform in the time domain -function _fluence_DA_Nlay_cylinder_Laplace(ρ, μa, μsp, n_ext, n_med, l, a, z, s, besselroots) - ν = ν_coeff.(n_med) - μa_complex = μa .+ s ./ ν +##################################### +# Time-Domain Flux +##################################### +""" + flux_DA_Nlay_cylinder_TD(t, ρ, μa, μsp, n_ext, n_med, l, a, z, besselroots; N = 24) - return fluence_DA_Nlay_cylinder_CW(ρ, μa_complex, μsp, n_ext, n_med, l, a, z, besselroots) +Compute the time-domain flux using Fick's law D[1]*∂ϕ(ρ, t)/∂z for z = 0 (reflectance) and -D[end]*∂ϕ(ρ, t)/∂z for z = sum(l) (transmittance) in an N-layered cylinder. +Source is assumed to be located on the top middle of the cylinder. z must be equal to 0 or the total length sum(l) of cylinder. +Uses the hyperbola contour if t is an AbstractFloat and the fixed hyperbola contour if t is an AbstractVector. +∂ϕ(ρ, t)/∂z is calculated using forward mode auto-differentiation with ForwardDiff.jl + +# Arguments +- `t`: time point or vector (ns) +- `ρ`: source-detector separation in cylindrical coordinates (distance from middle z-axis of cylinder) (cm⁻¹) +- `μa`: absorption coefficient (cm⁻¹) +- `μsp`: reduced scattering coefficient (cm⁻¹) +- `n_ext`: the boundary's index of refraction (air or detector) +- `n_med`: the sample medium's index of refraction +- `l`: layer thicknesses (cm) +- `a`: cylinder radius (cm) +- `z`: source depth within cylinder +- `besselroots`: roots of bessel function of first kind zero order J0(x) = 0 +- `N`: number of Hankel-Laplace calculations + +# Examples +julia> `fluence_DA_Nlay_cylinder_TD(0.1:0.1:2.0, [0.1, 0.1], [10.0, 10.0], 1.0, [1.0, 1.0], [4.5, 4.5], 10.0, 0.0, besselroots)` +""" +function flux_DA_Nlay_cylinder_TD(t::AbstractVector, ρ, μa, μsp, n_ext, n_med, l, a, z, besselroots; N = 24) + @assert z == zero(eltype(z)) || z == sum(l) + D = D_coeff.(μsp, μa) + + ## need to know the type of z to preallocate vectors coorectly in hyper_fixed for autodiff + function _ILT(z::T) where {T} + return hyper_fixed(s -> _fluence_DA_Nlay_cylinder_Laplace(ρ, μa, μsp, n_ext, n_med, l, a, z, s, besselroots), t, N = N, T = T) + end + if z == zero(eltype(z)) + return D[1]*ForwardDiff.derivative(_ILT, z) + elseif z == sum(l) + return -D[end]*ForwardDiff.derivative(_ILT, z) + end +end +function flux_DA_Nlay_cylinder_TD(t::AbstractFloat, ρ, μa, μsp, n_ext, n_med, l, a, z, besselroots; N = 24) + @assert z == zero(eltype(z)) || z == sum(l) + D = D_coeff.(μsp, μa) + if z == zero(eltype(z)) + return D[1] * ForwardDiff.derivative(dz -> hyperbola(s -> _fluence_DA_Nlay_cylinder_Laplace(ρ, μa, μsp, n_ext, n_med, l, a, dz, s, besselroots), t, N = N), 0.0) + elseif z == sum(l) + return -D[end] * ForwardDiff.derivative(dz -> hyperbola(s -> _fluence_DA_Nlay_cylinder_Laplace(ρ, μa, μsp, n_ext, n_med, l, a, dz, s, besselroots), t, N = N), sum(l)) + end +end + +""" + flux_DA_Nlay_cylinder_TD(t, data; besselroots, N = 24) + +Wrapper to flux_DA_Nlay_cylinder_TD(t, ρ, μa, μsp, n_ext, n_med, l, a, z, besselroots; N = 24) with inputs given as a structure (data). + +# Examples +julia> data = Nlayer_cylinder(a = 10.0, l = [1.0, 1.0, 1.0, 2.0], z = 5.0) +julia> `flux_DA_Nlay_cylinder_TD(0.5:1.0:2.5, data, bessels = besselroots[1:600])` +""" +function flux_DA_Nlay_cylinder_TD(t, data; bessels = besselroots, N = 24) + return flux_DA_Nlay_cylinder_TD(t, data.ρ, data.μa, data.μsp, data.n_ext, data.n_med, data.l, data.a, data.z, bessels, N = N) end ##################################### @@ -171,6 +266,16 @@ function fluence_DA_Nlay_cylinder_FD(data, besselroots) return fluence_DA_Nlay_cylinder_FD(data.ρ, data.μa, data.μsp, data.n_ext, data.n_med, data.l, data.a, data.z, data.ω, besselroots) end +################################################################################ +# this function is the base for the laplace transform in the time-domain +################################################################################ +function _fluence_DA_Nlay_cylinder_Laplace(ρ, μa, μsp, n_ext, n_med, l, a, z, s, besselroots) + ν = ν_coeff.(n_med) + μa_complex = μa .+ s ./ ν + + return fluence_DA_Nlay_cylinder_CW(ρ, μa_complex, μsp, n_ext, n_med, l, a, z, besselroots) +end + ################################################################################ # _kernel_fluence_DA_Nlay_cylinder provides the main kernel of the program. # The inputs are similar to fluence_DA_Nlay_cylinder_CW. From e13dbb02fbe931ab50b52069c2651aa5b5848ec9 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Mon, 16 Aug 2021 14:26:56 -0400 Subject: [PATCH 10/13] export flux --- src/LightPropagation.jl | 28 +++++++++++++++++++++++++--- 1 file changed, 25 insertions(+), 3 deletions(-) diff --git a/src/LightPropagation.jl b/src/LightPropagation.jl index c984661..f46e04a 100644 --- a/src/LightPropagation.jl +++ b/src/LightPropagation.jl @@ -8,29 +8,51 @@ using FFTW using Parameters using SpecialFunctions using JLD +using ForwardDiff export TPSF_DA_semiinf_refl export TPSF_DA_slab_refl export TPSF_DA_paralpip_refl +### Diffusion Theory + +# Infinite export fluence_DA_inf_CW export fluence_DA_inf_TD + +# Semi-infinite export fluence_DA_semiinf_TD export fluence_DA_semiinf_CW + +export flux_DA_semiinf_CW +export flux_DA_semiinf_TD + +# Slab export fluence_DA_slab_CW export fluence_DA_slab_TD + +export flux_DA_slab_CW +export flux_DA_slab_TD + +# Parallelepiped export fluence_DA_paralpip_TD export fluence_DA_paralpip_CW -export Nlayer_cylinder +# Layered cylinder export fluence_DA_Nlay_cylinder_CW export fluence_DA_Nlay_cylinder_TD + +export flux_DA_Nlay_cylinder_CW +export flux_DA_Nlay_cylinder_TD + +# Structures +export Nlayer_cylinder + +# constants export besselroots export getfit - export load_asc_data - export fitDTOF, DTOF_fitparams include("forwardmodels/Diffusion Approximation/diffusionparameters.jl") From 5a7baee9c55f1d14312bf62891d26555cd5345d4 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Mon, 16 Aug 2021 14:27:56 -0400 Subject: [PATCH 11/13] update hyperbola contours for autodiff --- .../Diffusion Approximation/transforms.jl | 26 +++++++------------ 1 file changed, 9 insertions(+), 17 deletions(-) diff --git a/src/forwardmodels/Diffusion Approximation/transforms.jl b/src/forwardmodels/Diffusion Approximation/transforms.jl index b46b6d4..bcfa76e 100644 --- a/src/forwardmodels/Diffusion Approximation/transforms.jl +++ b/src/forwardmodels/Diffusion Approximation/transforms.jl @@ -68,10 +68,9 @@ function hyperbola(f::Function, t::AbstractFloat; N = 16) return imag(a) * h / pi end -# loop through an entire time array with multithreads -function hyperbola(f::Function, t::AbstractVector; N = 16) - out = similar(t) - Threads.@threads for ind in eachindex(t) +function hyperbola(f::Function, t::AbstractArray; N = 16) + out = fill(f(t[1]), length(t)) + Threads.@threads for ind in 2:length(t) out[ind] = hyperbola(f, t[ind], N = N) end return out @@ -93,10 +92,10 @@ s_fixed(θ, μ; ϕ = 1.09) = μ + im * μ * sinh(θ + im * ϕ) ds_fixed(θ, μ; ϕ = 1.09) = im * μ * cosh(θ + im * ϕ) # compute the function values over the fixed contour nodes -function fixed_sk(f::Function, N, t::AbstractVector) +function fixed_sk(f::Function, N, t::AbstractVector, T) μ, h = hyper_coef(N, t) - a = zeros(Complex{eltype(h)}, Int(N)) - sk = similar(a) + a = zeros(Complex{eltype(T)}, N) + sk = zeros(Complex{eltype(t)}, N) Threads.@threads for k in 0:Int(N)-1 sk[k+1] = s_fixed((k + 1/2) * h, μ) dsk = ds_fixed((k + 1/2) * h, μ) @@ -113,18 +112,11 @@ function hyper_fixed_points(a, sk, h, t::AbstractFloat) return imag(b) * h / π end - - -function hyper_fixed(f::Function, t::AbstractVector; N = 24) - N = convert(eltype(t), N) - a, sk, h = fixed_sk(f, N, t) - out = zeros(eltype(h), length(t)) +function hyper_fixed(f::Function, t::AbstractVector; N = 24, T = eltype(t)) + a, sk, h = fixed_sk(f, N, t, T) + out = zeros(T, length(t)) Threads.@threads for ind in eachindex(t) out[ind] = hyper_fixed_points(a, sk, h, t[ind]) end return out end - -# run as -# LT_hyper_fixed(s -> fluence_DA_inf_FD(1.0, [0.1, 10.0], s),28, t) - From 19fd102006409c8dc7e087040d2e27f6df30ea67 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Mon, 16 Aug 2021 14:34:08 -0400 Subject: [PATCH 12/13] add flux tests for layered --- .../DA_NlayerTests/DA_Nlayer_functionTests.jl | 25 +++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/test/DA_NlayerTests/DA_Nlayer_functionTests.jl b/test/DA_NlayerTests/DA_Nlayer_functionTests.jl index c6c1f64..cd1eda0 100644 --- a/test/DA_NlayerTests/DA_Nlayer_functionTests.jl +++ b/test/DA_NlayerTests/DA_Nlayer_functionTests.jl @@ -145,4 +145,29 @@ a96 = fluence_DA_Nlay_cylinder_TD(t, cylinder_data, bessels = besselroots[1:600] @test a72 ≈ slab @test a96 ≈ slab + +### Flux tests + +# CW +cylinder_data = Nlayer_cylinder(ρ = 1.0, μsp = [10.0, 10.0], μa = [0.1, 0.1], l = [10.0, 10.0], n_med = [1.0, 1.0], n_ext = 1.0, a = 40.0) +@test flux_DA_Nlay_cylinder_CW(cylinder_data, besselroots[1:10000]) ≈ flux_DA_semiinf_CW(1.0, 0.1, 10.0; n_ext = 1.0, n_med = 1.0) + +cylinder_data = Nlayer_cylinder(ρ = 1.0, μsp = [10.0, 10.0], μa = [0.1, 0.1], l = [2.0, 2.0], z = 4.0, n_med = [1.0, 1.0], n_ext = 1.0, a = 40.0) +@test flux_DA_Nlay_cylinder_CW(cylinder_data, besselroots[1:10000]) ≈ flux_DA_slab_CW(1.0, 0.1, 10.0; n_ext = 1.0, n_med = 1.0, s = 4.0, z = 4.0, xs = 15) + +# TD +t = 1.0 +cylinder_data = Nlayer_cylinder(ρ = 1.0, μsp = [10.0, 10.0], μa = [0.1, 0.1], l = [10.0, 10.0], n_med = [1.0, 1.0], n_ext = 1.0, a = 40.0) +@test flux_DA_Nlay_cylinder_TD(t, cylinder_data, bessels = besselroots[1:10000]) ≈ flux_DA_semiinf_TD(t, 1.0, 0.1, 10.0; n_ext = 1.0, n_med = 1.0) + +cylinder_data = Nlayer_cylinder(ρ = 1.0, μsp = [10.0, 10.0], μa = [0.1, 0.1], l = [2.0, 2.0], z = 4.0, n_med = [1.0, 1.0], n_ext = 1.0, a = 40.0) +@test flux_DA_Nlay_cylinder_TD(t, cylinder_data, bessels = besselroots[1:10000]) ≈ flux_DA_slab_TD(t, 1.0, 0.1, 10.0; n_ext = 1.0, n_med = 1.0, s = 4.0, z = 4.0, xs = 15) + +t = 0.5:0.1:1.5 +cylinder_data = Nlayer_cylinder(ρ = 1.0, μsp = [10.0, 10.0], μa = [0.1, 0.1], l = [10.0, 10.0], n_med = [1.0, 1.0], n_ext = 1.0, a = 40.0) +@test flux_DA_Nlay_cylinder_TD(t, cylinder_data, bessels = besselroots[1:10000]) ≈ flux_DA_semiinf_TD(t, 1.0, 0.1, 10.0; n_ext = 1.0, n_med = 1.0) + +cylinder_data = Nlayer_cylinder(ρ = 1.0, μsp = [10.0, 10.0], μa = [0.1, 0.1], l = [2.0, 2.0], z = 4.0, n_med = [1.0, 1.0], n_ext = 1.0, a = 40.0) +@test flux_DA_Nlay_cylinder_TD(t, cylinder_data, bessels = besselroots[1:10000]) ≈ flux_DA_slab_TD(t, 1.0, 0.1, 10.0; n_ext = 1.0, n_med = 1.0, s = 4.0, z = 4.0, xs = 15) + end # module From daea97a1cbf0a77300cb9cf76b402b9e579cf151 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Mon, 16 Aug 2021 14:49:57 -0400 Subject: [PATCH 13/13] update docs with flux --- docs/src/API.md | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/docs/src/API.md b/docs/src/API.md index 4ec8822..029dcd9 100644 --- a/docs/src/API.md +++ b/docs/src/API.md @@ -10,16 +10,29 @@ fluence_DA_inf_TD ```@docs fluence_DA_semiinf_CW fluence_DA_semiinf_TD +flux_DA_semiinf_CW +flux_DA_semiinf_TD ``` ## Slab ```@docs fluence_DA_slab_CW fluence_DA_slab_TD +flux_DA_slab_CW +flux_DA_slab_TD ``` ## Parallelepiped ```@docs fluence_DA_paralpip_CW fluence_DA_paralpip_TD -``` \ No newline at end of file +``` + +## Layered Cylinder +```@docs +fluence_DA_Nlay_cylinder_CW +fluence_DA_Nlay_cylinder_TD + +flux_DA_Nlay_cylinder_CW +flux_DA_Nlay_cylinder_TD +```