pkienzle
commented
1 year ago To evaluate the Hankel transform we are using $$G(ξ) = \int_0^\infty I(q) J_0(qξ)q dq = \int0^{q{min}} I(q) J0(qξ)q dq + \int{q{min}}^{q{max}} I(q) J0(qξ)q dq + \int{q_{max}}^\infty I(q) J0(qξ)q dq = ε{min} + ε{mid} + ε{max}$$ where we approximate $$ε_{mid} ≈ \sum I(q_k) J_0(q_k ξ) q_k$$
The goal is to find $q{min}$ and $q{max}$ so that we achieve a desired relative error $ε_{max}/G(ξ) < Δ$ for some small constant $Δ$. What follows are some musings on the topic.
Consider a power law $I(q) \propto q^{-4}$ as the limiting case for $I(q)$. Since $ε{min}$ diverges, we need to introduce a correlation length cutoff $R{max}$ where $I(q) = q{min}^{-4}$ is constant for $q < q{min}$. The cutoff will be determined by the coherence length of the neutron, and may be above the largest correlation length probed in the measurement. Analytically this gives $$ε{min}(ξ) = I(q{min}) J1(q{min} ξ)q_{min}/ξ $$
For the $ε_{max}$ term we can use the approximation $$J0(z) ≈ \sqrt\frac{2}{πz}\cos(z-π/4)$$ for $z = qξ > 1$. Ignoring the oscillation in the cosine we can get a weak upper bound $$ε{max} < \int{q{max}}^\infty I(q)\sqrt\frac{2}{πqξ} q dq = \sqrt\frac{2}{πξ} \frac{q{max}^{-2.5}}{2.5}$$ For $q{max} = 10$ / nm and $ξ{min} = 1$ nm this gives $ε{max} ≈ 0.0010$, whereas the value from the Hankel integral ≈-0.00012.
Even if we had a good bound for $ε{max}$ we can't use it to estimate error since we don't have a bound on $ε{mid}$. Let's assume that oscillations in $I(q)$ are small relative to $I(q)$ and we might be able to get away with using the envelope function for $J0$, which gives $$ε{mid}+ε{max} \sim \int{q{min}}^\infty I(q)\sqrt\frac{2}{πqξ} q dq = \sqrt\frac{2}{πξ} \frac{q{min}^{-2.5}}{2.5}$$ The hope is that the overestimate in this term will cancel the overestimate in the $ε{max}$ term. Then we can set $ε{max} / (ε{mid} + ε{max}) < Δ$. Solving this gives $$q{max} = q{min} Δ^{-1/2.5}$$ For $Δ=10^{-5}$ this gives $q{max} = 100 q{min}$. Testing for a variety $q_{min}$ and $ξ$ in the SANS/USANS range using the Hankel integral rather than the approximation the relative error is sometimes 10x too high.
This scheme only works when $q{min}$ is at the start of the power law region. For the sphere model this is at $q{min} = 1/R$. Given an $R{max}$ parameter, then, we would need $q{min} = 0.1 / R{max}$ and $q{max} = 100 / R_{max}$ for better than 0.1% accuracy. Checking $\int_0^y J_0(ξ q) (3 j_1(q R)/(q R))^2 q dq$ for $y = 10^n/R$ the range from $n=-1$ to $n=2$ gives good accuracy. I don't expect structure factor will have a big influence: the $J_0$ is oscillating so rapidly that the integral will be dominated by the $J_0$ envelope, especially at high $ξ$. Lower values of $ξ$ require a larger $q$ range to achieve the same accuracy.
Our $q_k$ are sampled far to sparsely to hit every oscillation in $J_0$. We should compute a Hankel integrals to high precision for a few select (model, ξ) sets to verify that our sampling is high enough.
Looking at the calculation of sesans $G(ξ) = \frac{1}{2\pi} \int I(q) J_0(qξ) q dq$, the value of G is going to be driven primarily by the rapid decrease in I(q), with higher q values contributing minimally to the integral.
The proposed qmin, qmax range for single point ξ in #563 will not work. Instead q min should be based on the optional parameter Rmax, which defaults to ξ max if it is not provided. q max can be set to a number of decades above q min, assuming scattering intensity follows a power law, or perhaps from an additional parameter Rmin whose default value comes from ξ min, with a minimal number of decades in q max/q min.
In principle the calculation should be independent of the particular ξ points, with a single point calculation giving the same result as that point in a vector of ξ values. That is, we should only use Rmin and Rmax to determine q, not the ξ density.
Regarding structure factors, there will be correlations beyond maximum particle dimension but these show up as structure factor peaks scaling the I(q) values. Again, since high q values contribute minimally to the integral, we may not need to increase Rmax to capture the correlation structure in ξ above Rmax. This is another reason the Rmax parameter should not be ignored.