Add to the report types (e.g. GEN, REF, CEL, SCT, TMO, etc) a new reportable event, NAN, and report and kill phonon in the following situations:

• Phonons asked to move a PathLength or TravelTime of NaN
• Phonon accrues negative travel time or path length
• Phonon has multiple consecutive "moves" of zero path length.

Note in latter case, it must be multiple consecutive zero-length moves, as reflection events, for example, produce legitimate "moves" where only the direction, not the location, are changed. (There's one known bug (#5) where a phonon gets stuck reflecting back and forth ad infinitum without actually moving anywhere...)

Additionally, add a counter so that at the end of simulation a report can be given of how many of these events occurred.

Related issues: #7 , #5

Dear Christopher,

Thank you for sharing your code with the community!
I had a question about how the user could define the desired scattering operators (differential scattering cross-sections). The analytical expressions for these operators being different based on the symmetry class of the background (average) and foreground (fluctuation) media, it would be interesting to be able to add that option to the code.

Best regards,
Shahram

Specifying model grids in code and compiling them in as hard-coded options is, of course, needlessly tedious. Write a basic/simple file ingestor that can take these grids in from an ascii file.

Compiling with 32-bit floats instead of 64-bit floats offers significant memory savings, and perhaps some speedup as well. However, there is a stuck-phonon issue that happens ONLY when using 32-bit floats and has not been observed when using 64-bit floats.

Track this down and fix it.

It is already possible to input Vs=0 into the model grid, and basic ray-tracing runs fine. However, the following issues exist with the scattering and Q dynamics:

• Specifying a finite, non-zero Qs via a QmQk() constructor, when Vs is zero, produces infinite Qp, and does not alert user that an illogical value was input, and a likely unintended Qp value was inferred. In this case, the user should inout Qp explicitly, rather than rely on automatic computation.

• Vs==0 results in scattering parameters el and gam0 being infinite, and MFPs for both P and S end up being NaNs. Not actually sure what effect this has during runtime. Probably effectively turns off scattering in these regions. Not a problem for S waves, since they won't be there, but for the P waves this is probably an unintended result.

Phonons trending exactly vertically (towards or away from model center) in SphereShell cells with radial quadratic velocities have infinite arc radius (they follow straight lines). This is actually handled by falling back to the linear handlers used for uniform-velocity SphereShells, which gives correct answers for distances, and then a special case solver is used for travel times. The math seems right, and gives correct answer for roughly 49 out of every 50 phonons simulated, but about 1 in 50 picks up a NaN somewhere and "gets stuck". Haven't figured out where this is happening yet except I think it has it's origins in SphereShell::GetRayArc_RD2(), where a check for verticality is (sometimes) failing to trigger, and and the intended sentinel value for Ray.Center isn't getting assigned. But again, no idea why.

The good news is these perfectly vertical phonons are exceedingly rare in practice, and are only statistically likely to produce a NaN. (In fact, I've only been able to test them by force injecting them — haven't actually seen one in the wild.) Still, I'm not quite comfortable with the chance that one of these could scuttle a simulation that is potentially hours or days into computing.

So... need to track down and squash.

Make the code multi-threaded.

The spherical shells work well with radial quadratic gradients — the ray paths are circular arc segments. But the solutions are computed under that assumption that the coefficient on r^2 is negative, and the constant term positive. (Meaning: velocities increase downwards.)

In some earth models, there a regions where velocity decreases as you go downwards. Although these can be approximated with stair-steps, it would be nice if the ray paths for downward-decreasing gradients could be handled directly, allowing much more modeling freedom.