GammaLib

It turned out that for the model

<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<source_library title="source library">
  <source name="GammaCygni" type="DiffuseSource" tscalc="1">
    <spectrum type="BrokenPowerLaw">
      <parameter name="Prefactor" value="1.08949" error="0" scale="1e-15" min="1e-07" max="1000" free="1" />
      <parameter name="Index1" value="2.05" scale="-1" min="0" max="5" free="0" />
      <parameter name="Index2" value="2.79" error="0" scale="-1" min="0.01" max="1000" free="1" />
      <parameter name="BreakEnergy" value="1" scale="100000" min="0.01" max="1000" free="0" />
    </spectrum>
    <spatialModel type="DiffuseMap" file="gammacygni_map.fits">
      <parameter name="Normalization" value="1" scale="1" min="0.001" max="1000" free="0" />
    </spatialModel>
  </source>
</source_library>

the flux is zero when the break energy and the energy boundary are identical:

>>> print(m[0].spectral().flux(gammalib.GEnergy(0.1,'TeV'),gammalib.GEnergy(100.0,'TeV')))
0.0
>>> print(m[0].spectral().flux(gammalib.GEnergy(0.09,'TeV'),gammalib.GEnergy(100.0,'TeV')))
7.30030511938e-11
>>> print(m[0].spectral().flux(gammalib.GEnergy(0.0999,'TeV'),gammalib.GEnergy(100.0,'TeV')))
6.09741642757e-11
>>> print(m[0].spectral().flux(gammalib.GEnergy(0.11,'TeV'),gammalib.GEnergy(100.0,'TeV')))
5.13186370633e-11
>>> print(m[0].spectral().flux(gammalib.GEnergy(0.10,'TeV'),gammalib.GEnergy(100.0,'TeV')))
0.0