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 39 : 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 39 : real(RK) :: ebrk, alphaPlusTwo
39 39 : real(RK) :: thisUpperLim, betaPlusOne
40 39 : real(RK) :: alphaMinusBeta, coef
41 39 : real(RK) :: abserr
42 : integer(IK) :: neval
43 : integer(IK) :: ierr
44 : !real(RK) :: alphaPlusOne, logGammaAlphaPlusOne
45 :
46 39 : Err%occurred = .false.
47 :
48 39 : if (lowerLim>=upperLim) then
49 6 : photonFluence = 0._RK
50 6 : return
51 : end if
52 :
53 : ! check if the photon indices are consistent with the mathematical rules
54 33 : if (alpha<beta .or. alpha<-2._RK) then
55 6 : photonFluence = -HUGE_RK
56 6 : Err%occurred = .true.
57 6 : Err%msg = MODULE_NAME // PROCEDURE_NAME // ": Error occurred: alpha<beta .or. alpha<-2._RK"
58 6 : return
59 : end if
60 :
61 : ! integrate the spectrum
62 27 : alphaPlusTwo = alpha + 2._RK
63 27 : alphaMinusBeta = alpha - beta
64 27 : ebrk = epk*alphaMinusBeta/alphaPlusTwo
65 :
66 27 : if (lowerLim>ebrk) then
67 :
68 : ! there is only the high energy component in the photonFluence
69 6 : betaPlusOne = beta + 1._RK
70 6 : coef = ebrk**alphaMinusBeta * exp(-alphaMinusBeta);
71 6 : photonFluence = coef * ( upperLim**betaPlusOne - lowerLim**betaPlusOne ) / betaPlusOne
72 6 : return
73 :
74 21 : elseif (lowerLim<ebrk) then
75 :
76 21 : mv_alpha = alpha
77 :
78 21 : mv_alphaPlusTwoOverEpk = alphaPlusTwo / epk
79 21 : 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 21 : )
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 21 : 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 21 : if (upperLim>ebrk) then
130 : ! add the remaining part of the photonFluence from the high-energy component
131 6 : betaPlusOne = beta + 1._RK
132 6 : alphaMinusBeta = alpha - beta
133 6 : coef = ebrk**alphaMinusBeta * exp(-alphaMinusBeta);
134 6 : photonFluence = photonFluence + coef * ( upperLim**betaPlusOne - ebrk**betaPlusOne ) / betaPlusOne
135 6 : return
136 : end if
137 :
138 : end if
139 :
140 39 : end subroutine getPhotonFluence
141 :
142 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
143 :
144 5985 : 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 5985 : bandCompLow = energy**mv_alpha * exp(-mv_alphaPlusTwoOverEpk*energy)
152 6024 : end function getBandCompLowPhoton
153 :
154 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
155 :
156 : end submodule PhotonFluence_smod ! LCOV_EXCL_LINE
|
