Line data Source code
1 : submodule (BandSpectrum_mod) EnergyFluence_smod
2 :
3 : implicit none
4 :
5 : real(RK) :: mv_alphaPlusOne
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 in units of the input energy.
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] energyFluence : The fluence in units of the input energy (keV) 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 36 : module subroutine getEnergyFluence(lowerLim,upperLim,epk,alpha,beta,tolerance,energyFluence,Err)
26 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
27 : !DEC$ ATTRIBUTES DLLEXPORT :: getEnergyFluence
28 : #endif
29 :
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) :: energyFluence
36 : type(Err_type), intent(out) :: Err
37 : character(*), parameter :: PROCEDURE_NAME = "@getEnergyFluence()"
38 36 : real(RK) :: ebrk, alphaPlusTwo
39 36 : real(RK) :: thisUpperLim, betaPlusTwo
40 36 : real(RK) :: alphaMinusBeta, coef
41 36 : real(RK) :: abserr
42 : integer(IK) :: neval
43 : integer(IK) :: ierr
44 :
45 36 : Err%occurred = .false.
46 :
47 36 : if (lowerLim>=upperLim) then
48 3 : energyFluence = 0._RK
49 3 : return
50 : end if
51 :
52 : ! check if the photon indices are consistent with the mathematical rules
53 33 : if (alpha<beta .or. alpha<-2._RK) then
54 6 : energyFluence = -HUGE_RK
55 6 : Err%occurred = .true.
56 6 : Err%msg = MODULE_NAME // PROCEDURE_NAME // ": Error occurred: alpha<beta .or. alpha<-2._RK"
57 6 : return
58 : end if
59 :
60 : ! integrate the spectrum
61 27 : alphaPlusTwo = alpha + 2._RK
62 27 : alphaMinusBeta = alpha - beta
63 27 : ebrk = epk*alphaMinusBeta/alphaPlusTwo
64 :
65 27 : if (lowerLim>ebrk) then
66 :
67 : ! there is only the high energy component in the energyFluence
68 6 : betaPlusTwo = beta + 2._RK
69 6 : coef = ebrk**alphaMinusBeta * exp(-alphaMinusBeta);
70 6 : energyFluence = coef * ( upperLim**betaPlusTwo - lowerLim**betaPlusTwo ) / betaPlusTwo
71 6 : return
72 :
73 : !#if defined OS_IS_WSL && defined CODECOV_ENABLED
74 : !! LCOV_EXCL_START
75 : !#endif
76 :
77 : elseif (lowerLim<ebrk) then
78 :
79 : mv_alphaPlusTwoOverEpk = alphaPlusTwo / epk
80 : thisUpperLim = min(upperLim,ebrk)
81 : mv_alphaPlusOne = alpha + 1._RK
82 : ! use brute-force integration
83 : call qag( f = getBandCompLowEnergy &
84 : , a = lowerLim &
85 : , b = thisUpperLim &
86 : , epsabs = 0._RK &
87 : , epsrel = tolerance &
88 : , key = 1_IK &
89 : , result = energyFluence &
90 : , abserr = abserr &
91 : , neval = neval &
92 : , ier = ierr &
93 : )
94 :
95 : if (ierr/=0_IK) then
96 : energyFluence = -HUGE_RK
97 : Err%occurred = .true.
98 : Err%stat = ierr
99 : Err%msg = MODULE_NAME // PROCEDURE_NAME // ": Error occurred at QuadPack routine. Check the error code to identify the root cause."
100 : return
101 : end if
102 :
103 : if (upperLim>ebrk) then ! add the remaining part of the energyFluence from the high-energy component
104 : betaPlusTwo = beta + 2._RK
105 : alphaMinusBeta = alpha - beta
106 : coef = ebrk**alphaMinusBeta * exp(-alphaMinusBeta);
107 : energyFluence = energyFluence + coef * ( upperLim**betaPlusTwo - ebrk**betaPlusTwo ) / betaPlusTwo
108 : return
109 : end if
110 :
111 : end if
112 :
113 : end subroutine getEnergyFluence
114 :
115 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
116 :
117 : pure function getBandCompLowEnergy(energy) result(bandCompLow)
118 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
119 : !DEC$ ATTRIBUTES DLLEXPORT :: getBandCompLowEnergy
120 : #endif
121 : implicit none
122 : real(RK), intent(in) :: energy
123 : real(RK) :: bandCompLow
124 : bandCompLow = energy**mv_alphaPlusOne * exp(-mv_alphaPlusTwoOverEpk*energy)
125 : end function getBandCompLowEnergy
126 :
127 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
128 :
129 : end submodule EnergyFluence_smod ! LCOV_EXCL_LINE
|
