Merge pull request #4 from heltonmc/flux

Calculate flux of N-layered solutions using automatic differentiation

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
-```
+```
+
+## Layered Cylinder
+```@docs
+fluence_DA_Nlay_cylinder_CW
+fluence_DA_Nlay_cylinder_TD
+
+flux_DA_Nlay_cylinder_CW
+flux_DA_Nlay_cylinder_TD
+```

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")

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.

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,51 +121,80 @@ 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)
 
     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