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 This module contains procedures for computing the Cosmic Star Formation Rate in the Universe.
44 : !> \author Amir Shahmoradi
45 :
46 : module StarFormation_mod
47 :
48 : use Constants_mod, only: IK, RK, PI, NEGINF_RK
49 :
50 : implicit none
51 :
52 : character(*), parameter :: MODULE_NAME = "@StarFormation_mod"
53 :
54 : real(RK) , parameter :: ZMIN = 0.1e0_RK
55 : real(RK) , parameter :: ZMAX = 1.0e2_RK
56 :
57 : abstract interface
58 : function getLogRate_proc(zplus1,logzplus1,twiceLogLumDisMpc) result(logRate)
59 : import :: RK, IK
60 : real(RK), intent(in) :: zplus1, logzplus1, twiceLogLumDisMpc
61 : real(RK) :: logRate
62 : end function getLogRate_proc
63 : end interface
64 :
65 : abstract interface
66 : function getRateDensity_proc(zplus1) result(rateDensity)
67 : import :: RK, IK
68 : real(RK), intent(in) :: zplus1
69 : real(RK) :: rateDensity
70 : end function getRateDensity_proc
71 : end interface
72 :
73 : abstract interface
74 : function getMergerDelayTimePDF_proc(mergerDelayTime) result(mergerDelayTimeProb)
75 : import :: RK, IK
76 : real(RK), intent(in) :: mergerDelayTime
77 : real(RK) :: mergerDelayTimeProb
78 : end function getMergerDelayTimePDF_proc
79 : end interface
80 :
81 : #if defined OS_IS_WSL
82 : procedure(getMergerDelayTimePDF_proc), pointer :: getMergerDelayTimePDF_WSL !< This madness bypasses the Microsoft Subsystem for Linux Internal Function call GFortran Segmentation Fault error.
83 : procedure(getRateDensity_proc), pointer :: getStarFormationRateDensity_WSL !< This madness bypasses the Microsoft Subsystem for Linux Internal Function call GFortran Segmentation Fault error.
84 : real(RK) :: maxRelativeErrorDefault_WSL !< This madness bypasses the Microsoft Subsystem for Linux Internal Function call GFortran Segmentation Fault error.
85 : integer(IK) :: nRefinementDefault_WSL !< This madness bypasses the Microsoft Subsystem for Linux Internal Function call GFortran Segmentation Fault error.
86 : real(RK) :: lookBackTimeRef_WSL !< This madness bypasses the Microsoft Subsystem for Linux Internal Function call GFortran Segmentation Fault error.
87 : #endif
88 :
89 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
90 :
91 : contains
92 :
93 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
94 :
95 : !> \brief
96 : !> Return GRBFR density based on the formation rate estimates of Petrosian et al (2015).
97 3 : pure function getLogRateDensityP15(logzplus1) result(logDensitySFR)
98 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
99 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogRateDensityP15
100 : #endif
101 : use Constants_mod, only: RK
102 : implicit none
103 : real(RK), intent(in) :: logzplus1
104 : real(RK), parameter :: logz1plus1 = log(1._RK+4.50_RK)
105 : real(RK), parameter :: exponentHighZ = -7.8_RK
106 : real(RK), parameter :: logNormFac2 = -exponentHighZ * logz1plus1
107 : real(RK) :: logDensitySFR
108 3 : if (logzplus1<0._RK) then
109 0 : logDensitySFR = NEGINF_RK
110 3 : elseif (logzplus1<logz1plus1) then
111 0 : logDensitySFR = 0._RK
112 : else
113 3 : logDensitySFR = logzplus1*exponentHighZ + logNormFac2
114 : end if
115 3 : end function getLogRateDensityP15
116 :
117 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
118 :
119 : !> \brief
120 : !> Return GRBFR density based on the formation rate estimates of Hopkins and Beacom (2007).
121 3 : pure function getLogRateDensityH06(logzplus1) result(logDensitySFR)
122 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
123 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogRateDensityH06
124 : #endif
125 3 : use Constants_mod, only: RK
126 : implicit none
127 : real(RK), intent(in) :: logzplus1
128 : real(RK), parameter :: logz0plus1 = log(1._RK+0.97_RK)
129 : real(RK), parameter :: logz1plus1 = log(1._RK+4.50_RK)
130 : real(RK), parameter :: g0 = +3.4_RK
131 : real(RK), parameter :: g1 = -0.3_RK
132 : real(RK), parameter :: g2 = -7.8_RK
133 : real(RK), parameter :: logNormFac1 = logz0plus1*(g0-g1)
134 : real(RK), parameter :: logNormFac2 = logz1plus1*(g1-g2) + logNormFac1
135 : real(RK) :: logDensitySFR
136 3 : if (logzplus1<0._RK) then
137 0 : logDensitySFR = NEGINF_RK
138 3 : elseif (logzplus1<logz0plus1) then
139 0 : logDensitySFR = logzplus1*g0
140 3 : elseif (logzplus1<logz1plus1) then
141 0 : logDensitySFR = logzplus1*g1 + logNormFac1
142 : else
143 3 : logDensitySFR = logzplus1*g2 + logNormFac2
144 : end if
145 3 : end function getLogRateDensityH06
146 :
147 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
148 :
149 : !> \brief
150 : !> Return GRBFR density based on the formation rate estimates of Li (2008).
151 3 : pure function getLogRateDensityL08(logzplus1) result(logDensitySFR)
152 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
153 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogRateDensityL08
154 : #endif
155 3 : use Constants_mod, only: RK
156 : implicit none
157 : real(RK), intent(in) :: logzplus1
158 : real(RK), parameter :: logz0plus1 = log(1._RK+0.993_RK)
159 : real(RK), parameter :: logz1plus1 = log(1._RK+3.800_RK)
160 : real(RK), parameter :: g0 = +3.3000_RK
161 : real(RK), parameter :: g1 = +0.0549_RK
162 : real(RK), parameter :: g2 = -4.4600_RK
163 : real(RK), parameter :: logNormFac1 = logz0plus1*(g0-g1)
164 : real(RK), parameter :: logNormFac2 = logz1plus1*(g1-g2) + logNormFac1
165 : real(RK) :: logDensitySFR
166 3 : if (logzplus1<0._RK) then
167 0 : logDensitySFR = NEGINF_RK
168 3 : elseif (logzplus1<logz0plus1) then
169 0 : logDensitySFR = logzplus1*g0
170 3 : elseif (logzplus1<logz1plus1) then
171 0 : logDensitySFR = logzplus1*g1 + logNormFac1
172 : else
173 3 : logDensitySFR = logzplus1*g2 + logNormFac2
174 : end if
175 3 : end function getLogRateDensityL08
176 :
177 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
178 :
179 : !> \brief
180 : !> Return GRBFR density based on the formation rate estimates of Butler et al. (2010).
181 : !>
182 : !> \remark
183 : !> The formation estimate of B10 takes into acount the metalicity-correction to the SFR.
184 3891 : pure function getLogRateDensityB10(logzplus1) result(logDensitySFR)
185 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
186 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogRateDensityB10
187 : #endif
188 3 : use Constants_mod, only: IK, RK
189 : implicit none
190 : real(RK), intent(in) :: logzplus1
191 : real(RK), parameter :: logz0plus1 = log(1._RK+0.97_RK)
192 : real(RK), parameter :: logz1plus1 = log(1._RK+4.00_RK)
193 : real(RK), parameter :: g0 = +3.14_RK
194 : real(RK), parameter :: g1 = +1.36_RK
195 : real(RK), parameter :: g2 = -2.92_RK
196 : real(RK), parameter :: logNormFac1 = logz0plus1*(g0-g1)
197 : real(RK), parameter :: logNormFac2 = logz1plus1*(g1-g2) + logNormFac1
198 : real(RK) :: logDensitySFR
199 3891 : if (logzplus1<0._RK) then
200 0 : logDensitySFR = NEGINF_RK
201 3891 : elseif (logzplus1<logz0plus1) then
202 0 : logDensitySFR = logzplus1*g0
203 3891 : elseif (logzplus1<logz1plus1) then
204 1194 : logDensitySFR = logzplus1*g1 + logNormFac1
205 : else
206 2697 : logDensitySFR = logzplus1*g2 + logNormFac2
207 : end if
208 3891 : end function getLogRateDensityB10
209 :
210 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
211 :
212 : !> \brief
213 : !> Return the Comoving Star Formation Rate Density according to Eqn 15 of Madau 2014: Cosmic Star-Formation History.\n
214 : !> `densitySFR(z) = 0.015 * (1+z)^2.7 / ( 1 + [(1+z)/2.9]^5.6 )`
215 3 : pure function getLogRateDensityM14(zplus1,logzplus1) result(logDensitySFR)
216 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
217 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogRateDensityM14
218 : #endif
219 3891 : use Constants_mod, only: RK
220 : implicit none
221 : real(RK), intent(in) :: zplus1, logzplus1
222 : real(RK) :: logDensitySFR
223 : real(RK), parameter :: logAmplitude = log(0.015_RK)
224 : real(RK), parameter :: lowerExp = 2.7_RK
225 : real(RK), parameter :: upperExp = 5.6_RK
226 : real(RK), parameter :: zplus1Break = 2.9_RK
227 : real(RK), parameter :: zplus1Coeff = 1._RK / (zplus1Break**upperExp)
228 3 : logDensitySFR = logAmplitude + lowerExp*logzplus1 - log( 1._RK + zplus1Coeff * zplus1**upperExp )
229 6 : end function getLogRateDensityM14
230 :
231 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
232 :
233 : !> \brief
234 : !> Return the Comoving Star Formation Rate Density according to Eqn 1 of Madau 2017: Cosmic Star-Formation History
235 : !> `densitySFR(z) = 0.01 * (1+z)^2.6 / ( 1 + [(1+z)/3.2]^6.2 )`
236 3 : pure function getLogRateDensityM17(zplus1,logzplus1) result(logDensitySFR)
237 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
238 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogRateDensityM17
239 : #endif
240 3 : use Constants_mod, only: RK
241 : implicit none
242 : real(RK), intent(in) :: zplus1, logzplus1
243 : real(RK) :: logDensitySFR
244 : real(RK), parameter :: logAmplitude = log(0.01_RK)
245 : real(RK), parameter :: lowerExp = 2.6_RK
246 : real(RK), parameter :: upperExp = 6.2_RK
247 : real(RK), parameter :: zplus1Break = 3.2_RK
248 : real(RK), parameter :: zplus1Coeff = 1._RK / (zplus1Break**upperExp)
249 3 : logDensitySFR = logAmplitude + lowerExp*logzplus1 - log( 1._RK + zplus1Coeff * zplus1**upperExp )
250 6 : end function getLogRateDensityM17
251 :
252 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
253 :
254 : !> \brief
255 : !> Return the Mordau Comoving Star Formation Rate Density with updated parameters from Fermi 2018.\n
256 : ! `densitySFR(z) = 0.013 * (1+z)^2.99 / ( 1 + [(1+z)/2.63]^6.19 )`
257 3 : pure function getLogRateDensityF18(zplus1,logzplus1) result(logDensitySFR)
258 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
259 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogRateDensityF18
260 : #endif
261 3 : use Constants_mod, only: RK
262 : implicit none
263 : real(RK), intent(in) :: zplus1, logzplus1
264 : real(RK) :: logDensitySFR
265 : real(RK), parameter :: logAmplitude = log(0.013_RK)
266 : real(RK), parameter :: lowerExp = 2.99_RK
267 : real(RK), parameter :: upperExp = 6.19_RK
268 : real(RK), parameter :: zplus1Break = 2.63_RK
269 : real(RK), parameter :: zplus1Coeff = 1._RK / (zplus1Break**upperExp)
270 3 : logDensitySFR = logAmplitude + lowerExp*logzplus1 - log( 1._RK + zplus1Coeff * zplus1**upperExp )
271 6 : end function getLogRateDensityF18
272 :
273 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
274 :
275 : !> \brief
276 : !> Return the cosmic formation rate according to the work of Petrosian et al. (2015).
277 3 : pure function getLogRateP15(zplus1,logzplus1,twiceLogLumDisMpc) result(logRate)
278 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
279 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogRateP15
280 : #endif
281 3 : use Cosmology_mod, only: LS2HC, OMEGA_DM, OMEGA_DE
282 : use Constants_mod, only: RK, PI
283 : implicit none
284 : real(RK), intent(in) :: zplus1, logzplus1, twiceLogLumDisMpc
285 : real(RK), parameter :: LOG_COEF = log(4._RK*PI*LS2HC)
286 : real(RK) :: logRate
287 : logRate = LOG_COEF + twiceLogLumDisMpc - ( 3._RK*logzplus1 + 0.5_RK*log(OMEGA_DM*zplus1**3+OMEGA_DE) ) &
288 3 : + getLogRateDensityP15(logzplus1)
289 6 : end function getLogRateP15
290 :
291 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
292 :
293 : !> \brief
294 : !> Return the cosmic star formation rate according to the work of Hopkins and Beacom (2007).
295 3 : pure function getLogRateH06(zplus1,logzplus1,twiceLogLumDisMpc) result(logRate)
296 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
297 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogRateH06
298 : #endif
299 3 : use Cosmology_mod, only: LS2HC, OMEGA_DM, OMEGA_DE
300 : use Constants_mod, only: RK, PI
301 : implicit none
302 : real(RK), intent(in) :: zplus1, logzplus1, twiceLogLumDisMpc
303 : real(RK), parameter :: LOG_COEF = log(4._RK*PI*LS2HC)
304 : real(RK) :: logRate
305 : logRate = LOG_COEF + twiceLogLumDisMpc - ( 3._RK*logzplus1 + 0.5_RK*log(OMEGA_DM*zplus1**3+OMEGA_DE) ) &
306 3 : + getLogRateDensityH06(logzplus1)
307 6 : end function getLogRateH06
308 :
309 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
310 :
311 : !> \brief
312 : !> Return the cosmic star formation rate according to the work of Li (2008).
313 3 : pure function getLogRateL08(zplus1,logzplus1,twiceLogLumDisMpc) result(logRate)
314 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
315 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogRateL08
316 : #endif
317 3 : use Cosmology_mod, only: LS2HC, OMEGA_DM, OMEGA_DE
318 : use Constants_mod, only: RK, PI
319 : implicit none
320 : real(RK), intent(in) :: zplus1, logzplus1, twiceLogLumDisMpc
321 : real(RK), parameter :: LOG_COEF = log(4._RK*PI*LS2HC)
322 : real(RK) :: logRate
323 : logRate = LOG_COEF + twiceLogLumDisMpc - ( 3._RK*logzplus1 + 0.5_RK*log(OMEGA_DM*zplus1**3+OMEGA_DE) ) &
324 3 : + getLogRateDensityL08(logzplus1)
325 6 : end function getLogRateL08
326 :
327 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
328 :
329 : !> \brief
330 : !> Return the cosmic GRB formation rate according to the work of Butler et al. (2010).
331 3 : pure function getLogRateB10(zplus1,logzplus1,twiceLogLumDisMpc) result(logRate)
332 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
333 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogRateB10
334 : #endif
335 3 : use Cosmology_mod, only: LS2HC, OMEGA_DM, OMEGA_DE
336 : use Constants_mod, only: RK, PI
337 : implicit none
338 : real(RK), intent(in) :: zplus1, logzplus1, twiceLogLumDisMpc
339 : real(RK), parameter :: LOG_COEF = log(4._RK*PI*LS2HC)
340 : real(RK) :: logRate
341 : logRate = LOG_COEF + twiceLogLumDisMpc - ( 3._RK*logzplus1 + 0.5_RK*log(OMEGA_DM*zplus1**3+OMEGA_DE) ) &
342 3 : + getLogRateDensityB10(logzplus1)
343 6 : end function getLogRateB10
344 :
345 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
346 :
347 : !> \brief
348 : !> Return the cosmic star formation rate according to the work of Madau et al. (2014).
349 3 : pure function getLogRateM14(zplus1,logzplus1,twiceLogLumDisMpc) result(logRate)
350 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
351 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogRateM14
352 : #endif
353 3 : use Cosmology_mod, only: LS2HC, OMEGA_DM, OMEGA_DE
354 : use Constants_mod, only: RK, PI
355 : implicit none
356 : real(RK), intent(in) :: zplus1, logzplus1, twiceLogLumDisMpc
357 : real(RK), parameter :: LOG_COEF = log(4._RK*PI*LS2HC)
358 : real(RK) :: logRate
359 : logRate = LOG_COEF + twiceLogLumDisMpc - ( 3._RK*logzplus1 + 0.5_RK*log(OMEGA_DM*zplus1**3+OMEGA_DE) ) &
360 3 : + getLogRateDensityM14(zplus1,logzplus1)
361 6 : end function getLogRateM14
362 :
363 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
364 :
365 : !> \brief
366 : !> Return the cosmic star formation rate according to the work of Madau et al. (2017).
367 3 : pure function getLogRateM17(zplus1,logzplus1,twiceLogLumDisMpc) result(logRate)
368 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
369 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogRateM17
370 : #endif
371 3 : use Cosmology_mod, only: LS2HC, OMEGA_DM, OMEGA_DE
372 : use Constants_mod, only: RK, PI
373 : implicit none
374 : real(RK), intent(in) :: zplus1, logzplus1, twiceLogLumDisMpc
375 : real(RK), parameter :: LOG_COEF = log(4._RK*PI*LS2HC)
376 : real(RK) :: logRate
377 : logRate = LOG_COEF + twiceLogLumDisMpc - ( 3._RK*logzplus1 + 0.5_RK*log(OMEGA_DM*zplus1**3+OMEGA_DE) ) &
378 3 : + getLogRateDensityM17(zplus1,logzplus1)
379 6 : end function getLogRateM17
380 :
381 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
382 :
383 : !> \brief
384 : !> Return the cosmic star formation rate according to the work of the Fermi collaboration (2018).
385 3 : pure function getLogRateF18(zplus1,logzplus1,twiceLogLumDisMpc) result(logRate)
386 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
387 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogRateF18
388 : #endif
389 3 : use Cosmology_mod, only: LS2HC, OMEGA_DM, OMEGA_DE
390 : use Constants_mod, only: RK, PI
391 : implicit none
392 : real(RK), intent(in) :: zplus1, logzplus1, twiceLogLumDisMpc
393 : real(RK), parameter :: LOG_COEF = log(4._RK*PI*LS2HC)
394 : real(RK) :: logRate
395 : logRate = LOG_COEF + twiceLogLumDisMpc - ( 3._RK*logzplus1 + 0.5_RK*log(OMEGA_DM*zplus1**3+OMEGA_DE) ) &
396 3 : + getLogRateDensityF18(zplus1,logzplus1)
397 6 : end function getLogRateF18
398 :
399 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
400 :
401 : ! Amir Shahmoradi, Wednesday 5:43 PM, December 25, 2013, IFS, UT Austin
402 : ! compute the binary merger rate as a function of z, according to Shahmoradi and Nemiroff (2015)
403 : ! equivalent to delayed_rate_Belz_Li(z) in S15
404 : ! returns 0, if z>6.501_RK or z<0.09_RK
405 12 : pure function getBinaryMergerRateS15(z) result(binaryMergerRateS15)
406 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
407 : !DEC$ ATTRIBUTES DLLEXPORT :: getBinaryMergerRateS15
408 : #endif
409 3 : use Cosmology_mod, only: LS2HC, OMEGA_DM, OMEGA_DE
410 : use Constants_mod, only: RK, PI
411 : implicit none
412 : real(RK) , intent(in) :: z
413 : real(RK) :: binaryMergerRateS15
414 12 : if (z>2.5_RK .and. z<=6.501_RK) then
415 : binaryMergerRateS15 = &
416 : - 2.09118024744342000_RK &
417 : + 5.15382361299299000_RK * z &
418 : - 5.46442271664195000_RK * z**2 &
419 : + 3.29445310883082000_RK * z**3 &
420 : - 1.24547016168265000_RK * z**4 &
421 : + 0.30628893690508400_RK * z**5 &
422 : - 0.04904403249641820_RK * z**6 &
423 : + 0.00493757380504717_RK * z**7 &
424 : - 2.84061971928750e-4_RK * z**8 &
425 3 : + 7.12674138757750e-6_RK * z**9
426 9 : elseif (z>1.0_RK .and. z<=2.5_RK) then
427 : binaryMergerRateS15 = &
428 : - 0.86022576265904100_RK &
429 : + 4.22669545558817000_RK * z &
430 : - 8.86086728534670000_RK * z**2 &
431 : + 10.4863792284648000_RK * z**3 &
432 : - 7.64722909221129000_RK * z**4 &
433 : + 3.51616699500767000_RK * z**5 &
434 : - 0.99555474471022000_RK * z**6 &
435 : + 0.15876893754371900_RK * z**7 &
436 3 : - 0.01092541997736420_RK * z**8
437 6 : elseif (z<=1._RK .and. z>=0.09_RK) then
438 : binaryMergerRateS15 = &
439 : + 1.92595299989370e-4_RK &
440 : - 0.00345273599582578_RK * z &
441 : + 0.03157500615320920_RK * z**2 &
442 : - 0.04470545521198460_RK * z**3 &
443 : + 0.06812481521281660_RK * z**4 &
444 0 : - 0.03846033416253570_RK * z**5
445 : else
446 6 : binaryMergerRateS15 = 0._RK
447 : end if
448 12 : end function getBinaryMergerRateS15
449 :
450 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
451 :
452 : ! computes the natural log of the binary merger rate as a function of logzplus1.
453 : ! However, note that the computed rate is dN/dz, even though the input is logzplus1.
454 : ! the merger delay time distribution is the same as that of Shahmoradi and Nemiroff (2015), but with the logMean:log(0.1_RK) and sigma: 0.9612813_RK..
455 : ! returns 0, if z>19.929999999999882_RK or z<0.03_RK
456 21 : pure function getLogBinaryMergerRateLognormH06(logzplus1) result(logBinaryMergerRate)
457 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
458 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogBinaryMergerRateLognormH06
459 : #endif
460 12 : use Cosmology_mod, only: LS2HC, OMEGA_DM, OMEGA_DE
461 : use Constants_mod, only: RK, PI
462 : implicit none
463 : real(RK) , intent(in) :: logzplus1
464 : real(RK) :: logBinaryMergerRate
465 21 : if (logzplus1>0.02955880224154443_RK .and. logzplus1<=0.1441003439737565_RK) then
466 : logBinaryMergerRate = &
467 : - 14.26464149493092000_RK &
468 : + 84.73477757043948000_RK * logzplus1 &
469 : - 488.5893985602366500_RK * logzplus1**2 &
470 3 : + 1154.414655194473900_RK * logzplus1**3
471 18 : elseif (logzplus1>0.1441003439737565_RK .and. logzplus1<=0.6575200029167926_RK) then
472 : logBinaryMergerRate = &
473 : - 11.19700066921606300_RK &
474 : + 20.46712963401572300_RK * logzplus1 &
475 : - 24.31794334813894300_RK * logzplus1**2 &
476 3 : + 12.21213317590724400_RK * logzplus1**3
477 15 : elseif (logzplus1>0.6575200029167926_RK .and. logzplus1<=1.5591966959973538_RK) then
478 : logBinaryMergerRate = &
479 : - 9.094912666461765000_RK &
480 : + 15.23119806754538900_RK * logzplus1 &
481 : - 18.77526325204311800_RK * logzplus1**2 &
482 : + 9.941360355936961000_RK * logzplus1**3 &
483 3 : - 2.077370913197473000_RK * logzplus1**4
484 12 : elseif (logzplus1>1.5591966959973538_RK .and. logzplus1<=1.7056567701746455_RK) then
485 : logBinaryMergerRate = &
486 : - 2392.907733171019000_RK &
487 : + 6210.872744126407000_RK * logzplus1 &
488 : - 6054.866136454215000_RK * logzplus1**2 &
489 : + 2622.628785434413700_RK * logzplus1**3 &
490 3 : - 426.0273477222719000_RK * logzplus1**4
491 9 : elseif (logzplus1>1.7056567701746455_RK .and. logzplus1<=3.0411835364579027_RK) then
492 : logBinaryMergerRate = &
493 : + 9.538876239886940000_RK &
494 : - 8.753418172517534000_RK * logzplus1 &
495 : - 0.159980818030374640_RK * logzplus1**2 &
496 3 : - 0.088551503657680930_RK * logzplus1**3
497 6 : elseif (logzplus1<=0.02955880224154443_RK .or. logzplus1>3.0411835364579027_RK) then
498 6 : logBinaryMergerRate = 0._RK
499 : end if
500 21 : end function getLogBinaryMergerRateLognormH06
501 :
502 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
503 :
504 : ! computes the natural log of the binary merger rate as a function of logzplus1.
505 : ! However, note that the computed rate is dN/dz, even though the input is logzplus1.
506 : ! the merger delay time distribution is the same as that of Shahmoradi and Nemiroff (2015), but with the logMean:log(0.1_RK) and sigma: 0.9612813_RK..
507 : ! returns 0, if z>19.929999999999882_RK or z<0.03_RK
508 21 : pure function getLogBinaryMergerRateLognormL08(logzplus1) result(logBinaryMergerRate)
509 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
510 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogBinaryMergerRateLognormL08
511 : #endif
512 21 : use Cosmology_mod, only: LS2HC, OMEGA_DM, OMEGA_DE
513 : use Constants_mod, only: RK, PI
514 : implicit none
515 : real(RK) , intent(in) :: logzplus1
516 : real(RK) :: logBinaryMergerRate
517 21 : if (logzplus1>0.02955880224154443_RK .and. logzplus1<=0.20701416938432557_RK) then
518 : logBinaryMergerRate = &
519 : - 14.53696144309043900_RK &
520 : + 94.70274747509626000_RK * logzplus1 &
521 : - 687.3663996060040000_RK * logzplus1**2 &
522 : + 2695.421036673770700_RK * logzplus1**3 &
523 3 : - 4077.601561165490000_RK * logzplus1**4
524 18 : elseif (logzplus1>0.20701416938432557_RK .and. logzplus1<=0.8241754429663476_RK) then
525 : logBinaryMergerRate = &
526 : - 13.5104005566057670_RK &
527 : + 49.6443928683743600_RK * logzplus1 &
528 : - 164.286063097338630_RK * logzplus1**2 &
529 : + 315.721394966368100_RK * logzplus1**3 &
530 : - 300.345052726248640_RK * logzplus1**4 &
531 3 : + 108.470535327547080_RK * logzplus1**5
532 15 : elseif (logzplus1>0.8241754429663476_RK .and. logzplus1<=1.4243124283074096_RK) then
533 : logBinaryMergerRate = &
534 : - 8.77634469738400500_RK &
535 : + 13.1999684738558810_RK * logzplus1 &
536 : - 15.8698236818922140_RK * logzplus1**2 &
537 : + 8.48676936452957000_RK * logzplus1**3 &
538 3 : - 1.83190451512279620_RK * logzplus1**4
539 12 : elseif (logzplus1>1.4243124283074096_RK .and. logzplus1<=1.6154199841116488_RK) then
540 : logBinaryMergerRate = &
541 : + 4158.29353781047900_RK &
542 : - 10954.1105856433040_RK * logzplus1 &
543 : + 10789.3451136201870_RK * logzplus1**2 &
544 : - 4713.80244702217800_RK * logzplus1**3 &
545 3 : + 770.488645040204600_RK * logzplus1**4
546 9 : elseif (logzplus1>1.6154199841116488_RK .and. logzplus1<=3.0411835364579027_RK) then
547 : logBinaryMergerRate = &
548 : + 0.37742655174185624_RK &
549 : + 0.30883738015163340_RK * logzplus1 &
550 : - 4.04937550957291800_RK * logzplus1**2 &
551 : + 1.11680537027038170_RK * logzplus1**3 &
552 3 : - 0.13770838345089523_RK * logzplus1**4
553 6 : elseif (logzplus1<=0.02955880224154443_RK .or. logzplus1>3.0411835364579027_RK) then
554 6 : logBinaryMergerRate = 0._RK
555 : end if
556 21 : end function getLogBinaryMergerRateLognormL08
557 :
558 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
559 :
560 : ! computes the natural log of the binary merger rate as a function of logzplus1.
561 : ! However, note that the computed rate is dN/dz, even though the input is logzplus1.
562 : ! the merger delay time distribution is the same as that of Shahmoradi and Nemiroff (2015), but with the logMean:log(0.1_RK) and sigma: 0.9612813_RK.
563 : ! returns 0, if z>19.929999999999882_RK or z<0.03_RK
564 21 : pure function getLogBinaryMergerRateLognormB10(logzplus1) result(logBinaryMergerRate)
565 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
566 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogBinaryMergerRateLognormB10
567 : #endif
568 21 : use Cosmology_mod, only: LS2HC, OMEGA_DM, OMEGA_DE
569 : use Constants_mod, only: RK, PI
570 : implicit none
571 : real(RK) , intent(in) :: logzplus1
572 : real(RK) :: logBinaryMergerRate
573 21 : if (logzplus1>0.02955880224154443_RK .and. logzplus1<=0.20701416938432557_RK) then
574 : logBinaryMergerRate = &
575 : - 15.27802857671202000_RK &
576 : + 94.54179164991284000_RK * logzplus1 &
577 : - 687.3676159275769000_RK * logzplus1**2 &
578 : + 2695.420977251770600_RK * logzplus1**3 &
579 3 : - 4077.601406650646000_RK * logzplus1**4
580 18 : elseif (logzplus1>0.20701416938432557_RK .and. logzplus1<=0.8241754429663476_RK) then
581 : logBinaryMergerRate = &
582 : - 13.5066182170954650_RK &
583 : + 40.1985219822299200_RK * logzplus1 &
584 : - 121.506350703598660_RK * logzplus1**2 &
585 : + 224.621285123736100_RK * logzplus1**3 &
586 : - 210.878836655472500_RK * logzplus1**4 &
587 3 : + 76.3335749498628400_RK * logzplus1**5
588 15 : elseif (logzplus1>0.8241754429663476_RK .and. logzplus1<=1.4243124283074096_RK) then
589 : logBinaryMergerRate = &
590 : - 10.0515447816113700_RK &
591 : + 12.6659826494097970_RK * logzplus1 &
592 : - 13.2268991886238200_RK * logzplus1**2 &
593 : + 6.84523627043807100_RK * logzplus1**3 &
594 3 : - 1.44645280124922220_RK * logzplus1**4
595 12 : elseif (logzplus1>1.4243124283074096_RK .and. logzplus1<=1.6104374127671848_RK) then
596 : logBinaryMergerRate = &
597 : - 1187.90539057029950_RK &
598 : + 3240.19327021926350_RK * logzplus1 &
599 : - 3330.70645904271000_RK * logzplus1**2 &
600 : + 1522.87499612399850_RK * logzplus1**3 &
601 3 : - 261.341408956542300_RK * logzplus1**4
602 9 : elseif (logzplus1>1.6104374127671848_RK .and. logzplus1<=3.0411835364579027_RK) then
603 : logBinaryMergerRate = &
604 : - 1.43934839576471260_RK &
605 : + 1.72951867017028120_RK * logzplus1 &
606 : - 4.06729555225025000_RK * logzplus1**2 &
607 : + 1.18253386764330200_RK * logzplus1**3 &
608 3 : - 0.15201156018584210_RK * logzplus1**4
609 6 : elseif (logzplus1<=0.02955880224154443_RK .or. logzplus1>3.0411835364579027_RK) then
610 6 : logBinaryMergerRate = 0._RK
611 : end if
612 21 : end function getLogBinaryMergerRateLognormB10
613 :
614 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
615 :
616 : ! computes the natural log of the binary merger rate as a function of logzplus1.
617 : ! However, note that the computed rate is dN/dz, even though the input is logzplus1.
618 : ! the merger delay time distribution is the same as that of Shahmoradi and Nemiroff (2015), but with the logMean:log(0.1_RK) and sigma: 0.9612813_RK..
619 : ! returns 0, if z>19.929999999999882_RK or z<0.03_RK
620 21 : pure function getLogBinaryMergerRateLognormM14(logzplus1) result(logBinaryMergerRate)
621 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
622 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogBinaryMergerRateLognormM14
623 : #endif
624 21 : use Cosmology_mod, only: LS2HC, OMEGA_DM, OMEGA_DE
625 : use Constants_mod, only: RK, PI
626 : implicit none
627 : real(RK) , intent(in) :: logzplus1
628 : real(RK) :: logBinaryMergerRate
629 21 : if (logzplus1>0.02955880224154443_RK .and. logzplus1<=0.16551443847757297_RK) then
630 : logBinaryMergerRate = &
631 : - 13.91129314580349600_RK &
632 : + 78.88963489621422000_RK * logzplus1 &
633 : - 420.9801740859396700_RK * logzplus1**2 &
634 3 : + 902.4783078800951000_RK * logzplus1**3
635 18 : elseif (logzplus1>0.16551443847757297_RK .and. logzplus1<=0.9282193027394269_RK) then
636 : logBinaryMergerRate = &
637 : - 11.00951046136480500_RK &
638 : + 21.38817515748999000_RK * logzplus1 &
639 : - 33.29451048508970000_RK * logzplus1**2 &
640 : + 29.32158835244860400_RK * logzplus1**3 &
641 3 : - 10.97377449040440000_RK * logzplus1**4
642 15 : elseif (logzplus1>0.9282193027394269_RK .and. logzplus1<=1.3937663759585892_RK) then
643 : logBinaryMergerRate = &
644 : - 8.254476015464371000_RK &
645 : + 3.620963332444886000_RK * logzplus1 &
646 : + 6.734585433384001000_RK * logzplus1**2 &
647 : - 9.151412394211048000_RK * logzplus1**3 &
648 3 : + 2.516171777428496000_RK * logzplus1**4
649 12 : elseif (logzplus1>1.3937663759585892_RK .and. logzplus1<=3.0411835364579027_RK) then
650 : logBinaryMergerRate = &
651 : - 6.539697727782377000_RK &
652 : + 8.522331572601950000_RK * logzplus1 &
653 : - 8.242990979412244000_RK * logzplus1**2 &
654 : + 2.316632169715435300_RK * logzplus1**3 &
655 6 : - 0.266462340853027450_RK * logzplus1**4
656 6 : elseif (logzplus1<=0.02955880224154443_RK .or. logzplus1>3.0411835364579027_RK) then
657 6 : logBinaryMergerRate = 0._RK
658 : end if
659 21 : end function getLogBinaryMergerRateLognormM14
660 :
661 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
662 :
663 : ! computes the natural log of the binary merger rate as a function of logzplus1.
664 : ! However, note that the computed rate is dN/dz, even though the input is logzplus1.
665 : ! the merger delay time distribution is the same as that of Shahmoradi and Nemiroff (2015), but with the logMean:log(0.1_RK) and sigma: 0.9612813_RK..
666 : ! returns 0, if z>19.929999999999882_RK or z<0.03_RK
667 21 : pure function getLogBinaryMergerRateLognormM17(logzplus1) result(logBinaryMergerRate)
668 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
669 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogBinaryMergerRateLognormM17
670 : #endif
671 21 : use Cosmology_mod, only: LS2HC, OMEGA_DM, OMEGA_DE
672 : use Constants_mod, only: RK, PI
673 : implicit none
674 : real(RK) , intent(in) :: logzplus1
675 : real(RK) :: logBinaryMergerRate
676 21 : if (logzplus1>0.02955880224154443_RK .and. logzplus1<=0.16551443847757297_RK) then
677 : logBinaryMergerRate = &
678 : - 14.01939141013502300_RK &
679 : + 78.80010843737509000_RK * logzplus1 &
680 : - 420.9593253775164000_RK * logzplus1**2 &
681 3 : + 902.5668042795056000_RK * logzplus1**3
682 18 : elseif (logzplus1>0.16551443847757297_RK .and. logzplus1<=0.9282193027394269_RK) then
683 : logBinaryMergerRate = &
684 : - 11.12915953695671500_RK &
685 : + 21.43217730985805500_RK * logzplus1 &
686 : - 33.75904206577289000_RK * logzplus1**2 &
687 : + 30.03916282499634200_RK * logzplus1**3 &
688 3 : - 11.12086545981264500_RK * logzplus1**4
689 15 : elseif (logzplus1>0.9282193027394269_RK .and. logzplus1<=1.3937663759585892_RK) then
690 : logBinaryMergerRate = &
691 : - 1.802362223155230800_RK &
692 : - 20.58526172567768200_RK * logzplus1 &
693 : + 38.93828966743146000_RK * logzplus1**2 &
694 : - 27.19863916580484500_RK * logzplus1**3 &
695 3 : + 6.138928143113263000_RK * logzplus1**4
696 12 : elseif (logzplus1>1.3937663759585892_RK .and. logzplus1<=3.0411835364579027_RK) then
697 : logBinaryMergerRate = &
698 : - 7.711815956299844000_RK &
699 : + 11.68891993486079700_RK * logzplus1 &
700 : - 10.62908897824095300_RK * logzplus1**2 &
701 : + 2.945685425778300700_RK * logzplus1**3 &
702 6 : - 0.327069839977957850_RK * logzplus1**4
703 6 : elseif (logzplus1<=0.02955880224154443_RK .or. logzplus1>3.0411835364579027_RK) then
704 6 : logBinaryMergerRate = 0._RK
705 : end if
706 21 : end function getLogBinaryMergerRateLognormM17
707 :
708 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
709 :
710 : ! computes the natural log of the binary merger rate as a function of logzplus1.
711 : ! However, note that the computed rate is dN/dz, even though the input is logzplus1.
712 : ! the merger delay time distribution is the same as that of Shahmoradi and Nemiroff (2015), but with the logMean:log(0.1_RK) and sigma: 0.9612813_RK..
713 : ! returns 0, if z>19.929999999999882_RK or z<0.03_RK
714 21 : pure function getLogBinaryMergerRateLognormF18(logzplus1) result(logBinaryMergerRate)
715 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
716 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogBinaryMergerRateLognormF18
717 : #endif
718 21 : use Cosmology_mod, only: LS2HC, OMEGA_DM, OMEGA_DE
719 : use Constants_mod, only: RK, PI
720 : implicit none
721 : real(RK) , intent(in) :: logzplus1
722 : real(RK) :: logBinaryMergerRate
723 21 : if (logzplus1>0.02955880224154443_RK .and. logzplus1<=0.16551443847757297_RK) then
724 : logBinaryMergerRate = &
725 : - 13.80128475140318000_RK &
726 : + 79.17963739241087000_RK * logzplus1 &
727 : - 420.9808813943490700_RK * logzplus1**2 &
728 3 : + 902.4149755380632000_RK * logzplus1**3
729 18 : elseif (logzplus1>0.16551443847757297_RK .and. logzplus1<=0.9282193027394269_RK) then
730 : logBinaryMergerRate = &
731 : - 10.89131941401880800_RK &
732 : + 21.59483273763076400_RK * logzplus1 &
733 : - 33.07662054750123000_RK * logzplus1**2 &
734 : + 29.23608723915102600_RK * logzplus1**3 &
735 3 : - 11.33984422193848700_RK * logzplus1**4
736 15 : elseif (logzplus1>0.9282193027394269_RK .and. logzplus1<=1.3937663759585892_RK) then
737 : logBinaryMergerRate = &
738 : - 14.02518830695731000_RK &
739 : + 24.92009885816913300_RK * logzplus1 &
740 : - 20.04762612951797000_RK * logzplus1**2 &
741 : + 4.885281389900379500_RK * logzplus1**3 &
742 3 : - 0.168402813838908260_RK * logzplus1**4
743 12 : elseif (logzplus1>1.3937663759585892_RK .and. logzplus1<=3.0411835364579027_RK) then
744 : logBinaryMergerRate = &
745 : - 4.348081430972656000_RK &
746 : + 4.815143234949144000_RK * logzplus1 &
747 : - 6.143880845780776000_RK * logzplus1**2 &
748 : + 1.738835623950871300_RK * logzplus1**3 &
749 6 : - 0.206972882929076480_RK * logzplus1**4
750 6 : elseif (logzplus1<=0.02955880224154443_RK .or. logzplus1>3.0411835364579027_RK) then
751 6 : logBinaryMergerRate = 0._RK
752 : end if
753 21 : end function getLogBinaryMergerRateLognormF18
754 :
755 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
756 :
757 24 : function getBinaryMergerRateDensity ( zplus1 &
758 : , zplus1Max &
759 : , nRefinement &
760 : , maxRelativeError &
761 : , getMergerDelayTimePDF &
762 : , getStarFormationRateDensity &
763 24 : ) result(binaryMergerRateDensity)
764 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
765 : !DEC$ ATTRIBUTES DLLEXPORT :: getBinaryMergerRateDensity
766 : #endif
767 : use, intrinsic :: iso_fortran_env, only: output_unit
768 21 : use Cosmology_mod, only: getLookBackTime
769 : use Constants_mod, only: RK, HUGE_RK
770 : use Integration_mod, only: doQuadRombOpen, ErrorMessage!, midinf
771 : use Integration_mod, only: midexp
772 : !use Integration_mod, only: midinf
773 : implicit none
774 : real(RK) , intent(in) :: zplus1
775 : real(RK) , intent(in), optional :: zplus1Max, maxRelativeError
776 : integer(IK) , intent(in), optional :: nRefinement
777 : procedure(getRateDensity_proc) :: getStarFormationRateDensity
778 : procedure(getMergerDelayTimePDF_proc) :: getMergerDelayTimePDF
779 : integer(IK) :: neval, ierr, nRefinementDefault
780 24 : real(RK) :: zplus1MaxDefault, maxRelativeErrorDefault, relerr
781 24 : real(RK) :: binaryMergerRateDensity, lookBackTimeRef
782 :
783 21 : nRefinementDefault = 7_IK; if (present(nRefinement)) nRefinementDefault = nRefinement
784 24 : zplus1MaxDefault = HUGE_RK; if (present(zplus1Max)) zplus1MaxDefault = zplus1Max
785 24 : maxRelativeErrorDefault = 1.e-6_RK; if (present(maxRelativeError)) maxRelativeErrorDefault = maxRelativeError
786 :
787 : lookBackTimeRef = getLookBackTime ( zplus1 = zplus1 &
788 : , maxRelativeError = maxRelativeErrorDefault &
789 : , nRefinement = nRefinementDefault &
790 24 : )
791 :
792 : #if defined OS_IS_WSL
793 : getMergerDelayTimePDF_WSL => getMergerDelayTimePDF
794 : getStarFormationRateDensity_WSL => getStarFormationRateDensity
795 : maxRelativeErrorDefault_WSL = maxRelativeErrorDefault
796 : nRefinementDefault_WSL = nRefinementDefault
797 : lookBackTimeRef_WSL = lookBackTimeRef
798 : #endif
799 :
800 : call doQuadRombOpen ( getFunc = getBinaryMergerRateDensityIntegrand &
801 : , integrate = midexp &
802 : !, integrate = midinf &
803 : , lowerLim = zplus1 &
804 : , upperLim = zplus1MaxDefault &
805 : , maxRelativeError = maxRelativeErrorDefault &
806 : , nRefinement = nRefinementDefault &
807 : , integral = binaryMergerRateDensity &
808 : , relativeError = relerr &
809 : , numFuncEval = neval &
810 : , ierr = ierr &
811 24 : )
812 24 : if (ierr/=0_IK) then
813 : ! LCOV_EXCL_START
814 : write(output_unit,"(A)") ErrorMessage(ierr)
815 : error stop
816 : ! LCOV_EXCL_STOP
817 : end if
818 :
819 : #if defined OS_IS_WSL
820 : nullify(getStarFormationRateDensity_WSL)
821 : nullify(getMergerDelayTimePDF_WSL)
822 : #else
823 : contains
824 :
825 3888 : function getBinaryMergerRateDensityIntegrand(zplus1) result(binaryMergerRateIntegrand)
826 :
827 24 : use Cosmology_mod, only: getUniverseAgeDerivative
828 : implicit none
829 : real(RK) , intent(in) :: zplus1
830 : real(RK) :: binaryMergerRateIntegrand !,lognormpdf
831 3888 : real(RK) :: mergerDelayTime
832 :
833 : ! note that zp<z always, so that delay>0.
834 : mergerDelayTime = getLookBackTime ( zplus1 = zplus1 &
835 : , maxRelativeError = maxRelativeErrorDefault &
836 : , nRefinement = nRefinementDefault &
837 3888 : )
838 3888 : mergerDelayTime = mergerDelayTime - lookBackTimeRef
839 3888 : if (mergerDelayTime<=0._RK) then
840 : ! LCOV_EXCL_START
841 : write(output_unit,"(A)") "The mergerDelayTime is non-positive in getBinaryMergerRateDensityIntegrand(): (zplus1, mergerDelayTime) = ", zplus1, mergerDelayTime
842 : error stop
843 : ! LCOV_EXCL_STOP
844 : end if
845 :
846 : binaryMergerRateIntegrand = getMergerDelayTimePDF(mergerDelayTime) &
847 : * getStarFormationRateDensity(zplus1) &
848 3888 : * getUniverseAgeDerivative(zplus1)
849 :
850 3888 : end function getBinaryMergerRateDensityIntegrand
851 : #endif
852 : end function getBinaryMergerRateDensity
853 :
854 : #if defined OS_IS_WSL
855 : ! This madness bypasses the Microsoft Subsystem for Linux Internal Function call GFortran Segmentation Fault error.
856 : function getBinaryMergerRateDensityIntegrand(zplus1) result(binaryMergerRateIntegrand)
857 : use, intrinsic :: iso_fortran_env, only: output_unit
858 : use Cosmology_mod, only: getUniverseAgeDerivative
859 : use Cosmology_mod, only: getLookBackTime
860 : implicit none
861 : real(RK) , intent(in) :: zplus1
862 : real(RK) :: binaryMergerRateIntegrand !,lognormpdf
863 : real(RK) :: mergerDelayTime
864 : ! note that zp<z always, so that delay>0.
865 : mergerDelayTime = getLookBackTime ( zplus1 = zplus1 &
866 : , maxRelativeError = maxRelativeErrorDefault_WSL &
867 : , nRefinement = nRefinementDefault_WSL &
868 : )
869 : mergerDelayTime = mergerDelayTime - lookBackTimeRef_WSL
870 : if (mergerDelayTime<=0._RK) then
871 : ! LCOV_EXCL_START
872 : write(output_unit,"(A)") "The mergerDelayTime is non-positive in getBinaryMergerRateDensityIntegrand(): (zplus1, mergerDelayTime) = ", zplus1, mergerDelayTime
873 : error stop
874 : end if
875 : ! LCOV_EXCL_START
876 :
877 : binaryMergerRateIntegrand = getMergerDelayTimePDF_WSL(mergerDelayTime) &
878 : * getStarFormationRateDensity_WSL(zplus1) &
879 : * getUniverseAgeDerivative(zplus1)
880 :
881 : end function getBinaryMergerRateDensityIntegrand
882 : #endif
883 :
884 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
885 :
886 : function getBinaryMergerRate( zplus1 &
887 : , zplus1Max &
888 : , nRefinement &
889 : , maxRelativeError &
890 : , getMergerDelayTimePDF &
891 : , getStarFormationRateDensity &
892 : ) result(binaryMergerRate)
893 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
894 : !DEC$ ATTRIBUTES DLLEXPORT :: getBinaryMergerRate
895 : #endif
896 : use Cosmology_mod, only: LS2HC, OMEGA_DM, OMEGA_DE, getLogLumDisWicMpc
897 : use Constants_mod, only: RK, PI
898 : implicit none
899 : real(RK) , intent(in) :: zplus1
900 : real(RK) , intent(in), optional :: zplus1Max, maxRelativeError
901 : integer(IK) , intent(in), optional :: nRefinement
902 : procedure(getRateDensity_proc) :: getStarFormationRateDensity
903 : procedure(getMergerDelayTimePDF_proc) :: getMergerDelayTimePDF
904 : real(RK) :: binaryMergerRate
905 : real(RK), parameter :: LOG_COEF = log(4._RK*PI*LS2HC)
906 :
907 : binaryMergerRate = exp( LOG_COEF + 2_IK*getLogLumDisWicMpc(zplus1) - ( 3._RK*log(zplus1) + 0.5_RK*log(OMEGA_DM*zplus1**3+OMEGA_DE) ) ) &
908 : * getBinaryMergerRateDensity( zplus1 &
909 : , zplus1Max &
910 : , nRefinement &
911 : , maxRelativeError &
912 : , getMergerDelayTimePDF &
913 : , getStarFormationRateDensity &
914 : )
915 :
916 : end function getBinaryMergerRate
917 :
918 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
919 :
920 : end module StarFormation_mod
|
