Deal with attenuation below 2500 meters

[clark2668]

Three of AraSim's ice attenuation length functions extrapolate into depths where the underlying fits aren't reliable, and the code that uses them divides by L without checking whether the value makes sense. This PR addresses both. Per D. Besson, the attenuation model isn't trustworthy beyond ~2.5 km depth, and downstream code in Report computes exp(-d/L) and 1/L without guarding against division by zero issues, so bad attenuation lengths can show up as NaNs or unphysical signal estimates in the output.

This has changes to both the ice model and report.

In Icemodel:

• Added two public static constexpr class constants:
  • BEDROCK_DEPTH_M = 2850. (South Pole nominal ice column thickness)
  • MIN_PHYSICAL_ATTEN_M = 1. (attenuation lengths below this are treated as unphysical)
• GetFreqDepIceAttenuLength, GetARAIceAttenuLength, and EffectiveAttenuationLength now return 0 for:
  • negative depths (preserving the existing VERBOSE_MODE warning)
  • depths > 2500 m (the unreliable regime)
  • depths > bedrock (fail-safe in case the 2.5 km cutoff is ever relaxed)
• temperature() has a new optional bool strict = false argument:
  • With strict=false (the default, which keeps the API ~the ame), the cubic fit is evaluated and then capped at 0 via std::min(temp, 0.) so it can't return positive (above-freezing) temperatures at depth.
  • With strict=true, out-of-ice depths (z < 0 or z > BEDROCK_DEPTH_M) throw a runtime error. Generally people should use strict mode, but this way they can plot if they want to.
• Removed the EffectiveAttenuationLength(const Position&, const int&) const overload. it was unused.

In report, the new depth cutoffs mean the attenuation functions can now return 0. This runs a risk of a nan in exp(-dl / 0). So I put some safe guards in. If the returned value is less than IceModel::MIN_PHYSICAL_ATTEN_M, IceAttenFactor is set to 0.

To demonstrate this works, I plotted attenuation length and temperature vs depth. The old one is first, the new one is second. Note in the new one that the attention length floors to 0 meters, and the temp floors to 0 deg C.

I also added helper scripts for making the plots in examples/atten_length.

Addresses #224.

[clark2668]

I should also note that this change affects the RNG state, and therefore the unit tests.

In these files, you’ll see the first 10 events before and after my change, with info looking like this.

Event 0 : depth 2254.85 : RNG seed 3343220404
     Delta View Ang Str 0, Ant 0, Idx 0 : -14.202060 deg
     Delta View Ang Str 0, Ant 0, Idx 1 : -16.942841 deg
     Delta View Ang Str 0, Ant 1, Idx 0 : -14.220758 deg
     Delta View Ang Str 0, Ant 1, Idx 1 : -16.929951 deg
....
     Ants with nonzero signal 0

That’s event number, the depth of the interaction, and one of the RNG state variables at that point in the code. I also print out the angle from the cherenkov angle, and how many antennas passing the noise check that AraSim sees. What you can see is that old and new agree through event 7, but change going into event 8. And the thing that’s different is that “old” thinks ev 7 has 16 signals, but “new” thinks it has none. And that’s because ev 7 is the first event in the new “no no” region that also was close enough to the cherenkov angle to pass the non-zero signal check and get the noise builder to run. And that’s where the RNGs diverge. Because where ev 7 used to trigger the noise construction in the “old”, in the “new” that doesn’t happen.

arasim_old.txt
arasim_new.txt

[clark2668]

(Updated weights will be calculated this weekend)