Strange steppy behavior of luminance after sunset

[10110111]

I've noticed an interesting pattern of change of luminance after sunset. Basically, as a scattering order comes to its dusk, luminance of the sky decreases fast, and only then does the next order take over and slow down the decrease. See the following plot:

It was done using 32 texels in r, 128 in μ, 2048 in μₛ and 4 in ν, and with max_sun_zenith_angle=π. Changing the tradeoff between resolution in μₛ and in ν doesn't influence much the result.

This effect can also be easily seen for general illuminance if you try plotting the irradiance.dat data. There the same steps are present.

I wonder, what could explain this? Is this result even physical? Or is this some known artifact of the model used here?

[gzotti]

That's interesting. You know the usual limits of Civil, Nautical and Astronomical twilight where the sun is at -6, -12 and -18 degrees. This plot could indicate that astronomical twilight ends approximately where second-order scattering falls off. Now it would be interesting to correlate this with measurements that also record things like natural airglow. If I see the log units in the graph, I wonder at which point adding more orders of scattering really does no longer contribute to improvements in simulation quality. E.g., if airglow lies at 10^-12 arb.units, there is no point in computing more than third-order. If if is at 10^-10, second-order scattering is probably enough. Not sure about conditions close to the horizon above the sun when it is still at -20 degrees, though, maybe instrumentation can detect something, but you would need to separate this out from e.g. Zodiacal light.

[ebruneton]

The code is "reasonably" physically accurate (see the comparison with libRadtran in https://arxiv.org/pdf/1612.04336.pdf), but has not been validated for such low Sun elevations. So I'm not sure the result is physical. You could try with libRadtran to have a comparison.

[10110111]

OK, I've tried with libRadtran, and it does exhibit the first step. The shape of the curve differs though. Couldn't test other steps, because, as I understood, the only solver in libRadtran capable of handling negative solar elevations is the Monte-Carlo solver—MYSTIC—and it converges very slowly for solar elevations less than about -10° (e.g. can give several zero radiances in a row for 10⁸ photons at elevation of -20°, while expected value is non-zero).

What's interesting, I've found that much of this steppy behavior (second and higher-order "twilight") is unphysical, and it stems from inadequate angular resolution of integration in ComputeScatteringDensity. Namely, increasing SAMPLE_COUNT in this function from 16 to 32, 64 and 96 gives the following results (pure Rayleigh atmosphere, no ozone, black ground, radiance at 680 nm only): Versions with 64 and 96 samples appear to give almost the same curves, so the curve of the former isn't visible behind the one of the latter. To make the better-sampled version easier to see, here's it alone. Note how the second step is almost invisible now:

And here's how the first step compares to libRadtran (elevations α>-10° were sampled randomly, others were hand-picked): A magnified view of the step: Not sure why libRadtran gives higher radiance for -7.5°≤α≤-5.5°.

[gzotti]

I don't know libRadtran, does it account for airglow? (recombination would not be found in pure scattering computations.)

[10110111]

I don't know libRadtran, does it account for airglow?

No, it doesn't, at least not by default. But anyway, the only region where it definitely disagrees with the best-precision version of this GPU model is -7.5°≤α≤-5.5°., where airglow should be less important than in the lower solar elevation cases.

[gzotti]

Or any other difference around Ozone distribution? Or is there a non-logarithmic density falloff? Sorry, just guessing from the outside of the blackbox.

[ebruneton]

Thanks, this is very interesting. The solar elevation at which the point on the "top of the atmosphere" at the zenith no longer receives any direct Sun light is asin(6360/6420)-90 = -7.84°, which seems to correspond to the step's inflection point (could be confirmed by playing with the planet and atmosphere radii). So this step could be "physical" after all, for an atmosphere that would abruptly end at some radius (which is not the case in practice).

[gzotti]

Superb, of course! This explains the "bend". Do you have any insight in the excess brightness of libRadtrans vs. scattering?

[10110111]

OK, apparently, 100 layers of atmosphere (value taken from clear-sky-models comparison) were not enough to simulate the exponential tail of atmospheric density. So libRadtran showed excessive radiance of this tail due to its (linear?) interpolation of density which always overestimates the exponent. After I split the density into 1000 layers, libRadtran now shows reasonable agreement with this GPU model:

[gzotti]

Very good, that is really close now! Does that mean, GPU model with 64 samples (line overlapped with 96 samples) provides all we need? Is this with "pure" Rayleigh scattering, or how does Mie scattering for various pollutants come into play?

[10110111]

Does that mean, GPU model with 64 samples (line overlapped with 96 samples) provides all we need?

I guess yes.

Is this with "pure" Rayleigh scattering, or how does Mie scattering for various pollutants come into play?

This was pure Rayleigh without any Mie or absorption, but I don't think we'll get any surprises if we enable aerosols & ozone. Best is to check, of course. Might give it a shot at some point.