Calculate flux of N-layered solutions using automatic differentiation

[heltonmc]

This creates and edits several new features:

1. (more straightforward) Allows for calculation of flux using fick's law for the steady-state case in the layered solution. This can be done fairly easily with ForwardDiff by a simple call like...
   ForwardDiff.derivative(dz -> fluence_DA_Nlay_cylinder_CW(ρ, μa, μsp, n_ext, n_med, l, a, dz, besselroots), z). This allows the fluence and flux to be calculated in the same computation time.
2. (more complicated) In the time domain an approach like for the steady-state could not be done quite as simply except for the case when t is just an AbstractFloat. However, if we want to calculate for many time points using hyper_fixed there were several errors to fix. The primary one being is that ForwardDiff uses Dual numbers that are embedded within the function arguments that are passed to hyper_fixed. There are several spots where hyper_fixed allocates before evaluated the function so did not know the type. There were a couple of options that were tried here. One was just to evaluate the anonymous function passed to hyper_fixed first for a single value then allocate based on those types. However, we want to be able to thread these solutions because the single evaluation is costly. Instead, I updated the hyper_fixed functions to create an optional input argument being the eltype to preallocate the vectors in hyper_fixed with. This is done with creating an intermediate function:
   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 where now the type parameter of z which is what we are differentiating with respect to is known at compile time an d fed to hyper_fixed.
3. Several other edits to the semi-inf and slab cases for flux were also exported.