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 1 : 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 1 : ebrk = epk * (alpha-beta) / (alpha+2._RK)
123 2 : 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 1 : 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 1 : 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 1 : real(RK) :: alphaMinusBeta
147 1 : alphaPlusTwo = alpha + 2._RK
148 1 : alphaMinusBeta = alpha - beta
149 1 : ebrk = epk * alphaMinusBeta / alphaPlusTwo
150 1 : coef = ebrk**alphaMinusBeta * exp(-alphaMinusBeta)
151 1 : 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 4 : 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 1 : 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 4 : if (alpha<beta .or. alpha<-2._RK) then
185 2 : photonFlux = -HUGE_RK
186 2 : return
187 : end if
188 : ! compute the spectrum
189 2 : if (energy<=ebrk) then
190 1 : photonFlux = energy**alpha * exp(-energy*alphaPlusTwo/epk)
191 : else
192 1 : photonFlux = coef * energy**beta
193 : end if
194 4 : 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 1 : 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 4 : use Constants_mod, only: IK, RK, HUGE_RK
215 : implicit none
216 : real(RK), intent(in) :: energy, alpha, alphaPlusTwoOverEpk
217 : real(RK) :: photonFluxLower
218 1 : photonFluxLower = energy**alpha * exp(-energy*alphaPlusTwoOverEpk) ! the lower-energy spectrum
219 2 : 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 8 : 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 1 : 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 8 : real(RK) :: lowLimNew, uppLimNew, amplitude
524 :
525 8 : Err%occurred = .false.
526 :
527 : ! check if the photon indices are consistent with the mathematical rules
528 :
529 8 : if (lowerLim>=upperLim) then
530 1 : photonFluence = 0._RK
531 1 : return
532 : end if
533 :
534 : ! check if the photon indices are consistent with the mathematical rules
535 7 : if (alpha<beta .or. alpha<-2._RK) then
536 2 : photonFluence = -HUGE_RK
537 2 : Err%occurred = .true.
538 2 : Err%msg = MODULE_NAME // PROCEDURE_NAME // ": Error occurred: alpha<beta .or. alpha<-2._RK"
539 2 : 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 5 : amplitude = energyFluence / amplitude
560 :
561 5 : call getPhotonFluence(lowLimNew,uppLimNew,epk,alpha,beta,tolerance,photonFluence,Err)
562 5 : 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 5 : photonFluence = amplitude * photonFluence
570 :
571 : !#if defined OS_IS_WSL && defined CODECOV_ENABLED
572 : !! LCOV_EXCL_STOP
573 : !#endif
574 :
575 8 : end subroutine getPhotonFluenceFromEnergyFluence
576 :
577 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
578 :
579 : end module BandSpectrum_mod ! LCOV_EXCL_LINE
|
