The ParaMonte Documentation Website
Current view: top level - kernel - BandSpectrum_mod@PhotonFluence_smod.f90 (source / functions) Hit Total Coverage
Test: ParaMonte 1.5.1 :: Serial Kernel - Code Coverage Report Lines: 38 38 100.0 %
Date: 2021-01-08 13:03:42 Functions: 2 2 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : submodule (BandSpectrum_mod) PhotonFluence_smod
       2             : 
       3             :     implicit none
       4             : 
       5             :     real(RK)                        :: mv_alpha
       6             :     real(RK)                        :: mv_alphaPlusTwoOverEpk
       7             : 
       8             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       9             : 
      10             : contains
      11             : 
      12             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      13             : 
      14             :     !> \brief
      15             :     !> Integrate the Band differential spectrum over the input energy range.
      16             :     !>
      17             :     !> \param[in]   lowerLim        :   The lower limit energy (in units of [keV]) of the integration.
      18             :     !> \param[in]   upperLim        :   The upper limit energy (in units of [keV]) of the integration.
      19             :     !> \param[in]   epk             :   The spectral peak energy in units of [keV].
      20             :     !> \param[in]   alpha           :   The lower spectral exponent of the Band model.
      21             :     !> \param[in]   beta            :   The upper spectral exponent of the Band model.
      22             :     !> \param[in]   tolerance       :   The relative accuracy tolerance of the integration.
      23             :     !> \param[out]  photonFluence   :   The fluence in units of photon counts within the input energy range.
      24             :     !> \param[out]  Err             :   An object of class [Err_type](@ref err_mod::err_type) containing error-handling information.
      25          13 :     module subroutine getPhotonFluence(lowerLim,upperLim,epk,alpha,beta,tolerance,photonFluence,Err)
      26             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
      27             :         !DEC$ ATTRIBUTES DLLEXPORT :: getPhotonFluence
      28             : #endif
      29             :        !use Integration_mod, only: doQuadRombClosed
      30             :         use Constants_mod, only: IK, RK, HUGE_RK
      31             :         use QuadPackSPR_mod, only: qag
      32             :         use Err_mod, only: Err_type
      33             :         implicit none
      34             :         real(RK), intent(in)            :: lowerLim, upperLim, epk, alpha, beta, tolerance
      35             :         real(RK), intent(out)           :: photonFluence
      36             :         type(Err_type), intent(out)     :: Err
      37             :         character(*), parameter         :: PROCEDURE_NAME = "@getPhotonFluence()"
      38          13 :         real(RK)                        :: ebrk, alphaPlusTwo
      39          13 :         real(RK)                        :: thisUpperLim, betaPlusOne
      40          13 :         real(RK)                        :: alphaMinusBeta, coef
      41          13 :         real(RK)                        :: abserr
      42             :         integer(IK)                     :: neval
      43             :         integer(IK)                     :: ierr
      44             :        !real(RK)                        :: alphaPlusOne, logGammaAlphaPlusOne
      45             : 
      46          13 :         Err%occurred = .false.
      47             : 
      48          13 :         if (lowerLim>=upperLim) then
      49           2 :             photonFluence = 0._RK
      50           2 :             return
      51             :         end if
      52             : 
      53             :         ! check if the photon indices are consistent with the mathematical rules
      54          11 :         if (alpha<beta .or. alpha<-2._RK) then
      55           2 :             photonFluence = -HUGE_RK
      56           2 :             Err%occurred = .true.
      57           2 :             Err%msg = MODULE_NAME // PROCEDURE_NAME // ": Error occurred: alpha<beta .or. alpha<-2._RK"
      58           2 :             return
      59             :         end if
      60             : 
      61             :         ! integrate the spectrum
      62           9 :         alphaPlusTwo = alpha + 2._RK
      63           9 :         alphaMinusBeta = alpha - beta
      64           9 :         ebrk = epk*alphaMinusBeta/alphaPlusTwo
      65             : 
      66           9 :         if (lowerLim>ebrk) then
      67             : 
      68             :             ! there is only the high energy component in the photonFluence
      69           2 :             betaPlusOne = beta + 1._RK
      70           2 :             coef = ebrk**alphaMinusBeta * exp(-alphaMinusBeta);
      71           2 :             photonFluence = coef * ( upperLim**betaPlusOne - lowerLim**betaPlusOne ) / betaPlusOne
      72           2 :             return
      73             : 
      74           7 :         elseif (lowerLim<ebrk) then
      75             : 
      76           7 :             mv_alpha = alpha
      77             : 
      78           7 :             mv_alphaPlusTwoOverEpk = alphaPlusTwo / epk
      79           7 :             thisUpperLim = min(upperLim,ebrk)
      80             :             !alphaPlusOne = alpha + 1._RK
      81             :             !if (alpha>-1._RK) then
      82             :             !    logGammaAlphaPlusOne = log_gamma( alphaPlusOne )
      83             :             !    ! use the analytical approach to compute the photonFluence:
      84             :             !    ! https://www.wolframalpha.com/input/?i=integrate+x%5Ea+*+exp(-b*x)
      85             :             !    photonFluence = getUpperGamma( exponent = alphaPlusOne &
      86             :             !                            , logGammaExponent = logGammaAlphaPlusOne &
      87             :             !                            , lowerLim = mv_alphaPlusTwoOverEpk * lowerLim &
      88             :             !                            , tolerance = tolerance &
      89             :             !                            ) &
      90             :             !             - getUpperGamma( exponent = alphaPlusOne &
      91             :             !                            , logGammaExponent = logGammaAlphaPlusOne &
      92             :             !                            , lowerLim = mv_alphaPlusTwoOverEpk * thisUpperLim &
      93             :             !                            , tolerance = tolerance &
      94             :             !                            )
      95             :             !    photonFluence = photonFluence / mv_alphaPlusTwoOverEpk**alphaPlusOne
      96             :             !else
      97             :                 ! use brute-force integration
      98             :                 call qag( f             = getBandCompLowPhoton  &
      99             :                         , a             = lowerLim              &
     100             :                         , b             = thisUpperLim          &
     101             :                         , epsabs        = 0._RK                 &
     102             :                         , epsrel        = tolerance             &
     103             :                         , key           = 1_IK                  &
     104             :                         , result        = photonFluence         &
     105             :                         , abserr        = abserr                &
     106             :                         , neval         = neval                 &
     107             :                         , ier           = ierr                  &
     108           7 :                         )
     109             :                 !write(*,*) neval
     110             :                 !call doQuadRombClosed   ( getFunc       = getBandCompLowPhoton &
     111             :                 !                        , xmin          = lowerLim             &
     112             :                 !                        , xmax          = thisUpperLim         &
     113             :                 !                        , tolerance     = 1.e-7_RK             &
     114             :                 !                        , nRefinement   = 10_IK                &
     115             :                 !                        , photonFluence = photonFluence        &
     116             :                 !                        , ierr          = ierr                 &
     117             :                 !                        )
     118           7 :                 if (ierr/=0_IK) then
     119             :                 ! LCOV_EXCL_START
     120             :                     photonFluence = -HUGE_RK
     121             :                     Err%occurred = .true.
     122             :                     Err%stat = ierr
     123             :                     Err%msg = MODULE_NAME // PROCEDURE_NAME // ": Error occurred at QuadPack routine. Check the error code to identify the root cause."
     124             :                     return
     125             :                 ! LCOV_EXCL_STOP
     126             :                 end if
     127             :             !end if
     128             : 
     129           7 :             if (upperLim>ebrk) then
     130             :                 ! add the remaining part of the photonFluence from the high-energy component
     131           2 :                 betaPlusOne = beta + 1._RK
     132           2 :                 alphaMinusBeta = alpha - beta
     133           2 :                 coef = ebrk**alphaMinusBeta * exp(-alphaMinusBeta);
     134           2 :                 photonFluence = photonFluence + coef * ( upperLim**betaPlusOne - ebrk**betaPlusOne ) / betaPlusOne
     135           2 :                 return
     136             :             end if
     137             : 
     138             :         end if
     139             : 
     140          13 :     end subroutine getPhotonFluence
     141             : 
     142             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     143             : 
     144        1995 :     pure function getBandCompLowPhoton(energy) result(bandCompLow)
     145             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     146             :         !DEC$ ATTRIBUTES DLLEXPORT :: getBandCompLowPhoton
     147             : #endif
     148             :         implicit none
     149             :         real(RK), intent(in)    :: energy
     150             :         real(RK)                :: bandCompLow
     151        1995 :         bandCompLow = energy**mv_alpha * exp(-mv_alphaPlusTwoOverEpk*energy)
     152        2008 :     end function getBandCompLowPhoton
     153             : 
     154             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     155             : 
     156             : end submodule PhotonFluence_smod ! LCOV_EXCL_LINE

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