The ParaMonte Documentation Website
Current view: top level - kernel - BandSpectrum_mod.f90 (source / functions) Hit Total Coverage
Test: ParaMonte 1.5.1 :: MPI Parallel Kernel - Code Coverage Report Lines: 41 41 100.0 %
Date: 2021-01-08 13:07:16 Functions: 5 5 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       2             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       3             : !!!!
       4             : !!!!   MIT License
       5             : !!!!
       6             : !!!!   ParaMonte: plain powerful parallel Monte Carlo library.
       7             : !!!!
       8             : !!!!   Copyright (C) 2012-present, The Computational Data Science Lab
       9             : !!!!
      10             : !!!!   This file is part of the ParaMonte library.
      11             : !!!!
      12             : !!!!   Permission is hereby granted, free of charge, to any person obtaining a
      13             : !!!!   copy of this software and associated documentation files (the "Software"),
      14             : !!!!   to deal in the Software without restriction, including without limitation
      15             : !!!!   the rights to use, copy, modify, merge, publish, distribute, sublicense,
      16             : !!!!   and/or sell copies of the Software, and to permit persons to whom the
      17             : !!!!   Software is furnished to do so, subject to the following conditions:
      18             : !!!!
      19             : !!!!   The above copyright notice and this permission notice shall be
      20             : !!!!   included in all copies or substantial portions of the Software.
      21             : !!!!
      22             : !!!!   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
      23             : !!!!   EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
      24             : !!!!   MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
      25             : !!!!   IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
      26             : !!!!   DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
      27             : !!!!   OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE
      28             : !!!!   OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
      29             : !!!!
      30             : !!!!   ACKNOWLEDGMENT
      31             : !!!!
      32             : !!!!   ParaMonte is an honor-ware and its currency is acknowledgment and citations.
      33             : !!!!   As per the ParaMonte library license agreement terms, if you use any parts of
      34             : !!!!   this library for any purposes, kindly acknowledge the use of ParaMonte in your
      35             : !!!!   work (education/research/industry/development/...) by citing the ParaMonte
      36             : !!!!   library as described on this page:
      37             : !!!!
      38             : !!!!       https://github.com/cdslaborg/paramonte/blob/main/ACKNOWLEDGMENT.md
      39             : !!!!
      40             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      41             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      42             : 
      43             : !> \brief
      44             : !> This module contains the procedures for computing the various Band Spectrum properties of GRBs.
      45             : !> See Eqn. A6 of [Shahmoradi and Nemiroff, 2015, MNRAS, 451:4645-4662](https://www.cdslab.org/pubs/2015_Shahmoradi_I.pdf).
      46             : !>
      47             : !> \author
      48             : !> Amir Shahmoradi, Tuesday April 30, 2019, 12:58 PM, SEIR, UTA
      49             : 
      50             : module BandSpectrum_mod
      51             : 
      52             :     use Constants_mod, only: IK, RK
      53             : 
      54             :     implicit none
      55             : 
      56             :     character(*), parameter :: MODULE_NAME = "@BandSpectrum_mod"
      57             : 
      58             :     real(RK), parameter     :: AMPLITUDE_DEFAULT = 1._RK    ! default amplitude of the Band function
      59             :     real(RK), parameter     :: ALPHA_DEFAULT = -1.1_RK      ! default low-energy index of the Band function
      60             :     real(RK), parameter     :: BETA_DEFAULT = -2.3_RK       ! default high-energy index of the Band function
      61             : 
      62             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      63             : 
      64             :     interface
      65             :     module subroutine getEnergyFluence(lowerLim,upperLim,epk,alpha,beta,tolerance,energyFluence,Err)
      66             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
      67             :         !DEC$ ATTRIBUTES DLLEXPORT :: getEnergyFluence
      68             : #endif
      69             :         use Constants_mod, only: IK, RK, HUGE_RK
      70             :         use QuadPackSPR_mod, only: qag
      71             :         use Err_mod, only: Err_type
      72             :         implicit none
      73             :         real(RK), intent(in)            :: lowerLim, upperLim, epk, alpha, beta, tolerance
      74             :         real(RK), intent(out)           :: energyFluence
      75             :         type(Err_type), intent(out)     :: Err
      76             :     end subroutine getEnergyFluence
      77             :     end interface
      78             : 
      79             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      80             : 
      81             :     interface
      82             :     module subroutine getPhotonFluence(lowerLim,upperLim,epk,alpha,beta,tolerance,photonFluence,Err)
      83             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
      84             :         !DEC$ ATTRIBUTES DLLEXPORT :: getPhotonFluence
      85             : #endif
      86             : 
      87             :        !use Integration_mod, only: doQuadRombClosed
      88             :         use Constants_mod, only: IK, RK, HUGE_RK
      89             :         use QuadPackSPR_mod, only: qag
      90             :         use Err_mod, only: Err_type
      91             :         implicit none
      92             :         real(RK), intent(in)            :: lowerLim, upperLim, epk, alpha, beta, tolerance
      93             :         real(RK), intent(out)           :: photonFluence
      94             :         type(Err_type), intent(out)     :: Err
      95             :     end subroutine getPhotonFluence
      96             :     end interface
      97             : 
      98             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      99             : 
     100             : contains
     101             : 
     102             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     103             : 
     104             :     !> \brief
     105             :     !> Return `epk*(alpha-beta)/(alpha+2)`, which is the normalization
     106             :     !> factor in the exponent of the first component of the Band model.
     107             :     !> \param[in]   epk     :   The spectral peak energy in units of [keV].
     108             :     !> \param[in]   alpha   :   The lower spectral exponent of the Band model.
     109             :     !> \param[in]   beta    :   The upper spectral exponent of the Band model.
     110             :     !>
     111             :     !> \return
     112             :     !> `ebrk` : The spectral break energy in units of [keV].
     113           3 :     pure function getEbreak(epk,alpha,beta) result(ebrk)
     114             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     115             :         !DEC$ ATTRIBUTES DLLEXPORT :: getEbreak
     116             : #endif
     117             :         use Constants_mod, only: RK
     118             :         implicit none
     119             :         real(RK), intent(in)            :: epk
     120             :         real(RK), intent(in), optional  :: alpha, beta
     121             :         real(RK)                        :: ebrk
     122           3 :         ebrk = epk * (alpha-beta) / (alpha+2._RK)
     123           6 :     end function getEbreak
     124             : 
     125             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     126             : 
     127             :     !> \brief
     128             :     !> Compute the break energy `ebrk` of the Band model, as well as the coefficient by which the high energy component of the
     129             :     !> Band model must be multiplied in order to get a smooth function (assuming the normalization coefficient of the lower
     130             :     !> energy component is unity).
     131             :     !>
     132             :     !> \param[in]   epk             :   The spectral peak energy in units of [keV].
     133             :     !> \param[in]   alpha           :   The lower spectral exponent of the Band model.
     134             :     !> \param[in]   beta            :   The upper spectral exponent of the Band model.
     135             :     !> \param[out]  ebrk            :   The spectral break energy in units of [keV].
     136             :     !> \param[out]  coef            :   The spectral continuity coefficient.
     137             :     !> \param[out]  alphaPlusTwo    :   The lower spectral exponent of the Band model plus two.
     138           3 :     pure subroutine getBandParam(epk,alpha,beta,ebrk,coef,alphaPlusTwo)
     139             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     140             :         !DEC$ ATTRIBUTES DLLEXPORT :: getBandParam
     141             : #endif
     142           3 :         use Constants_mod, only: IK, RK
     143             :         implicit none
     144             :         real(RK), intent(in)    :: epk,alpha,beta
     145             :         real(RK), intent(out)   :: ebrk, coef, alphaPlusTwo
     146           3 :         real(RK)                :: alphaMinusBeta
     147           3 :         alphaPlusTwo = alpha + 2._RK
     148           3 :         alphaMinusBeta = alpha - beta
     149           3 :         ebrk = epk * alphaMinusBeta / alphaPlusTwo
     150           3 :         coef = ebrk**alphaMinusBeta * exp(-alphaMinusBeta)
     151           3 :     end subroutine getBandParam
     152             : 
     153             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     154             : 
     155             :     !> \brief
     156             :     !> Compute the differential photon flux according to the Band differential spectrum at the given input `energy` value.
     157             :     !>
     158             :     !> \param[in]   energy          :   The energy (in units of [keV]) at which the differential photon flux must be computed.
     159             :     !> \param[in]   epk             :   The spectral peak energy in units of [keV].
     160             :     !> \param[in]   alpha           :   The lower spectral exponent of the Band model.
     161             :     !> \param[in]   beta            :   The upper spectral exponent of the Band model.
     162             :     !> \param[in]   ebrk            :   The spectral break energy in units of [keV].
     163             :     !> \param[in]   coef            :   The spectral continuity coefficient.
     164             :     !> \param[in]   alphaPlusTwo    :   The lower spectral exponent of the Band model plus two.
     165             :     !>
     166             :     !> \return
     167             :     !> `photonFlux` : The energy flux in units of photon counts.
     168             :     !>
     169             :     !> \warning
     170             :     !> A negative huge output value is used to signal error has occurred. Under normal conditions, the output is always positive.
     171             :     !>
     172             :     !> \warning
     173             :     !> The input energy values `energy`, `epk`, `ebrk`, must all have the same units.
     174             :     !> It is expected that the input energy is in units of keV, although it does not affect the accuracy of the results.
     175          12 :     pure function getPhotonFlux(energy,epk,alpha,beta,ebrk,coef,alphaPlusTwo) result(photonFlux)
     176             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     177             :         !DEC$ ATTRIBUTES DLLEXPORT :: getPhotonFlux
     178             : #endif
     179           3 :         use Constants_mod, only: IK, RK, HUGE_RK
     180             :         implicit none
     181             :         real(RK), intent(in)            :: energy, epk, ebrk, alpha, beta,coef,alphaPlusTwo
     182             :         real(RK)                        :: photonFlux
     183             :         ! check if the photon indices are consistent with the mathematical rules
     184          12 :         if (alpha<beta .or. alpha<-2._RK) then
     185           6 :             photonFlux = -HUGE_RK
     186           6 :             return
     187             :         end if
     188             :         ! compute the spectrum
     189           6 :         if (energy<=ebrk) then
     190           3 :             photonFlux = energy**alpha * exp(-energy*alphaPlusTwo/epk)
     191             :         else
     192           3 :             photonFlux = coef * energy**beta
     193             :         end if
     194          12 :     end function getPhotonFlux
     195             : 
     196             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     197             : 
     198             :     !> \brief
     199             :     !> Compute the differential photon flux according to the lower component of the Band spectrum at the given input `energy` value.
     200             :     !>
     201             :     !> \param[in]   energy              :   The energy (in units of [keV]) at which the differential photon flux must be computed.
     202             :     !> \param[in]   alpha               :   The lower spectral exponent of the Band model.
     203             :     !> \param[in]   alphaPlusTwoOverEpk :   The lower spectral exponent of the Band model plus two.
     204             :     !>
     205             :     !> \return
     206             :     !> `photonFluxLower` : The energy flux in units of photon counts.
     207             :     !>
     208             :     !> \warning
     209             :     !> The input `energy` values must be less than `ebrk`.
     210           3 :     pure function getPhotonFluxLower(energy,alpha,alphaPlusTwoOverEpk) result(photonFluxLower)
     211             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     212             :         !DEC$ ATTRIBUTES DLLEXPORT :: getPhotonFluxLower
     213             : #endif
     214          12 :         use Constants_mod, only: IK, RK, HUGE_RK
     215             :         implicit none
     216             :         real(RK), intent(in)            :: energy, alpha, alphaPlusTwoOverEpk
     217             :         real(RK)                        :: photonFluxLower
     218           3 :         photonFluxLower = energy**alpha * exp(-energy*alphaPlusTwoOverEpk)  ! the lower-energy spectrum
     219           6 :     end function getPhotonFluxLower
     220             : 
     221             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     222             : 
     223             : ! Bizarrely and frustratingly, Microsoft Windows Subsystem for Linux Ubuntu with GFortran yields Segmentation faults with internal procedure calls.
     224             : ! This is so unfortunate. To bypass this issue for now, the following subroutine is implemented as separate submodule 
     225             : ! so that the internal shared parameters can be safely passed as submodule parameters.
     226             : 
     227             : !WSL_GFORTRAN_BUG     !> \brief
     228             : !WSL_GFORTRAN_BUG     !> Integrate the Band differential spectrum over the input energy range.
     229             : !WSL_GFORTRAN_BUG     !>
     230             : !WSL_GFORTRAN_BUG     !> \param[in]   lowerLim        :   The lower limit energy (in units of [keV]) of the integration.
     231             : !WSL_GFORTRAN_BUG     !> \param[in]   upperLim        :   The upper limit energy (in units of [keV]) of the integration.
     232             : !WSL_GFORTRAN_BUG     !> \param[in]   epk             :   The spectral peak energy in units of [keV].
     233             : !WSL_GFORTRAN_BUG     !> \param[in]   alpha           :   The lower spectral exponent of the Band model.
     234             : !WSL_GFORTRAN_BUG     !> \param[in]   beta            :   The upper spectral exponent of the Band model.
     235             : !WSL_GFORTRAN_BUG     !> \param[in]   tolerance       :   The relative accuracy tolerance of the integration.
     236             : !WSL_GFORTRAN_BUG     !> \param[out]  photonFluence   :   The fluence in units of photon counts within the input energy range.
     237             : !WSL_GFORTRAN_BUG     !> \param[out]  Err             :   An object of class [Err_type](@ref err_mod::err_type) containing error-handling information.
     238             : !WSL_GFORTRAN_BUG     subroutine getPhotonFluence(lowerLim,upperLim,epk,alpha,beta,tolerance,photonFluence,Err)
     239             : !WSL_GFORTRAN_BUG #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     240             : !WSL_GFORTRAN_BUG         !DEC$ ATTRIBUTES DLLEXPORT :: getPhotonFluence
     241             : !WSL_GFORTRAN_BUG #endif
     242             : !WSL_GFORTRAN_BUG 
     243             : !WSL_GFORTRAN_BUG        !use Integration_mod, only: doQuadRombClosed
     244             : !WSL_GFORTRAN_BUG         use Constants_mod, only: IK, RK, HUGE_RK
     245             : !WSL_GFORTRAN_BUG         use QuadPackSPR_mod, only: qag
     246             : !WSL_GFORTRAN_BUG         use Err_mod, only: Err_type
     247             : !WSL_GFORTRAN_BUG         implicit none
     248             : !WSL_GFORTRAN_BUG         real(RK), intent(in)            :: lowerLim, upperLim, epk, alpha, beta, tolerance
     249             : !WSL_GFORTRAN_BUG         real(RK), intent(out)           :: photonFluence
     250             : !WSL_GFORTRAN_BUG         type(Err_type), intent(out)     :: Err
     251             : !WSL_GFORTRAN_BUG         character(*), parameter         :: PROCEDURE_NAME = "@getPhotonFluence()"
     252             : !WSL_GFORTRAN_BUG         real(RK)                        :: ebrk, alphaPlusTwo
     253             : !WSL_GFORTRAN_BUG         real(RK)                        :: thisUpperLim, alphaPlusTwoOverEpk, betaPlusOne
     254             : !WSL_GFORTRAN_BUG         real(RK)                        :: alphaMinusBeta, coef
     255             : !WSL_GFORTRAN_BUG         real(RK)                        :: abserr
     256             : !WSL_GFORTRAN_BUG         integer(IK)                     :: neval
     257             : !WSL_GFORTRAN_BUG         integer(IK)                     :: ierr
     258             : !WSL_GFORTRAN_BUG        !real(RK)                        :: alphaPlusOne, logGammaAlphaPlusOne
     259             : !WSL_GFORTRAN_BUG 
     260             : !WSL_GFORTRAN_BUG         Err%occurred = .false.
     261             : !WSL_GFORTRAN_BUG 
     262             : !WSL_GFORTRAN_BUG         if (lowerLim>=upperLim) then
     263             : !WSL_GFORTRAN_BUG             photonFluence = 0._RK
     264             : !WSL_GFORTRAN_BUG             return
     265             : !WSL_GFORTRAN_BUG         end if
     266             : !WSL_GFORTRAN_BUG 
     267             : !WSL_GFORTRAN_BUG         ! check if the photon indices are consistent with the mathematical rules
     268             : !WSL_GFORTRAN_BUG         if (alpha<beta .or. alpha<-2._RK) then
     269             : !WSL_GFORTRAN_BUG             photonFluence = -HUGE_RK
     270             : !WSL_GFORTRAN_BUG             Err%occurred = .true.
     271             : !WSL_GFORTRAN_BUG             Err%msg = MODULE_NAME // PROCEDURE_NAME // ": Error occurred: alpha<beta .or. alpha<-2._RK"
     272             : !WSL_GFORTRAN_BUG             return
     273             : !WSL_GFORTRAN_BUG         end if
     274             : !WSL_GFORTRAN_BUG 
     275             : !WSL_GFORTRAN_BUG         ! integrate the spectrum
     276             : !WSL_GFORTRAN_BUG         alphaPlusTwo = alpha + 2._RK
     277             : !WSL_GFORTRAN_BUG         alphaMinusBeta = alpha - beta
     278             : !WSL_GFORTRAN_BUG         ebrk = epk*alphaMinusBeta/alphaPlusTwo
     279             : !WSL_GFORTRAN_BUG 
     280             : !WSL_GFORTRAN_BUG         if (lowerLim>ebrk) then
     281             : !WSL_GFORTRAN_BUG 
     282             : !WSL_GFORTRAN_BUG             ! there is only the high energy component in the photonFluence
     283             : !WSL_GFORTRAN_BUG             betaPlusOne = beta + 1._RK
     284             : !WSL_GFORTRAN_BUG             coef = ebrk**alphaMinusBeta * exp(-alphaMinusBeta);
     285             : !WSL_GFORTRAN_BUG             photonFluence = coef * ( upperLim**betaPlusOne - lowerLim**betaPlusOne ) / betaPlusOne
     286             : !WSL_GFORTRAN_BUG             return
     287             : !WSL_GFORTRAN_BUG 
     288             : !WSL_GFORTRAN_BUG !#if defined OS_IS_WSL && defined CODECOV_ENABLED
     289             : !WSL_GFORTRAN_BUG !! LCOV_EXCL_START
     290             : !WSL_GFORTRAN_BUG !#endif
     291             : !WSL_GFORTRAN_BUG 
     292             : !WSL_GFORTRAN_BUG         elseif (lowerLim<ebrk) then
     293             : !WSL_GFORTRAN_BUG 
     294             : !WSL_GFORTRAN_BUG             alphaPlusTwoOverEpk = alphaPlusTwo / epk
     295             : !WSL_GFORTRAN_BUG             thisUpperLim = min(upperLim,ebrk)
     296             : !WSL_GFORTRAN_BUG             !alphaPlusOne = alpha + 1._RK
     297             : !WSL_GFORTRAN_BUG             !if (alpha>-1._RK) then
     298             : !WSL_GFORTRAN_BUG             !    logGammaAlphaPlusOne = log_gamma( alphaPlusOne )
     299             : !WSL_GFORTRAN_BUG             !    ! use the analytical approach to compute the photonFluence:
     300             : !WSL_GFORTRAN_BUG             !    ! https://www.wolframalpha.com/input/?i=integrate+x%5Ea+*+exp(-b*x)
     301             : !WSL_GFORTRAN_BUG             !    photonFluence = getUpperGamma( exponent = alphaPlusOne &
     302             : !WSL_GFORTRAN_BUG             !                            , logGammaExponent = logGammaAlphaPlusOne &
     303             : !WSL_GFORTRAN_BUG             !                            , lowerLim = alphaPlusTwoOverEpk * lowerLim &
     304             : !WSL_GFORTRAN_BUG             !                            , tolerance = tolerance &
     305             : !WSL_GFORTRAN_BUG             !                            ) &
     306             : !WSL_GFORTRAN_BUG             !             - getUpperGamma( exponent = alphaPlusOne &
     307             : !WSL_GFORTRAN_BUG             !                            , logGammaExponent = logGammaAlphaPlusOne &
     308             : !WSL_GFORTRAN_BUG             !                            , lowerLim = alphaPlusTwoOverEpk * thisUpperLim &
     309             : !WSL_GFORTRAN_BUG             !                            , tolerance = tolerance &
     310             : !WSL_GFORTRAN_BUG             !                            )
     311             : !WSL_GFORTRAN_BUG             !    photonFluence = photonFluence / alphaPlusTwoOverEpk**alphaPlusOne
     312             : !WSL_GFORTRAN_BUG             !else
     313             : !WSL_GFORTRAN_BUG                 ! use brute-force integration
     314             : !WSL_GFORTRAN_BUG                 call qag( f             = getBandCompLowPhoton  &
     315             : !WSL_GFORTRAN_BUG                         , a             = lowerLim              &
     316             : !WSL_GFORTRAN_BUG                         , b             = thisUpperLim          &
     317             : !WSL_GFORTRAN_BUG                         , epsabs        = 0._RK                 &
     318             : !WSL_GFORTRAN_BUG                         , epsrel        = tolerance             &
     319             : !WSL_GFORTRAN_BUG                         , key           = 1_IK                  &
     320             : !WSL_GFORTRAN_BUG                         , result        = photonFluence         &
     321             : !WSL_GFORTRAN_BUG                         , abserr        = abserr                &
     322             : !WSL_GFORTRAN_BUG                         , neval         = neval                 &
     323             : !WSL_GFORTRAN_BUG                         , ier           = ierr                  &
     324             : !WSL_GFORTRAN_BUG                         )
     325             : !WSL_GFORTRAN_BUG                 !write(*,*) neval
     326             : !WSL_GFORTRAN_BUG                 !call doQuadRombClosed   ( getFunc       = getBandCompLowPhoton &
     327             : !WSL_GFORTRAN_BUG                 !                        , xmin          = lowerLim             &
     328             : !WSL_GFORTRAN_BUG                 !                        , xmax          = thisUpperLim         &
     329             : !WSL_GFORTRAN_BUG                 !                        , tolerance     = 1.e-7_RK             &
     330             : !WSL_GFORTRAN_BUG                 !                        , nRefinement   = 10_IK                &
     331             : !WSL_GFORTRAN_BUG                 !                        , photonFluence = photonFluence        &
     332             : !WSL_GFORTRAN_BUG                 !                        , ierr          = ierr                 &
     333             : !WSL_GFORTRAN_BUG                 !                        )
     334             : !WSL_GFORTRAN_BUG                 if (ierr/=0_IK) then
     335             : !WSL_GFORTRAN_BUG                     photonFluence = -HUGE_RK
     336             : !WSL_GFORTRAN_BUG                     Err%occurred = .true.
     337             : !WSL_GFORTRAN_BUG                     Err%stat = ierr
     338             : !WSL_GFORTRAN_BUG                     Err%msg = MODULE_NAME // PROCEDURE_NAME // ": Error occurred at QuadPack routine. Check the error code to identify the root cause."
     339             : !WSL_GFORTRAN_BUG                     return
     340             : !WSL_GFORTRAN_BUG                 end if
     341             : !WSL_GFORTRAN_BUG             !end if
     342             : !WSL_GFORTRAN_BUG 
     343             : !WSL_GFORTRAN_BUG             if (upperLim>ebrk) then
     344             : !WSL_GFORTRAN_BUG                 ! add the remaining part of the photonFluence from the high-energy component
     345             : !WSL_GFORTRAN_BUG                 betaPlusOne = beta + 1._RK
     346             : !WSL_GFORTRAN_BUG                 alphaMinusBeta = alpha - beta
     347             : !WSL_GFORTRAN_BUG                 coef = ebrk**alphaMinusBeta * exp(-alphaMinusBeta);
     348             : !WSL_GFORTRAN_BUG                 photonFluence = photonFluence + coef * ( upperLim**betaPlusOne - ebrk**betaPlusOne ) / betaPlusOne
     349             : !WSL_GFORTRAN_BUG                 return
     350             : !WSL_GFORTRAN_BUG             end if
     351             : !WSL_GFORTRAN_BUG 
     352             : !WSL_GFORTRAN_BUG         end if
     353             : !WSL_GFORTRAN_BUG 
     354             : !WSL_GFORTRAN_BUG     contains
     355             : !WSL_GFORTRAN_BUG 
     356             : !WSL_GFORTRAN_BUG         pure function getBandCompLowPhoton(energy) result(bandCompLow)
     357             : !WSL_GFORTRAN_BUG             implicit none
     358             : !WSL_GFORTRAN_BUG             real(RK), intent(in)    :: energy
     359             : !WSL_GFORTRAN_BUG             real(RK)                :: bandCompLow
     360             : !WSL_GFORTRAN_BUG             bandCompLow = energy**alpha * exp(-alphaPlusTwoOverEpk*energy)
     361             : !WSL_GFORTRAN_BUG         end function getBandCompLowPhoton
     362             : !WSL_GFORTRAN_BUG 
     363             : !WSL_GFORTRAN_BUG !#if defined OS_IS_WSL && defined CODECOV_ENABLED
     364             : !WSL_GFORTRAN_BUG !! LCOV_EXCL_STOP
     365             : !WSL_GFORTRAN_BUG !#endif
     366             : !WSL_GFORTRAN_BUG 
     367             : !WSL_GFORTRAN_BUG     end subroutine getPhotonFluence
     368             : 
     369             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     370             : 
     371             : ! Bizarrely and frustratingly, Microsoft Windows Subsystem for Linux Ubuntu with GFortran yields Segmentation faults with internal procedure calls.
     372             : ! This is so unfortunate. To bypass this issue for now, the following subroutine is implemented as separate submodule 
     373             : ! so that the internal shared parameters can be safely passed as submodule parameters.
     374             : 
     375             : !WSL_GFORTRAN_BUG     !> \brief
     376             : !WSL_GFORTRAN_BUG     !> Integrate the Band differential spectrum over the input energy range in units of the input energy.
     377             : !WSL_GFORTRAN_BUG     !>
     378             : !WSL_GFORTRAN_BUG     !> \param[in]   lowerLim        :   The lower limit energy (in units of [keV]) of the integration.
     379             : !WSL_GFORTRAN_BUG     !> \param[in]   upperLim        :   The upper limit energy (in units of [keV]) of the integration.
     380             : !WSL_GFORTRAN_BUG     !> \param[in]   epk             :   The spectral peak energy in units of [keV].
     381             : !WSL_GFORTRAN_BUG     !> \param[in]   alpha           :   The lower spectral exponent of the Band model.
     382             : !WSL_GFORTRAN_BUG     !> \param[in]   beta            :   The upper spectral exponent of the Band model.
     383             : !WSL_GFORTRAN_BUG     !> \param[in]   tolerance       :   The relative accuracy tolerance of the integration.
     384             : !WSL_GFORTRAN_BUG     !> \param[out]  energyFluence   :   The fluence in units of the input energy (keV) within the input energy range.
     385             : !WSL_GFORTRAN_BUG     !> \param[out]  Err             :   An object of class [Err_type](@ref err_mod::err_type) containing error-handling information.
     386             : !WSL_GFORTRAN_BUG     subroutine getEnergyFluence(lowerLim,upperLim,epk,alpha,beta,tolerance,energyFluence,Err)
     387             : !WSL_GFORTRAN_BUG #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     388             : !WSL_GFORTRAN_BUG         !DEC$ ATTRIBUTES DLLEXPORT :: getEnergyFluence
     389             : !WSL_GFORTRAN_BUG #endif
     390             : !WSL_GFORTRAN_BUG 
     391             : !WSL_GFORTRAN_BUG         use Constants_mod, only: IK, RK, HUGE_RK
     392             : !WSL_GFORTRAN_BUG         use QuadPackSPR_mod, only: qag
     393             : !WSL_GFORTRAN_BUG         use Err_mod, only: Err_type
     394             : !WSL_GFORTRAN_BUG         implicit none
     395             : !WSL_GFORTRAN_BUG         real(RK), intent(in)            :: lowerLim, upperLim, epk, alpha, beta, tolerance
     396             : !WSL_GFORTRAN_BUG         type(Err_type), intent(out)     :: Err
     397             : !WSL_GFORTRAN_BUG         real(RK), intent(out)           :: energyFluence
     398             : !WSL_GFORTRAN_BUG         character(*), parameter         :: PROCEDURE_NAME = "@getEnergyFluence()"
     399             : !WSL_GFORTRAN_BUG         real(RK)                        :: ebrk, alphaPlusTwo
     400             : !WSL_GFORTRAN_BUG         real(RK)                        :: thisUpperLim, alphaPlusTwoOverEpk, betaPlusTwo
     401             : !WSL_GFORTRAN_BUG         real(RK)                        :: alphaMinusBeta, coef, alphaPlusOne
     402             : !WSL_GFORTRAN_BUG         real(RK)                        :: abserr
     403             : !WSL_GFORTRAN_BUG         integer(IK)                     :: neval
     404             : !WSL_GFORTRAN_BUG         integer(IK)                     :: ierr
     405             : !WSL_GFORTRAN_BUG 
     406             : !WSL_GFORTRAN_BUG         Err%occurred = .false.
     407             : !WSL_GFORTRAN_BUG 
     408             : !WSL_GFORTRAN_BUG         if (lowerLim>=upperLim) then
     409             : !WSL_GFORTRAN_BUG             energyFluence = 0._RK
     410             : !WSL_GFORTRAN_BUG             return
     411             : !WSL_GFORTRAN_BUG         end if
     412             : !WSL_GFORTRAN_BUG 
     413             : !WSL_GFORTRAN_BUG         ! check if the photon indices are consistent with the mathematical rules
     414             : !WSL_GFORTRAN_BUG         if (alpha<beta .or. alpha<-2._RK) then
     415             : !WSL_GFORTRAN_BUG             energyFluence = -HUGE_RK
     416             : !WSL_GFORTRAN_BUG             Err%occurred = .true.
     417             : !WSL_GFORTRAN_BUG             Err%msg = MODULE_NAME // PROCEDURE_NAME // ": Error occurred: alpha<beta .or. alpha<-2._RK"
     418             : !WSL_GFORTRAN_BUG             return
     419             : !WSL_GFORTRAN_BUG         end if
     420             : !WSL_GFORTRAN_BUG 
     421             : !WSL_GFORTRAN_BUG         ! integrate the spectrum
     422             : !WSL_GFORTRAN_BUG         alphaPlusTwo = alpha + 2._RK
     423             : !WSL_GFORTRAN_BUG         alphaMinusBeta = alpha - beta
     424             : !WSL_GFORTRAN_BUG         ebrk = epk*alphaMinusBeta/alphaPlusTwo
     425             : !WSL_GFORTRAN_BUG 
     426             : !WSL_GFORTRAN_BUG         if (lowerLim>ebrk) then
     427             : !WSL_GFORTRAN_BUG 
     428             : !WSL_GFORTRAN_BUG             ! there is only the high energy component in the energyFluence
     429             : !WSL_GFORTRAN_BUG             betaPlusTwo = beta + 2._RK
     430             : !WSL_GFORTRAN_BUG             coef = ebrk**alphaMinusBeta * exp(-alphaMinusBeta);
     431             : !WSL_GFORTRAN_BUG             energyFluence = coef * ( upperLim**betaPlusTwo - lowerLim**betaPlusTwo ) / betaPlusTwo
     432             : !WSL_GFORTRAN_BUG             return
     433             : !WSL_GFORTRAN_BUG 
     434             : !WSL_GFORTRAN_BUG !#if defined OS_IS_WSL && defined CODECOV_ENABLED
     435             : !WSL_GFORTRAN_BUG !! LCOV_EXCL_START
     436             : !WSL_GFORTRAN_BUG !#endif
     437             : !WSL_GFORTRAN_BUG 
     438             : !WSL_GFORTRAN_BUG         elseif (lowerLim<ebrk) then
     439             : !WSL_GFORTRAN_BUG 
     440             : !WSL_GFORTRAN_BUG             alphaPlusTwoOverEpk = alphaPlusTwo / epk
     441             : !WSL_GFORTRAN_BUG             thisUpperLim = min(upperLim,ebrk)
     442             : !WSL_GFORTRAN_BUG             alphaPlusOne = alpha + 1._RK
     443             : !WSL_GFORTRAN_BUG             ! use brute-force integration
     444             : !WSL_GFORTRAN_BUG             call qag( f             = getBandCompLowEnergy  &
     445             : !WSL_GFORTRAN_BUG                     , a             = lowerLim              &
     446             : !WSL_GFORTRAN_BUG                     , b             = thisUpperLim          &
     447             : !WSL_GFORTRAN_BUG                     , epsabs        = 0._RK                 &
     448             : !WSL_GFORTRAN_BUG                     , epsrel        = tolerance             &
     449             : !WSL_GFORTRAN_BUG                     , key           = 1_IK                  &
     450             : !WSL_GFORTRAN_BUG                     , result        = energyFluence         &
     451             : !WSL_GFORTRAN_BUG                     , abserr        = abserr                &
     452             : !WSL_GFORTRAN_BUG                     , neval         = neval                 &
     453             : !WSL_GFORTRAN_BUG                     , ier           = ierr                  &
     454             : !WSL_GFORTRAN_BUG                     )
     455             : !WSL_GFORTRAN_BUG 
     456             : !WSL_GFORTRAN_BUG             if (ierr/=0_IK) then
     457             : !WSL_GFORTRAN_BUG                 energyFluence = -HUGE_RK
     458             : !WSL_GFORTRAN_BUG                 Err%occurred = .true.
     459             : !WSL_GFORTRAN_BUG                 Err%stat = ierr
     460             : !WSL_GFORTRAN_BUG                 Err%msg = MODULE_NAME // PROCEDURE_NAME // ": Error occurred at QuadPack routine. Check the error code to identify the root cause."
     461             : !WSL_GFORTRAN_BUG                 return
     462             : !WSL_GFORTRAN_BUG             end if
     463             : !WSL_GFORTRAN_BUG 
     464             : !WSL_GFORTRAN_BUG             if (upperLim>ebrk) then ! add the remaining part of the energyFluence from the high-energy component
     465             : !WSL_GFORTRAN_BUG                 betaPlusTwo = beta + 2._RK
     466             : !WSL_GFORTRAN_BUG                 alphaMinusBeta = alpha - beta
     467             : !WSL_GFORTRAN_BUG                 coef = ebrk**alphaMinusBeta * exp(-alphaMinusBeta);
     468             : !WSL_GFORTRAN_BUG                 energyFluence = energyFluence + coef * ( upperLim**betaPlusTwo - ebrk**betaPlusTwo ) / betaPlusTwo
     469             : !WSL_GFORTRAN_BUG                 return
     470             : !WSL_GFORTRAN_BUG             end if
     471             : !WSL_GFORTRAN_BUG 
     472             : !WSL_GFORTRAN_BUG         end if
     473             : !WSL_GFORTRAN_BUG 
     474             : !WSL_GFORTRAN_BUG     contains
     475             : !WSL_GFORTRAN_BUG 
     476             : !WSL_GFORTRAN_BUG         pure function getBandCompLowEnergy(energy) result(bandCompLow)
     477             : !WSL_GFORTRAN_BUG             implicit none
     478             : !WSL_GFORTRAN_BUG             real(RK), intent(in)    :: energy
     479             : !WSL_GFORTRAN_BUG             real(RK)                :: bandCompLow
     480             : !WSL_GFORTRAN_BUG             bandCompLow = energy**alphaPlusOne * exp(-alphaPlusTwoOverEpk*energy)
     481             : !WSL_GFORTRAN_BUG         end function getBandCompLowEnergy
     482             : !WSL_GFORTRAN_BUG 
     483             : !WSL_GFORTRAN_BUG !#if defined OS_IS_WSL && defined CODECOV_ENABLED
     484             : !WSL_GFORTRAN_BUG !! LCOV_EXCL_STOP
     485             : !WSL_GFORTRAN_BUG !#endif
     486             : !WSL_GFORTRAN_BUG 
     487             : !WSL_GFORTRAN_BUG     end subroutine getEnergyFluence
     488             : 
     489             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     490             : 
     491             :     !> \brief
     492             :     !> Convert an **input energy fluence** in `[lowerLim, upperLim]` energy window, all in units of keV (or ALL in some other units),
     493             :     !> **to photon fluence**, within the same energy range, or if `[lowerLimNew, upperLimNew]` is provided, then in that range.
     494             :     !> The input `tolerance` is passed to integrators as a measure of the desired accuracy for the computation of `photonFluence`.
     495             :     !>
     496             :     !> \param[in]   energyFluence   :   The input energy fluence (in units of [keV]) to be converted to the output photon fluence.
     497             :     !> \param[in]   lowerLim        :   The lower limit energy (in units of [keV]) of the integration.
     498             :     !> \param[in]   upperLim        :   The upper limit energy (in units of [keV]) of the integration.
     499             :     !> \param[in]   epk             :   The spectral peak energy in units of [keV].
     500             :     !> \param[in]   alpha           :   The lower spectral exponent of the Band model.
     501             :     !> \param[in]   beta            :   The upper spectral exponent of the Band model.
     502             :     !> \param[in]   tolerance       :   The relative accuracy tolerance of the integration.
     503             :     !> \param[out]  photonFluence   :   The output fluence in units of photon counts.
     504             :     !> \param[out]  Err             :   An object of class [Err_type](@ref err_mod::err_type) containing error-handling information.
     505             :     !> \param[in]   lowerLimNew     :   The lower limit of energy windows (in keV) to be used for the computation of `photonFluence` (**optional**).
     506             :     !> \param[in]   upperLimNew     :   The upper limit of energy windows (in keV) to be used for the computation of `photonFluence` (**optional**).
     507             :     !>
     508             :     !> \remark
     509             :     !> If the optional `[lowerLimNew, upperLimNew]` are provided, each will replace the
     510             :     !> corresponding input `[lowerLim, upperLim]` in the computation of the output `photonFluence`.
     511          24 :     subroutine getPhotonFluenceFromEnergyFluence( energyFluence, lowerLim, upperLim, epk, alpha, beta, tolerance, photonFluence, Err, lowerLimNew, upperLimNew )
     512             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     513             :         !DEC$ ATTRIBUTES DLLEXPORT :: getPhotonFluenceFromEnergyFluence
     514             : #endif
     515           3 :         use Constants_mod, only: IK, RK, HUGE_RK
     516             :         use Err_mod, only: Err_type
     517             :         implicit none
     518             :         real(RK), intent(in)            :: energyFluence, lowerLim, upperLim, epk, alpha, beta, tolerance
     519             :         real(RK), intent(out)           :: photonFluence
     520             :         type(Err_type), intent(out)     :: Err
     521             :         character(*), parameter         :: PROCEDURE_NAME = "@getPhotonFluenceFromEnergyFluence()"
     522             :         real(RK), intent(in), optional  :: lowerLimNew, upperLimNew
     523          24 :         real(RK)                        :: lowLimNew, uppLimNew, amplitude
     524             : 
     525          24 :         Err%occurred = .false.
     526             : 
     527             :         ! check if the photon indices are consistent with the mathematical rules
     528             : 
     529          24 :         if (lowerLim>=upperLim) then
     530           3 :             photonFluence = 0._RK
     531           3 :             return
     532             :         end if
     533             : 
     534             :         ! check if the photon indices are consistent with the mathematical rules
     535          21 :         if (alpha<beta .or. alpha<-2._RK) then
     536           6 :             photonFluence = -HUGE_RK
     537           6 :             Err%occurred = .true.
     538           6 :             Err%msg = MODULE_NAME // PROCEDURE_NAME // ": Error occurred: alpha<beta .or. alpha<-2._RK"
     539           6 :             return
     540             :         end if
     541             : 
     542             : !#if defined OS_IS_WSL && defined CODECOV_ENABLED
     543             : !! LCOV_EXCL_START
     544             : !#endif
     545             : 
     546             :         lowLimNew = lowerLim
     547             :         if (present(lowerLimNew)) lowLimNew = lowerLimNew
     548             :         uppLimNew = upperLim
     549             :         if (present(upperLimNew)) uppLimNew = upperLimNew
     550             : 
     551             :         call getEnergyFluence(lowerLim,upperLim,epk,alpha,beta,tolerance,amplitude,Err)
     552             :         if (Err%occurred) then
     553             :         ! LCOV_EXCL_START
     554             :             photonFluence = -HUGE_RK
     555             :             Err%msg = MODULE_NAME // PROCEDURE_NAME // Err%msg
     556             :             return
     557             :         end if
     558             :         ! LCOV_EXCL_STOP
     559          15 :         amplitude = energyFluence / amplitude
     560             : 
     561          15 :         call getPhotonFluence(lowLimNew,uppLimNew,epk,alpha,beta,tolerance,photonFluence,Err)
     562          15 :         if (Err%occurred) then
     563             :         ! LCOV_EXCL_START
     564             :             photonFluence = -HUGE_RK
     565             :             Err%msg = MODULE_NAME // PROCEDURE_NAME // Err%msg
     566             :             return
     567             :         end if
     568             :         ! LCOV_EXCL_STOP
     569          15 :         photonFluence = amplitude * photonFluence
     570             : 
     571             : !#if defined OS_IS_WSL && defined CODECOV_ENABLED
     572             : !! LCOV_EXCL_STOP
     573             : !#endif
     574             : 
     575          24 :     end subroutine getPhotonFluenceFromEnergyFluence
     576             : 
     577             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     578             : 
     579             : end module BandSpectrum_mod ! LCOV_EXCL_LINE

ParaMonte: Plain Powerful Parallel Monte Carlo Library 
The Computational Data Science Lab
© Copyright 2012 - 2021