Olga Sergijenko / gammalib

Commit 1790bed6c915c47ff5e06ed47f03d7166b86bd51

Authored by Jürgen Knödlseder
1 parent 0b175024

Fix broken power law flux and eflux computation bug (#3279)

Inline Side-by-side
Showing 4 changed files with 19 additions and 13 deletions
ChangeLog
  View file @1790bed
1 -2020-07-20 1 +2020-07-22
2 2
3 * Bug fix version 1.7.1 released 3 * Bug fix version 1.7.1 released
4 ================================ 4 ================================
5 5
  6 + Fix broken power law flux and eflux computation bug (#3279)
  7 +
6 8
7 2020-07-02 9 2020-07-02
8 10
1 New Features and Important Changes in GammaLib 1.7.1 1 New Features and Important Changes in GammaLib 1.7.1
2 2
3 -20 July 2020 3 +22 July 2020
4 4
5 5
6 1. Introduction 6 1. Introduction
@@ -63,7 +63,10 @@ None @@ -63,7 +63,10 @@ None
63 63
64 8. Model module 64 8. Model module
65 --------------- 65 ---------------
66 -None 66 +The flux() and eflux() methods of the GModelSpectralBrokenPlaw class returned
  67 +a value of zero if one of the energy boundaries was identical to the break
  68 +energy. The reason was an incomplete coverage of the possible conditions that
  69 +left this special case out of the computation. This is fixed now (#3279).
67 70
68 71
69 9. Numerics module 72 9. Numerics module
doc/source/admin/release_history/1.7.1.rst
  View file @1790bed
1 .. _1.7.1: 1 .. _1.7.1:
2 2
3 -GammaLib 1.7.1 release (20 July 2020) 3 +GammaLib 1.7.1 release (22 July 2020)
4 ===================================== 4 =====================================
5 5
6 Overview 6 Overview
@@ -10,4 +10,5 @@ GammaLib 1.7.1 is a bug fix release for GammaLib release 1.7. @@ -10,4 +10,5 @@ GammaLib 1.7.1 is a bug fix release for GammaLib release 1.7.
10 10
11 The following bugs were fixed: 11 The following bugs were fixed:
12 12
13 -* ... 13 +* [`3279 <https://cta-redmine.irap.omp.eu/issues/3279>`_] -
  14 + Fix broken power law flux and eflux computation bug
src/model/GModelSpectralBrokenPlaw.cpp
  View file @1790bed
@@ -452,8 +452,8 @@ double GModelSpectralBrokenPlaw::flux(const GEnergy&amp; emin, @@ -452,8 +452,8 @@ double GModelSpectralBrokenPlaw::flux(const GEnergy&amp; emin,
452 m_index2.value())); 452 m_index2.value()));
453 } 453 }
454 454
455 - // Second case: breakenergy > emax:  
456 - else if (m_breakenergy.value() > emax.MeV()) { 455 + // Second case: breakenergy >= emax
  456 + else if (m_breakenergy.value() >= emax.MeV()) {
457 flux = m_norm.value() * 457 flux = m_norm.value() *
458 gammalib::plaw_photon_flux(emin.MeV(), 458 gammalib::plaw_photon_flux(emin.MeV(),
459 emax.MeV(), 459 emax.MeV(),
@@ -461,8 +461,8 @@ double GModelSpectralBrokenPlaw::flux(const GEnergy&amp; emin, @@ -461,8 +461,8 @@ double GModelSpectralBrokenPlaw::flux(const GEnergy&amp; emin,
461 m_index1.value()); 461 m_index1.value());
462 } 462 }
463 463
464 - // Third case breakenergy < emin:  
465 - else if (m_breakenergy.value() < emin.MeV()) { 464 + // Third case: breakenergy <= emin
  465 + else if (m_breakenergy.value() <= emin.MeV()) {
466 flux = m_norm.value() * 466 flux = m_norm.value() *
467 gammalib::plaw_photon_flux(emin.MeV(), 467 gammalib::plaw_photon_flux(emin.MeV(),
468 emax.MeV(), 468 emax.MeV(),
@@ -517,8 +517,8 @@ double GModelSpectralBrokenPlaw::eflux(const GEnergy&amp; emin, @@ -517,8 +517,8 @@ double GModelSpectralBrokenPlaw::eflux(const GEnergy&amp; emin,
517 m_index2.value())); 517 m_index2.value()));
518 } 518 }
519 519
520 - // Second case: breakenergy > emax:  
521 - else if (m_breakenergy.value() > emax.MeV()) { 520 + // Second case: breakenergy >= emax
  521 + else if (m_breakenergy.value() >= emax.MeV()) {
522 eflux = m_norm.value() * 522 eflux = m_norm.value() *
523 gammalib::plaw_energy_flux(emin.MeV(), 523 gammalib::plaw_energy_flux(emin.MeV(),
524 emax.MeV(), 524 emax.MeV(),
@@ -526,8 +526,8 @@ double GModelSpectralBrokenPlaw::eflux(const GEnergy&amp; emin, @@ -526,8 +526,8 @@ double GModelSpectralBrokenPlaw::eflux(const GEnergy&amp; emin,
526 m_index1.value()); 526 m_index1.value());
527 } 527 }
528 528
529 - // Third case breakenergy < emin:  
530 - else if (m_breakenergy.value() < emin.MeV()) { 529 + // Third case: breakenergy <= emin
  530 + else if (m_breakenergy.value() <= emin.MeV()) {
531 eflux = m_norm.value() * 531 eflux = m_norm.value() *
532 gammalib::plaw_energy_flux(emin.MeV(), 532 gammalib::plaw_energy_flux(emin.MeV(),
533 emax.MeV(), 533 emax.MeV(),