Inconsistencies between radiance calculations, unexpected impact of changing ε

[gerritholl]

Differences between radiances calculated directly and through the symbolic measurement equation are too large when either emissivity or non-linearity are switched off individually. With a codebase that includes a15aeed and 9e67af7 (such that #336 should not be the cause anymore), with a breakpoint at

https://github.com/FIDUCEO/FCDR_HIRS/blob/9e67af790da5a9f541ad3b7ad69702d77a33509a/FCDR_HIRS/fcdr.py#L1722

I then prepare debugging with

import collections
locals().update({str(k): ureg.Quantity(typhon.math.common.promote_maximally((v.sel(calibrated_channel=4) if "calibrated_channel" in v.dims else v)), v.units).ravel()[0] for (k, v) in self._make_adict_dims_consistent_if_needed(collections.ChainMap(self._quantities, self._other_quantities), me.symbols["R_e"]).items()})
R_e_full = me.expressions[me.symbols["R_e"]]
R_e_simp = me.expression_Re_simplified
from .measurement_equation import expressions, symbols
R_e_full2 = me.recursive_substitution(R_e_full, expressions=me.expressions, stop_at={symbols["T_IWCT"], symbols["h"], symbols["c"], symbols["k_b"]})

and compare:

• between complete radiances, -0.058% difference:

In : x = eval(str(R_e_full.subs({"R_selfIWCT": 0, "R_selfs": 0, "O_RIWCT": 0, "O_Re": 0, "R_refl": 0}))); y=R_e; print(x, y, x-y, "{:%}".format(((x-y)/x).m), sep="\n")
1.427965404115273e-12 watt / hertz / meter ** 2 / steradian
1.4287890253991072e-12 watt / hertz / meter ** 2 / steradian
-8.236212838342122e-16 watt / hertz / meter ** 2 / steradian
-0.057678%

• assuming linearity only, 6.40% difference:

In : x = eval(str(R_e_full.subs({"R_selfIWCT": 0, "R_selfs": 0, "O_RIWCT": 0, "O_Re": 0, "R_refl": 0, "a_2": 0}))); y=rad_wn_linear; print(x, y, x-y, "{:%}".format(((x-y)/x).m), sep="\n")
1.5131952299227614e-12 watt / hertz / meter ** 2 / steradian
1.4163937849327702e-12 watt / hertz / meter ** 2 / steradian
9.680144498999117e-14 watt / hertz / meter ** 2 / steradian
6.397155%

• assuming no emissivity correction only, -7.6% difference:

In : x = eval(str(R_e_full.subs({"R_selfIWCT": 0, "R_selfs": 0, "O_RIWCT": 0, "O_Re": 0, "R_refl": 0, "a_3": 0}))); y=rad_wn_noεcorr; print(x, y, x-y, "{:%}".format(((x-y)/x).m), sep="\n")
1.3041246534944941e-12 watt / hertz / meter ** 2 / steradian
1.4026476662769981e-12 watt / hertz / meter ** 2 / steradian
-9.852301278250398e-14 watt / hertz / meter ** 2 / steradian
-7.554724%

• neither emissivity correction nor nonlinearity, -0.064% difference:

In : x = eval(str(R_e_full.subs({"R_selfIWCT": 0, "R_selfs": 0, "O_RIWCT": 0, "O_Re": 0, "R_refl": 0, "a_2": 0, "a_3": 0}))); y=rad_wn_linearnoεcorr; print(x, y, x-y, "{:%}".format(((x-y)/x).m), sep="\n")
1.3893544793019825e-12 watt / hertz / meter ** 2 / steradian
1.390252425810661e-12 watt / hertz / meter ** 2 / steradian
-8.97946508678595e-16 watt / hertz / meter ** 2 / steradian
-0.064630%

Differences in the order of 0.06% may be compatible with the difference in band correction or direct calculations, but there is clearly a serious bug when I get a 6% difference. I trust the symbolic version better than the direct version, so I suspect there is a bug in the latter.

[gerritholl]

After commit f524535, two of the problems have swapped… now it goes wrong whenever the nonlinearity is taken out, but right for the case where only the emissivity correction is taken out…

In : x = eval(str(R_e_full.subs({"R_selfIWCT": 0, "R_selfs": 0, "O_RIWCT": 0, "O_Re": 0, "R_refl": 0}))); y=R_e; print(x, y, x-y, "{:%}".format(((x-y)/x).m), sep="\n")
1.427965404115273e-12 watt / hertz / meter ** 2 / steradian
1.4287890253991072e-12 watt / hertz / meter ** 2 / steradian
-8.236212838342122e-16 watt / hertz / meter ** 2 / steradian
-0.057678%
In : x = eval(str(R_e_full.subs({"R_selfIWCT": 0, "R_selfs": 0, "O_RIWCT": 0, "O_Re": 0, "R_refl": 0, "a_2": 0}))); y=rad_wn_linear; print(x, y, x-y, "{:%}".format(((x-y)/x).m), sep="\n")
1.5131952299227614e-12 watt / hertz / meter ** 2 / steradian
1.4287890253991072e-12 watt / hertz / meter ** 2 / steradian
8.440620452365414e-14 watt / hertz / meter ** 2 / steradian
5.578012%
In : x = eval(str(R_e_full.subs({"R_selfIWCT": 0, "R_selfs": 0, "O_RIWCT": 0, "O_Re": 0, "R_refl": 0, "a_3": 0}))); y=rad_wn_noεcorr; print(x, y, x-y, "{:%}".format(((x-y)/x).m), sep="\n")
1.3041246534944941e-12 watt / hertz / meter ** 2 / steradian
1.3049416770324787e-12 watt / hertz / meter ** 2 / steradian
-8.170235379846112e-16 watt / hertz / meter ** 2 / steradian
-0.062649%
In : x = eval(str(R_e_full.subs({"R_selfIWCT": 0, "R_selfs": 0, "O_RIWCT": 0, "O_Re": 0, "R_refl": 0, "a_2": 0, "a_3": 0}))); y=rad_wn_linearnoεcorr; print(x, y, x-y, "{:%}".format(((x-y)/x).m), sep="\n")
1.3893544793019825e-12 watt / hertz / meter ** 2 / steradian
1.3049416770324787e-12 watt / hertz / meter ** 2 / steradian
8.441280226950374e-14 watt / hertz / meter ** 2 / steradian
6.075685%