Smarter computation of spherical Bessel functions for real arguments
Calculating the Mie scattering coefficients requires the computation of spherical Bessel functions j_l(x) and h_^1_l(x) (or, equivalently, Riccatti-Bessel functions) for real arguments. Right now, the Fortran extensions for scattering also need these special functions for the radial dependence of the scattered field, and obtain them via A. Ross Barnett's SBESJY.F.
scipy.special.
http://
These use upwards recursion for y_l and downwards recursion for j_l, which is good numerical practice. However, the Barnett codes use a more sophisticated approach based on Lentz's continued fraction method. Numerical experimentation shows that there exist a range of x (>5000) where the scipy.special routines fail by giving infinity or nan, while the Barnett routines still work. Should we remove dependence on scipy.special and replace with Barnett?
N.B: We handle the considerably more dangerous case of complex arguments by avoiding direct computation and using logarithmic derivatives.
[VNM] I would suggest we first see if this can be fixed upstream, by emailing Jin or seeing if scipy would be willing to use the Barnett version. Jerome, do you want to try to contact Jin about this?
Blueprint information
- Status:
- Not started
- Approver:
- None
- Priority:
- Undefined
- Drafter:
- Jerome Fung
- Direction:
- Needs approval
- Assignee:
- None
- Definition:
- Discussion
- Series goal:
- None
- Implementation:
-
Unknown
- Milestone target:
- None
- Started by
- Completed by
Related branches
Related bugs
Sprints
Whiteboard
The main issue we need to determine is the degree to which we will require Fortran extensions. Currently, it is at least possible in principle (the way the build system and packages work would make this a kluge) to calculate Mie coefficients and radiometric quantities without invoking any Fortran code. Given that Fortran is required to calculate any position-dependent scattered quantity, do we just decide to require Fortran for all scattering? Or is it important to keep as much as we can in pure Python?
